Egyptian fractions [Le Monde puzzle #922]
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
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.