Greedy algorithm for Egyptian fractions

Greedy algorithm for Egyptian fractions is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

definition

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 (unity) 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

  Proper fractions   are of the form   a/b   where   a   and   b   are positive integers, such that   a < b.

Improper fractions are of the form   a/b   where   a   and   b   are positive integers, such that   a ≥ b.

For improper fractions, the integer part of any improper fraction should be first isolated and shown preceding the Egyptian unit fractions, and be surrounded by square brackets [n].

task requirements

• show the Egyptian fractions for:   ${\tfrac {43}{48}}$   and   ${\tfrac {5}{121}}$   and   ${\tfrac {2014}{59}}$
• for all proper fractions,   a/b   where   a   and   b   are positive one-or two-digit (decimal) integers, find and show an Egyptian fraction that has:

• the largest number of terms,  
• the largest denominator.

• for all one-, two-, and three-digit integers (extra credit), find and show (as above).

Also see

• Wolfram MathWorld^TM entry: Egyptian fraction

Perl 6

<lang perl6>role Egyptian {

   method gist {

join ' + ', (self.abs >= 1 ?? "[{self.floor}]" !! Nil), map {"1/$_"}, self.denominators;

   }
   method denominators {

my ($x, $y) = self.nude; $x %= $y; gather ($x, $y) = (-$y) % $x, $y * take ($y / $x).ceiling while $x;

   }

}

say .nude.join('/'), " = ", $_ but Egyptian for 43/48, 5/121, 2014/59;

my @sample = map { $_ => ($_ but Egyptian).denominators },

   grep * < 1, 
   map {$_ but Egyptian}, 
   ((2 .. 99) X/ (2 .. 99));

say .key.nude.join("/"), " has max denominator, namely ", .value.max given max :by(*.value.max), @sample;

say .key.nude.join("/"), " has max number of denominators, namely ", .value.elems given max :by(*.value.elems), @sample;</lang>

Output:

43/48 = 1/2 + 1/3 + 1/16
5/121 = 1/25 + 1/757 + 1/763309 + 1/873960180913 + 1/1527612795642093418846225
2014/59 = [34] + 1/8 + 1/95 + 1/14947 + 1/670223480
8/97 has max denominator, namely 579504587067542801713103191859918608251030291952195423583529357653899418686342360361798689053273749372615043661810228371898539583862011424993909789665
8/97 has max number of denominators, namely 8

Because the harmonic series diverges (albeit very slowly), it is possible to write even improper fractions as a sum of distinct unit fractions. Here is a code to do that:

<lang perl6>role Egyptian {

   method gist { join ' + ', map {"1/$_"}, self[] }
   method list {

gather { my @h = ([\+] map 1/*, 2 .. *) ...^ * > self; take (2 .. *)[^@h]; my ($x, $y) = (self - (0, @h)[* - 1]).nude; ($x, $y) = (-$y) % $x, $y * take ($y / $x).ceiling while $x; }

   }

}

say 5/4 but Egyptian;</lang>

Output:

1/2 + 1/3 + 1/4 + 1/6

The list of terms grows exponentially with the value of the fraction, though.

Python

(Awaiting task clarification before completion).

<lang python>from fractions import Fraction from math import ceil

class Fr(Fraction):

   def __repr__(self):
       return '%s/%s' % (self.numerator, self.denominator)

def ef(fr):

   ans = []
   if fr >= 1:
       if fr.denominator == 1:
           return [[int(fr)], Fr(0, 1)]
       intfr = int(fr)
       ans, fr = intfr, fr - intfr
   x, y = fr.numerator, fr.denominator
   while x != 1:
       ans.append(Fr(1, ceil(1/fr)))
       fr = Fr(-y % x, y* ceil(1/fr))
       x, y = fr.numerator, fr.denominator
   ans.append(fr)
   return ans

if __name__ == '__main__':

   for fr in [Fr(43, 48), Fr(5, 121), Fr(2014, 59)]:
       print('%r ─► %s' % (fr, ' '.join(str(x) for x in ef(fr))))
   lenmax = denommax = (0, None) 
   for fr in set(Fr(a, b) for a in range(1,100) for b in range(1, 100)):
       e = ef(fr)
       elen, edenom = len(e), e[-1].denominator
       if elen > lenmax[0]:
           lenmax = (elen, fr, e)
       if edenom > denommax[0]:
           denommax = (edenom, fr, e)
   print('Term max is %r with %i terms' % (lenmax[1], lenmax[0]))
   dstr = str(denommax[0])
   print('Denominator max is %r with %i digits %s...%s' %
         (denommax[1], len(dstr), dstr[:5], dstr[-5:]))</lang>

Output:

43/48 ─► 1/2 1/3 1/16
5/121 ─► 1/25 1/757 1/763309 1/873960180913 1/1527612795642093418846225
2014/59 ─► [34] 1/8 1/95 1/14947 1/670223480
Term max is 97/53 with 9 terms
Denominator max is 8/97 with 150 digits 57950...89665

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