Greedy algorithm for Egyptian fractions

definition

From Wikipedia, an Egyptian fraction is the sum of distinct unit fractions such as: ${\tfrac {1}{2}}+{\tfrac {1}{3}}+{\tfrac {1}{16}}$ .

That is, each fraction in the expression has a numerator equal to 1 and a denominator that is a positive integer, and all the denominators differ from each other (that is, none are repeated).

The (above) Egyptian fraction sums to ${\tfrac {43}{48}}$ .

Fibonacci's Greedy algorithm for Egyptian fractions expands the fraction ${\tfrac {x}{y}}$ to be represented by repeatedly performing the replacement

{\frac {x}{y}}={\frac {1}{\lceil y/x\rceil }}+{\frac {(-y)modx}{y\lceil y/x\rceil }}

(simplifying the 2^nd term in this replacement as necessary, and

where

\lceil x\rceil

is the ceiling function.

This Rosetta Code task will be using the Fibonacci greedy algorithm for computing Egyptian fractions; however, if different method is used instead, please note which method is being employed. Having all the programming examples use the Fibonacci greedy algorithm will make for easier comparison and compatible results.

Proper and improper fractions must be able to be expressed. (See the REXX programming example to view one method of expressing the whole number part of an improper fraction.)

Improper fractions are also known as vulgar or top-heavy fractions. See the Wikipedia article: simple, common, vulgar fractions

task requirements

show the Egyptian fraction for: ${\tfrac {43}{48}}$ and ${\tfrac {5}{121}}$ and ${\tfrac {2014}{59}}$
for all 1- and 2-digit integers, find and show an Egyptian fraction that has:

the largest number of terms
the largest denominator

for all 1-, 2-, and 3-digit integers (extra credit), find and show (as above).

Also see

Wolfram Mathworld^TM entry: Egyptian fraction

REXX

<lang rexx>/*REXX pgm converts a fraction (can be improper) to an Egyptian fraction*/ parse arg fract; z=$egyptF(fract) /*compute the Egyptian fraction. */ say fract ' ───► ' z /*show Egyptian fraction from CL.*/ return z /*stick a fork in it, we're done.*/ /*────────────────────────────────$EGYPTF subroutine────────────────────*/ $egyptF: parse arg z 1 zn '/' zd,,$; if zd== then zd=1 /*whole #?*/ if z= then call erx "no fraction was specified." if zd==0 then call erx "denominator can't be zero:" zd if zn==0 then call erx "numerator can't be zero:" zn if zd<0 | zn<0 then call erx "fraction can't be negative" z if \datatype(zn,'W') then call erx "numerator must be an integer:" zn if \datatype(zd,'W') then call erx "denominator must be an integer:" zd _=zn%zd /*check if it's an improper fract*/ if _>=1 then do /*if improper fract, then append.*/

             $='['_"]"                /*append the whole# part of fract*/
             zn=zn-_*zd               /*now, just use the proper fract.*/
             if zn==0  then return $  /*if no fraction, we're done.    */
             end

if zd//zn==0 then do; zd=zd%zn; zn=1; end

 do  forever
 if zn==1 & datatype(zd,'W')  then return $ '1/'zd   /*append E. fract.*/
 nd=zd%zn+1;      $=$ '1/'nd          /*add unity to int fract., append*/
 z=$fractSub(zn'/'zd,  "-",  1'/'nd)  /*go and subtract the two fracts.*/
 parse var z zn '/' zd                /*extract the numerator & denom. */
 L=2*max(length(zn),length(zd))       /*calculate if we need more digs.*/
 if L>=digits()  then numeric digits L*2   /*yes, then bump the digits.*/
 end   /*forever*/                    /* [↑]  loop ends when  zn==1.   */

/*────────────────────────────────$FRACTSUB subroutine──────────────────*/ $fractSub: procedure; parse arg z.1,,z.2 1 zz.2; arg ,op

                        do j=1  for 2;   z.j=translate(z.j,'/',"_");  end

if z.1== then z.1=(op\=="+" & op\=='-') /*unary +,- first fract.*/ if z.2== then z.2=(op\=="+" & op\=='-') /*unary +.- second fract.*/

 do j=1  for 2                             /*process both fractions.   */
 if pos('/',z.j)==0     then z.j=z.j"/1";     parse var  z.j  n.j '/' d.j
 if \datatype(n.j,'N')  then call erx "numerator isn't an integer:"   n.j
 if \datatype(d.j,'N')  then call erx "denominator isn't an integer:" d.j
 n.j=n.j/1;   d.j=d.j/1                    /*normalize numerator/denom.*/

      do  while \datatype(n.j,'W');   n.j=n.j*10/1;   d.j=d.j*10/1;   end
                                           /* [↑] normalize both nums. */
 if d.j=0  then call erx "denominator can't be zero:"  z.j
 g=gcd(n.j,d.j);   if g=0  then iterate;     n.j=n.j/g;         d.j=d.j/g
 end    /*j*/

l=lcm(d.1 d.2); do j=1 for 2; n.j=l*n.j/d.j; d.j=l; end if op=='-' then n.2=-n.2 t=n.1+n.2; u=l if t==0 then return 0; g=gcd(t,u); t=t/g; u=u/g if u==1 then return t

             return t'/'u

/*═════════════════════════════general 1-line subs══════════════════════*/ erx: say; say '***error!***' arg(1); say; exit 13 gcd:procedure;$=;do i=1 for arg();$=$ arg(i);end;parse var $ x z .;if x=0 then x=z;x=abs(x);do j=2 to words($);y=abs(word($,j));if y=0 then iterate;do until _==0;_=x//y;x=y;y=_;end;end;return x lcm:procedure;y=;do j=1 for arg();y=y arg(j);end;x=word(y,1);do k=2 to words(y);!=abs(word(y,k));if !=0 then return 0;x=x*!/gcd(x,!);end;return x p: return word(arg(1),1)</lang>

output when the input used is: 43/48


43/48  ───►   1/2 1/3 1/16

output when the input used is: 5/121

5/121  ───►   1/25 1/757 1/763309 1/873960180913 1/1527612795642093418846225

output when the input used is: 2014/59

2014/59  ───►   [34] 1/8 1/95 1/14947 1/670223480

The following is a driver program to address the requirements to find the largest number of terms for a 1- or 2-digit integer, and the largest denominator. Also, the same program is used for the 1-, 2-, and 3-digit extra credit task. <lang rexx> /*REXX pgm runs the EGYPTIAN program to find bigest denominator & #terms*/ parse arg top . /*get optional parameter from CL.*/ if top== then top=99 /*Not specified? Then use default*/ bigD=; bigT=; maxT=0; maxD=0 /*initialize some REXX variables.*/

                                      /* [↓]  determine biggest,longest*/
     do n=2      to top               /*traipse through the numerators.*/
         do d=n+1  to top             /*   "       "     " denominators*/
         fract=n'/'d                  /*create the fraction to be used.*/
         y='EGYPTIAN'(fract)          /*invoke the other REXX program. */
         t=words(y)                   /*find out how many terms in E.F.*/
         if t>maxT  then bigT=fract   /*is this a new high for # terms?*/
         maxT=max(maxT,T)             /*find the maximum number terms. */
         b=substr(word(y,t),3)        /*get the denominator from the EF*/
         if b>maxD  then bigD=fract   /*is this a new denominator high?*/
         maxD=max(maxD,b)             /*find the maximum denominator.  */
         end   /*d*/                  /* [↑]  only use proper fractions*/
     end       /*n*/                  /* [↑]  ignore  1/n  fractions.  */
                                      /* [↑]  display longest, biggest.*/

@='in the Egyptian fractions used is' /*literal to make a shorter SAY.*/ say 'largest number of terms' @ maxT "terms for" bigT say 'highest denominator' @ length(maxD) "digits is for" bigD':' say maxD /*stick a fork in it, we're done.*/</lang> output for all 1- and 2-digit integers (when using the default input):

largest number of terms in the Egyptian fractions used is 8 terms for 8/97
largest denominator in the Egyptian fractions is 150 digits is for 8/97
579504587067542801713103191859918608251030291952195423583529357653899418686342360361798689053273749372615043661810228371898539583862011424993909789665

output for all 1-, 2-, and 3-digit integers (when using for input: -999

largest number of terms in the Egyptian fractions used is 13 terms for 529/914
largest denominator in the Egyptian fractions is 2847 digits is for 36/457