```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 programverifies 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 wouldbe 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
```