Archive for pracma

Le Monde puzzle [#840]

Posted in Books, Kids, R with tags , , , , on November 23, 2013 by xi'an

Another number theory Le Monde mathematical puzzles:

Find 2≤n≤50 such that the sequence {1,…,n} can be permuted into a sequence such that the sum of two consecutive terms is a prime number. 

Now this is a problem with an R code solution:

library(pracma)
foundsol=TRUE
N=2
while (foundsol){

  N=N+1
  noseq=TRUE
  uplim=10^6
  t=0
  while ((t<uplim)&&(noseq)){

    randseq=sample(1:N)
    sumseq=randseq[-1]+randseq[-N]
    noseq=min(isprime(sumseq))==0
    t=t+1
    }

  foundsol=!noseq
  if (!noseq){
   lastsol=randseq}else{ N=N-1}
  }

which returns the solution as

> N
[1] 12
> lastsol
 [1]  6  7 12 11  8  5  2  1  4  3 10  9

and so it seems there is no solution beyond N=12…

However, reading the solution in the next edition of Le Monde, the authors claim there are solutions up to 50. I wonder why the crude search above fails so suddenly, between 12 and 13! So instead I tried a recursive program that exploits the fact that subchains are also verifying  the same property:

findord=function(ens){

  if (length(ens)==2){
    sol=ens
    foundsol=isprime(sum(ens))}
  else{
    but=sample(ens,1)
    nut=findord(ens[ens!=but])
    foundsol=FALSE
    sol=ens
    if (nut$find){
      tut=nut$ord
      foundsol=max(isprime(but+tut[1]),
         isprime(but+tut[length(tut)]))
      sol=c(tut,but)
      if (isprime(but+tut[1]))
         sol=c(but,tut)
      }
  }
  list(find=foundsol,ord=sol)
}

And I ran the R code for N=13,14,…

> stop=TRUE
> while (stop){
+   a=findord(1:N)
+   stop=!(a$find)}

until I reached N=20 for which the R code would not return a solution. Maybe the next step would be to store solutions in N before moving to N+1. This is just getting  me too far from a mere Saturday afternoon break.

Le Monde puzzle [#825]

Posted in Books, Kids, pictures, R with tags , , , , , on June 19, 2013 by xi'an

Yet another puzzle which first part does not require R programming, even though it is a programming question in essence:

Given five real numbers x1,…,x5, what is the minimal number of pairwise comparisons needed to rank them? Given 33 real numbers, what is the minimal number of pairwise comparisons required to find the three largest ones?

I do not see a way out of considering the first question as the quickest possible sorting of a real sample. Using either quicksort or heapsort, I achieve sorting the 5 numbers in exactly 6 comparisons for any order of the initial sample. (Now, there may be an even faster way based on comparing partial sums first… I just do not see how!) Update: Oops! I realised I made my reasoning based on a reasonable case, the correct answer is indeed 7!!!

For the second part, let us start from the remark that 32 comparisons are needed to find the largest number, then at most 31 for the second largest, and at most 30 for the third largest (since we can take advantage of the partial ordering resulting from the determination of the largest number). This is poor. If I instead use a heap algorithm, I need O(n log{n}) comparisons to build this binary tree whose parents are always larger than their siblings, as in the above example. (I can produce a sort of heap structure, although non-binary, in an average 16×2.5=40 steps. And a maximum 16×3=48 steps.) The resulting tree provides the largest number (100 in the above example) and at least the second largest number (36 in the above). To get the third largest number, I first need a comparison between the one-before-last terms of the heap (19 vs. 36 in the above), and one or two extra comparisons (25 vs. 19 and maybe 25 vs. 1 in the above). (This would induce an average 1.5 extra comparison and a maximum 2 extra comparisons, resulting in a total of 41.5 average and 49.5 maximum comparisons with my sub-optimal heap construction.)  Once again, using comparisons of sums may help in speeding up the process, for instance comparing numbers by groups of 3, but I did not pursue this solution…

quick3If instead I try to adapt quicksort to this problem, I can have a dynamic pivot that always keep at most two terms above it, providing the three numbers as a finale result. Here is an R code to check its performances:

quick3=function(x){

comp=0
i=1
lower=upper=NULL
pivot=x[1]

for (i in 2:length(x)){

if (x[i]<pivot){ lower=c(lower,x[i])
}else{

upper=c(upper,x[i])
if (length(upper)>1) comp=comp+1}

comp=comp+1

if (length(upper)==3){

pivot=min(upper)
upper=sort(upper)[-1]
}}

if (length(upper)<3) upper=c(pivot,upper)

list(top=upper,cost=comp)
}

When running this R code on 10⁴ random sequences of 33 terms, I obtained the following statistics, I obtained the following statistics on the computing costs

</p>
<p style="text-align: justify;">> summary(costs)
Min.   1st Qu.  Median  Mean    3rd Qu.  Max.
32.00   36.00   38.00   37.89   40.00   49.00

and the associated histogram represented above. Interestingly, the minimum is the number of comparisons needed to produce the maximum!

Reading the solution in Le Monde in the train to London and Bayes 250, I discovered that the author suggests a 7 comparison solution in the first case [compare A and B, C and D = 2 comparisons; if A>B and C>D, compare A and C and get, say A>C>D = 1 comparison; insert E in this ordering by comparing with C and then A or D = 2 comparisons, obtaining e.g. A>E>C>D; conclude by inserting B by first comparing with C then with D or E = 2 comparisons] and a 41 comparison solution in the second case. He is (or I am!) I was clearly mistaken in the first case while the quick3 algorithm does 41 or less most of the time (90%)  but not always.

Le Monde puzzle [#820]

Posted in Books, Kids, R with tags , , , , , on May 15, 2013 by xi'an

The current puzzle is… puzzling:

Given the set {1,…,N} with N<61, one iterates the following procedure: take (x,y) within the set and replace the pair with the smallest divider of x+y (bar 1). What are the values of N such that the final value in the set is 61?

I find it puzzling because the way the pairs are selected impacts the final value. Or not, depending upon N. Using the following code (with factors() from the pracma package):

library(pracma)
endof=function(N){
  coll=1:N
  for (t in 1:(N-1)){

    pair=sample(1:length(coll),2)
    dive=min(factors(sum(coll[pair])))
    coll=coll[-pair]
    coll=c(coll,dive)
    }
  print(dive)
  }

I got:

> for (t in 1:10) endof(10)
[1] 5
[1] 3
[1] 3
[1] 5
[1] 7
[1] 5
[1] 5
[1] 7
[1] 3
[1] 3> for (t in 1:10) endof(16)
[1] 2
[1] 2
[1] 2
[1] 2
[1] 2
[1] 2
[1] 2
[1] 2
[1] 2
[1] 2

For N of the form 4k or 4k-1, the final number is always 2 while for N‘s of the form 4k-2 and 4k-3, the final number varies, sometimes producing 61’s. Although I could not find solutions for N less than 17…  Looking more closely into the sequence leading to 61, I could not see a pattern, apart from producing prime numbers as, in, e.g.

61 = 2 + [12 +  (4 + {14 + [13 + 16]})]

for N=17.  (Another puzzle is that 61 plays no particular role: a long run of random calls to endof() return all prime numbers up to 79…)

Udate: Looking at the solution in today’s edition, there exist a solution for N=13 and a solution for N=14. Even though my R code fails to spot it. Of course, an exhaustive search would be feasible in these two cases.  (I had also eliminated values below as not summing up to 61.) The argument for eliminating 4k and 4k-1 is that there must be an odd number of odd numbers in the collection, otherwise, the final number is always 2.