Egyptian fractions [Le Monde puzzle #922]

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