Archive for Liber Abaci

puzzles & riddles

Posted in Books, Kids, R, Statistics with tags , , , , , , , , on January 3, 2021 by xi'an

A rather simplistic game on the Riddler of 18 December:

…two players, each of whom starts with a whole number of points. Players take turns “attacking” each other, which involves subtracting their own number of points from their opponent’s until one of the players is out of points.

Easy to code in R:

g=function(N,M)ifelse(N<M,-g(M-N,N),1)

f=function(N,M=N){
  while(g(N,M)>0)M=M+1
  return(M)}

which converges to the separating ratio 1.618. If decomposing the actions until one player wins, one gets a sequence of upper and lower bounds associated with the Fibonacci sequence: 1⁻, 2⁺, 3/2⁻, 5/3⁺, 8/5⁻, &tc, converging to the “golden ratio” φ.

As an aside, I also solved a relatively quick codegolf challenge, where the question was to find the sum of all possible binary values from a bitmask. Meaning that for a binary input, e.g., 101X0XX0…01X, with some entries masked by X’s, one had to find the sum of all binary numbers compatible with the input. Which can be solved succinctly by counting the number of X’s, k, and adding the visible bits 2^k times and replacing the invisible ones by  2^{k-1}. With some help, and using 2 instead of X, my R code moved from 158 bytes to 50:

function(x)2^sum(a<-x>1)*rev(x/4^a)%*%2^(seq(x)-1)

riddles on Egyptian fractions and Bernoulli factories

Posted in Books, Kids, R with tags , , , , , , , , , , , , , on June 11, 2019 by xi'an

Two fairy different riddles on the weekend Riddler. The first one is (in fine) about Egyptian fractions: I understand the first one as

Find the Egyptian fraction decomposition of 2 into 11 distinct unit fractions that maximises the smallest fraction.

And which I cannot solve despite perusing this amazing webpage on Egyptian fractions and making some attempts at brute force  random exploration. Using Fibonacci’s greedy algorithm. I managed to find such decompositions

2 = 1 +1/2 +1/6 +1/12 +1/16 +1/20 +1/24 +1/30 +1/42 +1/48 +1/56

after seeing in this short note

2 = 1 +1/3 +1/5 +1/7 +1/9 +1/42 +1/15 +1/18 +1/30 +1/45 +1/90

And then Robin came with the following:

2 = 1 +1/4 +1/5 +1/8 +1/10 +1/12 +1/15 +1/20 +1/21 +1/24 +1/28

which may prove to be the winner! But there is even better:

2 = 1 +1/5 +1/6 +1/8 +1/9 +1/10 +1/12 +1/15 +1/18 +1/20 +1/24

The second riddle is a more straightforward Bernoulli factory problem:

Given a coin with a free-to-choose probability p of head, design an experiment with a fixed number k of draws that returns three outcomes with equal probabilities.

For which I tried a brute-force search of all possible 3-partitions of the 2-to-the-k events for a range of values of p from .01 to .5 and for k equal to 3,4,… But never getting an exact balance between the three groups. Reading later the solution on the Riddler, I saw that there was an exact solution for 4 draws when

p=\frac{3-\sqrt{3(4\sqrt{9}-6)}}{6}

Augmenting the precision of my solver (by multiplying all terms by 100), I indeed found a difference of

> solver((3-sqrt(3*(4*sqrt(6)-9)))/6,ba=1e5)[1]
[1] 8.940697e-08

which means an error of 9 x 100⁻⁴ x 10⁻⁸, ie roughly 10⁻¹⁵.

Egyptian fractions [Le Monde puzzle #922]

Posted in Books, Kids, R with tags , , , , , , , on July 28, 2015 by xi'an

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