Date: Mar 27, 2000 12:25 AM
Author: r3769
Subject: Re: unitary (Egyptian) fractions




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.

Regards,

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