Egyptian Fractions on the WP-34s. Just for the heck of it.



#15

This program does a brute force search for the shortest Egyptian Fraction that is equal to the fraction that you enter. It returns the first solution it finds (there might be many).

This is my first attempt at a meaningful program on the WP-34s. It feels clunky and way too long, but for a first attempt at a 34s program, I was just happy to get it to work.

I suspect Valentin could whip this out in 4 or 5 lines on the 71b.

This operates in 64 bit integer mode. It will also work in 32 bit integer mode (just change the one instruction at the beginning), but it will overflow much easier. It currently makes no attempt to detect overflows - something I really should have put in.


		##########################################################################	
# Calculate Egyptian Fraction
# Invoke with:
# Numerator in Y
# Denominator in X
#
# 1) Only works if numerator < denominator (i.e. Y/X < 1)
# 2) Only works on positive numbers
# 3) Does not check for overflow (which can happen)
# 4) Runs in 64 bit integer mode; restores your original mode when done.
# 5) Sets stack size to 4. Current code will not work with stack size 8.
#
# When it has found a solution, the program will stop and display
# the denominator of the first fraction. Hit R/S to display the next
# denominator, and so on. It will display 0 after the last fraction.
#
##########################################################################
1 LBL 'EGF'
2 LocR 004
3 STOM .03
4 BASE 10 ; Calculations depend on integer math
5 WSIZE 64 ; Word size limits the range
6 SSIZE4 ; Code needs modification for SSIZE8
7 CF 00
8 STO .01 ; Denominator stored in .01
9 R DOWN
10 STO .00 ; Numerator stored in .00
11 1
12 STO 20 ; Holds the search length (number of terms)
13 LBL 01 ; Start with one term and hunt for shortest sequence
14 1
15 0
16 X<=? 20 ; This code only supports a max of 10 terms
17 GTO 02
18 RCL .00
19 RCL .01
20 GCD
21 STO .02
22 0
23 RCL .01
24 RCL .02
25 /
26 RCL .00
27 RCL .02
28 /
29 2
30 R DOWN
31 XEQ 03 ; Call the recursion that does the search
32 FS? 00 ; If the flag is set, we found a sequence
33 GTO 02
34 1
35 STO+ 20
36 GTO 01
37 LBL 02
38 1
39 STO .00
40 LBL 06 ; Ouput the denominators 1 by 1
41 0
42 X=? 20
43 GTO 07
44 RCL->.00
45 STOP
46 2
47 STO+ .00
48 1
49 STO- 20
50 GTO 06
51 LBL 07
52 RCLM .03
53 0
54 RTN
##########################################################################
# Recursively search for a sequence of terms that works
##########################################################################
55 LBL 03
56 LocR 011
57 STOS .00
58 RCL .02
59 X<>? 20 ; If we are the right number of terms
60 GTO 04 ; We have either found something or need to leave
61 0
62 X<>? .00 ; If there is a remainder
63 RTN ; We didn't find anything
64 SF 00 ; We have found something
65 RTN
66 LBL 04 ; Continue hunting
67 0
68 X=? .00
69 RTN
70 RCL .01
71 RCL .00
72 +
73 1
74 -
75 RCL .00
76 /
77 STO .04
78 STO .05
79 x>=? .03
80 GTO 05
81 RCL .03
82 STO .05
83 LBL 05
84 RCL 20
85 RCL .02
86 -
87 RCL .04
88 *
89 X<? .05
90 RTN
91 RCL .02
92 2
93 *
94 STO .06
95 1
96 +
97 STO .07
98 1
99 STO->.06
100 RCL .05
101 STO->.07
102 RCL .00
103 RCL->.07
104 *
105 RCL->.06
106 RCL .01
107 *
108 -
109 STO .08
110 RCL .01
111 RCL->.07
112 *
113 STO .09
114 GCD
115 STO .10
116 1
117 STO+ .05
118 RCL .02
119 1
120 +
121 RCL .09
122 RCL .10
123 /
124 RCL .08
125 RCL .10
126 /
127 RCL .05
128 R DOWN
129 XEQ 03
130 FS? 00
131 RTN
132 GTO 05
##########################################################################
133 END


#16

Thanks Marcel! For those (like me) who wonder what's an Egyptian fraction at all, however: See here a short explanation in German or a more elaborated one in English first.

