r3769
Posts:
30
Registered:
12/12/04


Re: unitary (Egyptian) fractions
Posted:
Mar 27, 2000 12:25 AM


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 >sixterm 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 mipsyear. 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 sixterm (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 I1 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*X0D0 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*X1D1 200 if X2=0 then 370 ' y=p/a0+1/a1 210 D2=A1*D1 220 A2=ceil(D2//X2) 230 X3=A2*X2D2 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*X3D3 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*X4D4 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*X5D5 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

