Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

For its summer edition, Le Monde mathematical puzzle switched to a lighter version with immediate solution. This #922 considers Egyptian fractions which only have distinct denominators (meaning the numerator is always 1) and can be summed. This means 3/4 is represented as ½+¼. Each denominator only appears once. As I discovered when looking on line, a lot of people are fascinated with this representation and have devised different algorithms to achieve decompositions with various properties. Including Fibonacci who devised a specific algorithm called the greedy algorithm in 1202 in the Liber Abaci. In the current Le Monde edition, the questions were somewhat modest and dealt with the smallest decompositions of 2/5, 5/12, and 50/77 under some additional constraint.

Since the issue was covered in so many places, I just spent one hour or so constructing a basic solution à la Fibonacci and then tried to improve it against a length criterion. Here are my R codes (using the numbers library):

osiris=function(a,b){
#can the fraction a/b be simplified
diva=primeFactors(a)
divb=primeFactors(b)
divc=c(unique(diva),unique(divb))
while (sum(duplicated(divc))>0){
n=divc[duplicated(divc)]
for (i in n){a=div(a,i);b=div(b,i)}
diva=primeFactors(a)
divb=primeFactors(b)
divc=c(unique(diva),unique(divb))
}
return(list(a=a,b=b))
}


presumably superfluous for simplifying fractions

horus=function(a,b,teth=NULL){
#simplification
anubis=osiris(a,b)
a=anubis$a;b=anubis$b
#decomposition by removing 1/b
isis=NULL
if (!(b %in% teth)){
a=a-1
isis=c(isis,b)
teth=c(teth,b)}
if (a>0){
#simplification
anubis=osiris(a,b)
bet=b;a=anubis$a;b=anubis$b
if (bet>b){ isis=c(isis,horus(a,b,teth))}else{
# find largest integer
k=ceiling(b/a)
while (k %in% teth) k=k+1
a=k*a-b
b=k*b
isis=c(isis,k,horus(a,b,teth=c(teth,k)))
}}
return(isis)}


which produces a Fibonacci solution (with the additional inclusion of the original denominator) and

nut=20
seth=function(a,b,isis=NULL){
#simplification
anubis=osiris(a,b)
a=anubis$a;b=anubis$b
if ((a==1)&(!(b %in% isis))){isis=c(isis,b)}else{
ra=hapy=ceiling(b/a)
if (max(a,b)<1e5) hapy=horus(a,b,teth=isis)
k=unique(c(hapy,ceiling(ra/runif(nut,min=.1,max=1))))
propa=propb=propc=propd=rep(NaN,le=length((k %in% isis)))
bastet=1
for (i in k[!(k %in% isis)]){
propa[bastet]=i*a-b
propb[bastet]=i*b
propc[bastet]=i
propd[bastet]=length(horus(i*a-b,i*b,teth=c(isis,i)))
bastet=bastet+1
}
k=propc[order(propd)[1]]
isis=seth(k*a-b,k*b,isis=c(isis,k))
}
return(isis)}


which compares solutions against their lengths. When calling those functions for the three fractions above the solutions are

> seth(2,5)
[1] 15 3
> seth(5,12)
[1] 12  3
> seth(50,77)
[1]   2 154   7


with no pretension whatsoever to return anything optimal (and with some like crashes when the magnitude of the entries grows, try for instance 5/121). For this latest counter-example, the alternative horus works quite superbly:

> horus(5,121)
[1] 121 31 3751 1876 7036876


Filed under: Books, Kids, R Tagged: Egyptian fractions, Fibonacci, greedy algorithm, Le Monde, Liber Abaci, mathematical puzzle, numerics, Rhind papyrus