d:-)

#17

  • Could you elaborate a little the brute force algorithm that is used?
  • What's the result for 5/121?
  • Unfortunately the comments are ending here:
    Quote:
    Continue hunting

    That's where the interesting stuff begins.

Kind regards

Thomas


#18

5/121 = 1/33 + 1/121 + 1/363 ?

#19

The Wikipedia reference Walter gave provides some heuristics for simplifying into Egyptian fractions. I don't know if these are complete in of themselves, but I'm sure some number theorist has come up with a good selection. The 34S can test for primality quickly, the integer factorisation function never made it in for pace reasons, although it was coded as an internal function.

- Pauli

#20

I knew the rest of my weekend was going to be very busy with family obligations and I posted what I had. In retrospect I guess I should have waited and commented first.

The wp-34s program was derived from this source code:

http://www.ics.uci.edu/~eppstein/numth/egypt/efrac.c

I was not able to find the original paper by the author that this code references, but the code is easier to read than the wp-34s program.

I will repost a commented version.

The answer it gives to 5/121 is:

1/25 + 1/759 + 1/208725

As I mentioned, it returns the first solution it finds. My initial plan was to find all the solutions of the shortest length and return the one that minimized the sum of the denominators. However, that might have been a bit slow.

#21

Below is an updated version of the program. The code is a bit tighter and has better commenting.

A summary of the brute force approach is as follows:

There is a recursive routine that does the search for the terms, but it requires knowledge of how many terms there will be.

In order to make proper use of the recursive routine (LBL 03), we have a main routine which first searches for solutions that are 1 term long, 2 terms long, 3 terms long, etc...

The recursion is started with the value we are searching for. On entry into the recursion, it takes an educated guess at the next term, calculates the delta between the solution and that term, and then recursively calls itself with that delta as the solution being searched for.

If the recursion encounters a delta of 0, then it has found a solution.

If it does not find a solution, then it increments the test denominator by 1 and tries again.

The recursion is limited by the fact that it searches the denominators in increasing order with each term. Thus, for a given delta, if there are N terms remaining to find, the next term cannot be less than delta/N. Once the test denominator gets large enough, the recursion can be stopped because it will not be possible to find an answer.

Thus it converges on a solution relatively quickly.

		##########################################################################	
# Calculate Egyptian Fraction
# Invoke with:
# Numerator in Y
# Denominator in X
#
# 1) Only works if numerator < denominator (i.e. Y/X < 1)
# 2) Only works on positive numbers
# 3) Does not check for overflow (which can happen)
# 4) Runs in 64 bit integer mode; restores your original mode when done.
# 5) Sets stack size to 4. Current code will not work with stack size 8.
# 6) Maximum number of terms is hardwired to 20.
#
# When it has found a solution, the program will stop and display
# the denominator of the first fraction. Hit R/S to display the next
# denominator, and so on. It will display 0 after the last fraction.
#
# Registers 00 to 19: the solution terms that have been found
# Register 20 Current length of solution being searched.
#
##########################################################################
1 LBL 'EGF'
2 LocR 004
3 STOM .03 ; Save current mode so we can restore it on completion
4 BASE 10 ; Calculations depend on integer math
5 WSIZE 64 ; Overflows easily with 32 bit words
6 SSIZE4 ; Code needs minor modification for SSIZE8
7 CF 00 ; Flag 00 indicates a solution was found (exit from recursion)
8 STO .01 ; Denominator stored in .01
9 R DOWN
10 STO .00 ; Numerator stored in .00
11 1
12 STO 20 ; Start searching for solutions of length 1
13 LBL 01 ; Loop with R20 incrementing the number of terms we hunting for
14 2
15 0
16 X<=? 20 ; This code only supports a max of 20 terms
17 GTO 07
18 RCL .00
19 RCL .01
20 GCD
21 STO .02 ; GCD of input numerator & denominator
22 0 ; 3rd param for XEQ 03 – which term we are looking for
23 RCL .01
24 RCL .02
25 / ; 2nd param for XEQ 03 – solution (input) denominator/GCD
26 RCL .00
27 RCL .02
28 / ; 1st param for XEQ 03 – solution (input) numerator/GCD
29 2 ; 4th param for XEQ 03 – denominator to start searching for next term
30 R DOWN
31 XEQ 03 ; Call the recursion that does the search
32 FS? 00 ; If the flag is set, we found a sequence
33 GTO 02
34 1
35 STO+ 20 ; Didn't find solution, increment the search length
36 GTO 01 ; and try again

