Archive for R

Le Monde puzzle [#1164]

Posted in Books, Kids, R with tags , , , , , , , , , , , , , on November 16, 2020 by xi'an

The weekly puzzle from Le Monde is quite similar to older Diophantine episodes (I find myself impossible to point out):

Give the maximum integer that cannot be written as 105x+30y+14z. Same question for 105x+70y+42z+30w.

These are indeed Diophantine equations and the existence of a solution is linked with Bézout’s Lemma. Take the first equation. Since 105 and 30 have a greatest common divisor equal to 3×5=15, there exists a pair (x⁰,y⁰) such that

105 x⁰ + 30 y⁰ = 15

hence a solution to every equation of the form

105 x + 30 y = 15 a

for any relative integer a. Similarly, since 14 and 15 are co-prime,

there exists a pair (a⁰,b⁰) such that

15 a⁰ + 14 b⁰ = 1

hence a solution to every equation of the form

15 a⁰ + 14 b⁰ = c

for every relative integer c. Meaning 105x+30y+14z=c can be solved in all cases. The same result applies to the second equation. Since algorithms for Bézout’s decomposition are readily available, there is little point in writing an R code..! However, the original question must impose the coefficients to be positive, which of course kills the Bézout’s identity argument. Stack Exchange provides the answer as the linear Diophantine problem of Frobenius! While there is no universal solution for three and more base integers, Mathematica enjoys a FrobeniusNumber solver. Producing 271 and 383 as the largest non-representable integers. Also found by my R code

o=function(i,e,x){
  if((a<-sum(!!i))==sum(!!e))sol=(sum(i*e)==x) else{sol=0
    for(j in 0:(x/e[a+1]))sol=max(sol,o(c(i,j),e,x))}
  sol}
a=(min(e)-1)*(max(e)-1)#upper bound
M=b=((l<-length(e)-1)*prod(e))^(1/l)-sum(e)#lower bound
for(x in a:b){sol=0
for(i in 0:(x/e[1]))sol=max(sol,o(i,e,x))
M=max(M,x*!sol)}

(And this led me to recover the earlier ‘Og entry on the coin problem! As of last November.) The published solution does not bring any useful light as to why 383 is the solution, except for demonstrating that 383 is non-representable and any larger integer is representable.

on arithmetic derivations of square roots

Posted in Books, Kids, pictures, R, Statistics with tags , , , , , , , , , on November 13, 2020 by xi'an

An intriguing question made a short-lived appearance on the CodeGolf section of Stack Exchange, before being removed, namely the (most concise possible) coding of an arithmetic derivation of the square root of an integer, S, with a 30 digit precision and using only arithmetic operators. I was not aware of the myriad of solutions available, as demonstrated on the dedicated WIkipedia page. And ended playing with three of them during a sleepless pre-election night!

The first solution for finding √S is based on a continued fraction representation of the root,

\sqrt{S}=a+\cfrac{r}{2a+\cfrac{r}{2a+\ddots}}

with a²≤S and r=S-a². It is straightforward to code-golf:

while((r<-S-T*T)^2>1e-9)T=(F<-2*T+r/(2*T+F))-T;F

but I found it impossible to reach the 30 digit precision (even when decreasing the error bound from 10⁻⁹). Given the strict rules of the game, this would have been considered to be a failure.

The second solution is Goldschmidt’s algorithm

b=S
T=1/sum((1:S)^2<S) 
while((1-S*prod(T)^2)^2>1e-9){
  b=b*T[1]^2
  T=c((3-b)/2,T)}
S*prod(T)

which is longer for code-golfing but produces both √S and 1/√S (and is faster than the Babylonian method and Newton-Raphson). Again no luck with high precision and almost surely unacceptable for the game.

The third solution is the most interesting [imho] as it mimicks long division, working two digits at a time (and connection with Napier’s bones)

`~`=length
D=~S
S=c(S,0*(1:30))
p=d=0
a=1:9
while(~S){ 
  F=c(F,x<-sum(a*(20*p+a)<=(g<-100*d+10*S[1]+S[2])))
  d=g-x*(20*p+x)
  p=x+10*p
  S=S[-1:-2]}
sum(10^{1+D/2-1:~F}*F)

plus providing an arbitrary number of digits with no error. This code requires S to be entered as a sequence of digits (with a possible extra top digit 0 to make the integer D even). Returning one digit at a time, it would further have satisfied the constraints of the question (if in a poorly condensed manner).

Le Monde puzzle [#1158]

Posted in Books, Kids, R with tags , , , , , on November 10, 2020 by xi'an

A weekly puzzle from Le Monde on umbrella sharing:

Four friends, Antsa, Cyprien, Domoina and Fy, are leaving school to return to their common housing. It is raining and they only have one umbrella with only room for two. Given walking times, x¹, x², x³ and x⁴, find the fastest time by which all of the four will be home, assuming they all agree to come back with the umbrella if need be.

A recursive R function produces the solution

bez=function(starz=rexp(4),finiz=rep(0,4),rtrn=F){
  if((!rtrn)&(sum(starz>0)==2)){return(max(starz))
    }else{
      tim=1e6
      if(rtrn){
        for(i in (1:4)[finiz>0]){
          nstart=starz;nstart[i]=finiz[i]
          nfini=finiz;nfini[i]=0
          targ=finiz[i]+bez(nstart,nfini,FALSE)
          if(targ<tim){tim=targ}} 
          }else{
          for(i in (1:4)[starz>0])
          for(j in (1:4)[starz>0]){
            if(i!=j){
              nstar=starz;nstar[i]=nstar[j]=0
              nfini=finiz;nfini[i]=starz[i];nfini[j]=starz[j]
              targ=max(starz[i],starz[j])+bez(nstar,nfini,TRUE)
              if (targ<tim){tim=targ}
            }}}
      return(tim)}

which gives for instance

> bez()
[1] 3.297975
> bez(1:4)
[1] 11
> bez(rep(3,4))
[1] 15

asymmetric information

Posted in Kids, R with tags , , , , , on November 4, 2020 by xi'an

The Riddler of 16 October had the following puzzle:

Take a real number θ uniformly distributed over (0,100). Among three players, the winner is whoever guessed the closest price without going over θ. In the event all guesses exceeded θ, the contestant with the lowest (and therefore closest) guess is declared the winner. The second player knows the first player’s guess and the third player knows both other guesses. What is the optimal guess for the first player, assuming all players maximise their probability of winning?

Looking at the optimal solution z for the third player leads to six possible choices, depending on the connection between the other guesses, x and y. Which translates in the R code

topz=function(x,y){
  if((2*y>=x)&(y>=1-x))  z=y-.001
  if(max(4*y,1+y)<=2*x)  z=y+.001
  if((2*x<=1+y)&(x<=1-y))z=x+.001
  z}
  
third=function(x,y) ifelse(y<x,topz(x,y),topz(y,x))

For there, the optimal choice y for the second player follows and happens on a boundary of one of the six regions, which itself returns that the optimal choice for the first player is x=2/3, leading to equal chances of winning (although there is some uncertainty on the boundaries). It is thus feasible to beat the asymmetric information. The picture above was my attempt at representing the probabilities of gain for all three players, some of the six regions being clearly visible, with first axis being x and second being y [and z is one of x⁻,x⁺,y⁻,y⁺]. The R code is too pedestrian to be reproduced!

sampling w/o replacement except when replacing

Posted in Books, Kids, R with tags , , , , , , , on November 3, 2020 by xi'an

Another Riddle(r), considering a box with M myrtle balls and D dandelion balls. Drawing balls without replacement while they stay of the same color as the initial draw, else put back the last ball and repeat the process until all balls are drawn. The funny thing is that, unless M=0 or D=0, the probability to draw a myrtle ball at the end is always ½..! This can be easily checked by simulation (when M=2 and D=8)

r=function()sample(0:1,1,p=c(d,m))
for(t in 1:1e6){
  m=2;d=8
  i=r();m=m-!!i;d=d-!i
  while(!!m*d){
    j=r();i=ifelse(i==j,j,r())
    m=m-!!i;d=d-!i}
  F=F+(m>0)}
F/1e6

Now the proof that the probability is ½ is quite straightforward, for M=1 (or D=1). But I cannot find a quick fix for larger values. I thus reasoned by recursion, with the probability of emptying a given colour first is d!m!/(d+m)!, whatever the colour and whatever d>0,m>0. Hence half a chance to finish with myrtle. Any shorter sequence of a given colour reduces the value of either d or m, at which point we are using the recursion assumption that the probability is ½…