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