Greedy algorithm for Egyptian fractions

An Egyptian fraction is the sum of distinct unit fractions such as: ${\tfrac {1}{2}}+{\tfrac {1}{3}}+{\tfrac {1}{16}}\,(={\tfrac {43}{48}})$ . Each fraction in the expression has a numerator equal to $1$ and a denominator that is a positive integer, and all the denominators are distinct (i.e., no repetitions).

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)\!\!\!\!\mod x}{y\lceil y/x\rceil }}

(simplifying the 2^nd term in this replacement as necessary, and where $\lceil x\rceil$ is the ceiling function).

Proper and improper fractions must be able to be expressed. Proper fractions are of the form ${\tfrac {a}{b}}$ where $a$ and $b$ are positive integers, such that $a<b$ , and improper fractions are of the form ${\tfrac {a}{b}}$ where $a$ and $b$ are positive integers, such that a ≥ b. (See the REXX programming example to view one method of expressing the whole number part of an improper fraction.)

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, ${\tfrac {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™ 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 { $_ => .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

<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

Tcl

<lang tcl># Just compute the denominator terms, as the numerators are always 1 proc egyptian {num denom} {

   set result {}
   while {$num} {

# Compute ceil($denom/$num) without floating point inaccuracy set term [expr {$denom / $num + ($denom/$num*$num < $denom)}] lappend result $term set num [expr {-$denom % $num}] set denom [expr {$denom * $term}]

   }
   return $result

}</lang> Demonstrating:

Works with: Tcl version 8.6

<lang tcl>package require Tcl 8.6

proc efrac {fraction} {

   scan $fraction "%d/%d" x y
   set prefix ""
   if {$x > $y} {

set whole [expr {$x / $y}] set x [expr {$x - $whole*$y}] set prefix "\[$whole\] + "

   }
   return $prefix[join [lmap y [egyptian $x $y] {format "1/%lld" $y}] " + "]

}

foreach f {43/48 5/121 2014/59} {

   puts "$f = [efrac $f]"

} set maxt 0 set maxtf {} set maxd 0 set maxdf {} for {set d 1} {$d < 100} {incr d} {

   for {set n 1} {$n < $d} {incr n} {

set e [egyptian $n $d] if {[llength $e] >= $maxt} { set maxt [llength $e] set maxtf $n/$d } if {[lindex $e end] > $maxd} { set maxd [lindex $e end] set maxdf $n/$d }

} puts "$maxtf has maximum number of terms = [efrac $maxtf]" puts "$maxdf has maximum denominator = [efrac $maxdf]"</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 maximum number of terms = 1/13 + 1/181 + 1/38041 + 1/1736503177 + 1/3769304102927363485 + 1/18943537893793408504192074528154430149 + 1/538286441900380211365817285104907086347439746130226973253778132494225813153 + 1/579504587067542801713103191859918608251030291952195423583529357653899418686342360361798689053273749372615043661810228371898539583862011424993909789665
8/97 has maximum denominator = 1/13 + 1/181 + 1/38041 + 1/1736503177 + 1/3769304102927363485 + 1/18943537893793408504192074528154430149 + 1/538286441900380211365817285104907086347439746130226973253778132494225813153 + 1/579504587067542801713103191859918608251030291952195423583529357653899418686342360361798689053273749372615043661810228371898539583862011424993909789665

Note also that ${\tfrac {44}{53}}$ also has 8 terms.

{\tfrac {1}{2}}+{\tfrac {1}{4}}+{\tfrac {1}{13}}+{\tfrac {1}{307}}+{\tfrac {1}{120871}}+{\tfrac {1}{20453597227}}+{\tfrac {1}{697249399186783218655}}+{\tfrac {1}{1458470173998990524806872692984177836808420}}