Greedy algorithm for Egyptian fractions: Difference between revisions
(→Tcl: Added implementation) |
m (Note about another 8-term find) |
||
Line 301: | Line 301: | ||
8/97 has maximum denominator = 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 |
||
</pre> |
</pre> |
||
Note also that <math>\tfrac{44}{53}</math> also has 8 terms. |
|||
:<math>\tfrac{1}{2} + \tfrac{1}{4} + \tfrac{1}{13} + \tfrac{1}{307} + \tfrac{1}{120871} + \tfrac{1}{20453597227} + \tfrac{1}{697249399186783218655} + \tfrac{1}{1458470173998990524806872692984177836808420}</math> |
Revision as of 10:55, 17 April 2014
An Egyptian fraction is the sum of distinct unit fractions such as: . Each fraction in the expression has a numerator equal to 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 to be represented by repeatedly performing the replacement
(simplifying the 2nd term in this replacement as necessary, and where is the ceiling function).
Proper and improper fractions must be able to be expressed. Proper fractions are of the form where and are positive integers, such that , and improper fractions are of the form where and 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: and and
- for all proper fractions, where and 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:
<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 also has 8 terms.