The Math Forum

Search All of the Math Forum:

Views expressed in these public forums are not endorsed by NCTM or The Math Forum.

Math Forum » Discussions » sci.math.* » sci.math

Notice: We are no longer accepting new posts, but the forums will continue to be readable.

Topic: unitary (Egyptian) fractions
Replies: 17   Last Post: Mar 31, 2000 6:38 AM

Advanced Search

Back to Topic List Back to Topic List Jump to Tree View Jump to Tree View   Messages: [ Previous | Next ]

Posts: 30
Registered: 12/12/04
Re: unitary (Egyptian) fractions
Posted: Mar 27, 2000 12:25 AM
  Click to see the message monospaced in plain text Plain Text   Click to reply to this topic Reply

Dan Hoey wrote >
>I dusted off some old Egyptian fraction code I wrote a couple of years
>ago. It uses Common Lisp's exact integer arithmetic, so roundoff is
>not an issue. I verified that all fractions up to 732/733 have
>six-term representations and that 732/733 requires seven terms.
>It took about 13 hours on a 296 Mhz Ultra. If a mips is a
>Mhz, that's almost half a mips-year. I don't how much of the
>remaining speed difference is due to the programming language/system
>and how much is due to the algorithm.

It could be the algorithm. For example, the following UBASIC program
verifies all a/b [2<=b<=733, 1<=a<b, (a,b)=1] except
{429/431,562/563,640/643,631/647,653/659,732/733} have a six-term (or less)
expansion in less than 4 MINUTES on my PC. A civilized implementation would
be even faster.


R. Burge

10 point 400
20 clr time
30 for I=2 to 733
40 for J=1 to I-1
50 if gcd(I,J)<>1 then 410
60 Y=J//I
70 P=prmdiv(J)
80 K=0
90 Cnt=0
100 repeat
110 X0=Y//P
120 D0=1
130 A0=ceil(D0//X0)+K
140 if A0@P<>0 then 380
150 X1=A0*X0-D0
160 if X1=0 then 370 ' y=p/a0=1/something (by 150)
170 D1=A0\P
180 A1=ceil(D1//X1)
190 X2=A1*X1-D1
200 if X2=0 then 370 ' y=p/a0+1/a1
210 D2=A1*D1
220 A2=ceil(D2//X2)
230 X3=A2*X2-D2
240 if X3=0 then 370 ' y=p/a0+1/a1+1/a2
250 D3=A2*D2
260 A3=ceil(D3//X3)
270 X4=A3*X3-D3
280 if X4=0 then 370 ' y=p/a0+1/a1+...+1/a3
290 D4=A3*D3
300 A4=ceil(D4//X4)
310 X5=A4*X4-D4
320 if X5=0 then 370 ' y=p/a0+1/a1+...1/a4
330 D5=A4*D4
340 A5=ceil(D5//X5)
350 X6=A5*X5-D5
360 if X6<>0 then 380 ' y=p/a0+1/a1+...1/a5
370 Cnt=1
380 K=K+1
390 until K>100000 or Cnt>0
400 if Cnt=0 then print Y; ' there are some cases where the algo fails
410 next J
420 next I
430 print time

Point your RSS reader here for a feed of the latest messages in this topic.

[Privacy Policy] [Terms of Use]

© The Math Forum at NCTM 1994-2018. All Rights Reserved.