37 LBL 02 ; Prepare too loop outputing the result
38 0
39 STO .00
40 LBL 06 ; Ouput the denominators 1 by 1
41 0
42 X=? 20
43 GTO 07
44 RCL->.00
45 STOP
46 1
47 STO+ .00
48 STO- 20
49 GTO 06
50 LBL 07
51 RCLM .03
52 0
53 RTN
##########################################################################
# Recursively search for a sequence of terms that works
# X: Numerator of solution to search for
# Y: Denominator of solution to search for
# Z: Term index in solution we are hunting for (0 = first term)
# T: minimum denominator for next term (since we don't allow duplicates and search in increasing order)
#
# Try the next possible denominator for the next term. Computes the error and recursively
# call ourselves with the error as the next solution we are looking for. If the error is 0
# we have found a solution, so set a flag and pop out of the recursion.
#
# R.00 (input) Numerator of remaining solution to search for
# R.01 (input) Denominator of remaining solution to search for
# R.02 (input) Index of term we are modifying
# R.03 (input) Minimum denominator for next term
# R.04 (computed) max bound for denominator for next term
# R.05 (computed) the denominator we are trying for the next term
# R.06 GCD of next solution (i.e. the remainder) we will recursively search for
# R.07 Denominator of next solution (i.e. the remainder) we will recursively search for
# R.08 Numerator of next solution (i.e. the remainder) we will recursively search for
##########################################################################
54 LBL 03
55 LocR 009
56 STOS .00
57 RCL .02
58 X<>? 20 ; If we are the right number of terms
59 GTO 04 ; We have either found something or need to leave
60 0
61 X<>? .00 ; If there is a remainder
62 RTN ; We didn't find anything
63 SF 00 ; We have found something
64 RTN
65 LBL 04 ; Continue hunting
66 RCL .01 ; Computes an upper bound for the next denominator
67 RCL .00 ; Works because we know how many terms we have left
68 +
69 1
70 -
71 RCL .00
72 /
73 STO .04
74 STO .05 ; start trying denominators with the computed bound
75 x>=? .03 ; but override it with the input parameter if that was greater.
76 GTO 05
77 RCL .03
78 STO .05
79 LBL 05 ; Loop, incrementing the test denominator
80 RCL 20 ; Test to see if the bounds have been exceeded
81 RCL .02
82 -
83 RCL .04
84 *
85 X<? .05
86 RTN ; pop out of the recursion if the bounds were exceeded
87 RCL .05 ; R.05 contains the next denominator to try
88 STO->.02 ; enter it as the next term in the solution
89 RCL .00 ; This section computes the delta
90 * ; between the solution we are searching for
91 RCL .01 ; and the term we just enetered into the solution.
92 - ; That will become the next solution we
93 STO .08 ; search for
94 RCL .01
95 RCL .05
96 *
97 STO .07
98 GCD ; Compute GCD so we can reduce the fraction
99 STO .06
100 INC .05 ; Next denominator must be at least 1 larger
101 RCL .02
102 1
103 + ; 3rd param for XEQ 03 – which term we are looking for
104 RCL .07
105 RCL .06
106 / ; 2nd param for XEQ 03 – solution (input) denominator/GCD
107 RCL .08
108 RCL .06
109 / ; 1st param for XEQ 03 – solution (input) numerator/GCD
110 RCL .05 ; 4th param for XEQ 03 – denominator to start searching for next term
111 R DOWN
112 XEQ 03 ; Recursively call ourselves searching for the delta
113 FS? 00 ; If the flag was set in the recursion (we found a solution)...
114 RTN ; .. then don't loop any more; we are done
115 GTO 05
##########################################################################
116 END

Edited: 1 Apr 2013, 2:18 p.m.


#22

Many thanks for your effort! Much appreciated! Now I'm looking forward to read and try to understand your code.

Kind regards

Thomas


#23

I had a few beers while writing that code. It is possible that you may need to be suitably lubricated to read it.

If you do make the effort, please remember that I wrote it as an exercise in trying to familiarize myself with the 34s. If you come across places where the features of the 34s could have been better taken advantage of, I would love to hear it.


#24

A couple of thoughts:

  • Recall arithmetic: e.g. lines 89 and 90 could be replaced by RCL* .00 There are quite a few instances where this would save a step.

  • INC and DEC to add and subtract one without changing LastX. e.g. lines 102 and 103 could be replaced by INC X and lines 46 through 48 by INC .00 DEC 20

  • Local flags. Once you execute the LocR command, flags .00 - .15 become available. These are initialised to clear by the LocR command. You get a separate set of flags for each LocR command so this mightn't be useful in this instance.

  • Lines 14 and 15 can be replaced by # 020 (in the constants catalogue as '#'). Saving a step and more time.

  • Double precision mode (DBLON) uses real numbers with 34 decimal digits. This would allow a greater range than sixty four bit integer arithmetic. The stack is converted automatically on mode change. Failing this, forcing unsigned integer mode would give you another bit (you don't set the integer sign mode at all which is probably a bug).

  • The stack shuffle <-> command doesn't seem relevant here, but it is worthwhile keeping in mind. It lets you rearrange X, Y, Z & T however you desire.

- Pauli


#25

Thanks for the suggestions!

The double precision idea is an interesting one and I will give that a try. I had though of doing a real version so that I could compare the relative speed of real versus integer modes for this particular application. However it never occurred to me that I could actually increase the range.

The stack shuffle command is definitely the one that most stood out in my mind when I first started poking around the WP-34s. I am curious as to whether it is original to this platform or are there other platforms that have it?

One thing that I wanted to do was eliminate any usage of global registers. My inner C programmer wanted the main routine to allocate local registers for the result and then have the recursive routine stuff those at is is going along. However that would require a mechanism by which a subroutine could be given access to the local data of its caller. Did I miss something - is that possible?


#26

Quote:
The stack shuffle command is definitely the one that most stood out in my mind when I first started poking around the WP-34s. I am curious as to whether it is original to this platform or are there other platforms that have it?

AFAIK there's no other calculator featuring a shuffle command.
Quote:
One thing that I wanted to do was eliminate any usage of global registers. My inner C programmer wanted the main routine to allocate local registers for the result and then have the recursive routine stuff those at is is going along. However that would require a mechanism by which a subroutine could be given access to the local data of its caller. Did I miss something - is that possible?

Please see pp. 79f of the printed manual.

d:-)


#27

I ordered my printed manual last week and it is scheduled to arrive tomorrow. I will wait with bated breath to see what page 79f will tell me.

#28

Quote:
However that would require a mechanism by which a subroutine could be given access to the local data of its caller. Did I miss something - is that possible?

Yes, but only if the subroutine doesn't allocate any local registers itself. Or at least before it allocates its own locals.


- Pauli


Possibly Related Threads...
Thread Author Replies Views Last Post
  [WP-34S] Unfortunate key damage with update to V3 :( svisvanatha 5 351 12-10-2013, 11:37 PM
Last Post: Les Bell
  WP-34S (Emulator Program Load/Save) Barry Mead 1 193 12-09-2013, 05:29 PM
Last Post: Marcus von Cube, Germany
  DIY HP 30b WP 34s serial flash/programming cable Richard Wahl 2 260 12-04-2013, 11:14 AM
Last Post: Barry Mead
  WP 34S/43 ?'s Richard Berler 3 285 11-10-2013, 02:27 AM
Last Post: Walter B
  My FrankenCulator (wp-34s) FORTIN Pascal 4 298 11-09-2013, 06:18 PM
Last Post: FORTIN Pascal
  WP 34S Owner's Handbook Walter B 5 439 11-09-2013, 05:34 PM
Last Post: Harald
  wp 34s overlay and programming. FORTIN Pascal 6 369 11-08-2013, 01:28 PM
Last Post: Nick_S
  m.dy in display of WP-34S Harold A Climer 3 210 11-05-2013, 11:28 AM
Last Post: Andrew Nikitin
  WP 34s : An old problem coming back Miguel Toro 2 212 11-05-2013, 03:00 AM
Last Post: Marcus von Cube, Germany
  [WP 34s] Pressure Conversion Factors Timothy Roche 8 429 11-04-2013, 07:17 PM
Last Post: Dave Shaffer (Arizona)

Forum Jump: