## shortened iterations [code golf]

Posted in Kids, pictures, Statistics, Travel with tags , , , , , , , , on October 29, 2019 by xi'an A codegolf lazy morning exercise towards finding the sequence of integers that starts with an arbitrary value n and gets updated by blocks of four as $a_{4k+1} = a_{4k} \cdot(4k+1)\\ a_{4k+2} = a_{4k+1} + (4k+2)\\ a_{4k+3} = a_{4k+2} - (4k+3)\\ a_{4k+4} = a_{4k+3} / (4k+4)$

until the last term is not an integer. While the update can be easily implemented with the appropriate stopping rule, a simple congruence analysis shows that, depending on n, the sequence is 4, 8 or 12 values long when $n\not\equiv 1(4)\\ n\equiv 1(4)\ \text{and}\ 3(n-1)+4\not\equiv 0(32)\\ 3(n-1)+4\equiv 0(32)$

respectively. But sadly the more interesting fixed length solution

~=rep #redefine function
b=(scan()-1)*c(32~4,8,40~4,1,9~3)/32+c(1,1,3,0~3,6,-c(8,1,9,-71,17)/8)
b[!b%%1] #keep integers only


ends up being longer than the more basic one:

a=scan()
while(!a[T]%%1)a=c(a,d<-a[T]*T,d+T+1,e<-d-1,e/((T<-T+4)-1))
a[-T]


where Robin’s suggestion of using T rather than length is very cool as T has double meaning, first TRUE (and 1) then the length of a…

## an attempt at code golf

Posted in Kids, R with tags , , , , , , , , , on May 15, 2019 by xi'an Having discovered codegolf on Stack Exchange a few weeks ago, I spotted a few interesting puzzles since then but only got the opportunity at a try over a quiet and rainy weekend (and Robin being on vacation)! The challenge was to write an R code for deciding whether or not a given integer n is congruent or not, when congruent means that it is the surface of a rectangle triangle with all three sides rational. The question included a pointer to the Birch and Swinnerton-Dyer conjecture as a mean to check congruence although the real solution was provided by Tunnell’s Theorem, which states that n is congruent if and only if the number of integer solutions to 2x²+y²+8z²=n is twice as much as the number of integer solutions to 2x²+y²+32z²=n if n is odd and  the number of integer solutions to 8x²+y²+16z²=n is twice as much as the number of integer solutions to 8x²+y²+64z²=n if n is even. Although this is only true for squared-free integers. (I actually spent more time on figuring out the exact wording of the theorem than on optimising the R code!)

My original solution

p=function(n){
for (i in(n:2)^2)if(n%%i<1)n=n/i
if(n%%2){d=8;f=2;g=16}else{d=2;f=1;g=8}
A=0;b=(-n:n)^2
for(x in d*b)for(y in x+f*b)for(z in g*b)
A=A+(y+z==n)-2*(y+4*z==n)
A==0}


was quite naïve, as shown by the subsequent improvements by senior players, like the final (?) version of Guiseppe:

function(n){b=(-n:n)^2
for(i in b[b>0])n=n/i^(!n%%i)
P=2^(n%%2)
o=outer
!sum(!o(y<-o(8/P*b,2*b,"+")/P-n,z<-16/P*b,"+"),-2*!o(y,4*z,"+"))}


exhibiting a load of code golf tricks, from using an anonymous function to renaming functions with a single letter, to switching from integers to booleans and back with the exclamation mark.

## Le Monde puzzle [#960]

Posted in Kids, R with tags , , , , , on April 28, 2016 by xi'an An arithmetic Le Monde mathematical puzzle:

Given an integer k>1, consider the sequence defined by F(1)=1+1 mod k, F²(1)=F(1)+2 mod k, F³(1)=F²(1)+3 mod k, &tc. [With this notation, F is not necessarily a function.] For which value of k is the sequence the entire {0,1,…,k-1} set?

This leads to an easy brute force resolution, for instance writing the R function

crkl<-function(k) return(unique(cumsum(1:(2*k))%%k))


where 2k is a sufficient substitute for ∞. Then the cases where the successive images of 1 visit the entire set {0,1,…,k-1} are given by

> for (i in 2:550) if (length(crkl(i))==i) print(i)
 2
 4
 8
 16
 32
 64
 128
 256
 512


which suspiciously looks like the result that only the powers of 2 k=2,2²,2³,… lead to a complete exploration of the set {0,1,…,k-1}. Checking a few series in the plane back from Warwick, I quickly found that when k is odd, (1) the sequence is of period k and (2) there is symmetry in the sequence, which means it only takes (k-1)/2 values. For k even, there is a more complicated symmetry, with the sequence being of period 2k, symmetric around its two middle values, and taking the values 1,2,..,1+k(2k+1)/4,..,1+k(k+1)/2. Those values cannot cover the set {0,1,…,k-1} if two are equal, which means an i(i+1)/2 congruent to zero modulo k, hence equal to k. This is clearly impossible when k is a power of 2 because i and i+1 cannot both divide a power of 2. I waited for the published solution as of yesterday’s and the complete argument was to show that when N=2p, the corresponding sequence [for N] is made (modulo p) of the sequence for p plus the same sequence translated by p. The one for N is complete only if the one for p is complete, which by recursion eliminates all cases but the powers of 2…

## Le Monde puzzle [#29]

Posted in R, Statistics with tags , , on July 29, 2011 by xi'an

This week, the puzzle from the weekend edition of Le Monde was easy to state: in the sequence (8+17n), is there a 6th power? a 7th? an 8th? If so, give the first occurrence. So I first wrote an R code for a function testing whether an integer is any power:

ispower=function(x){
ispo=FALSE
logx=log(x)
i=trunc(logx/log(2))
while((i>1)&&(!ispo)){
j=t=trunc(exp(logx/i))
while (t<x) t=j*t
ispo=(x==t)
if (!ispo){
j=t=j+1
while (t<x) t=j*t
ispo=(x==t)}
i=i-1}
list(is=ispo,pow=j)}


(The function returns the highest possible power.) Then I ran the thing over the first million of values of the sequence:

fib=8
for (j in 1:10^6){
fib=fib+17
tes=ispower(fib)
if (tes$is) print(c(fib,tes$pow,log(fib)/log(tes\$pow)))}


only to find that only the powers 2,3,6,10,11,19 were present among the first terms. Continue reading

## Parallel computation [permutations]

Posted in R, Statistics, University life with tags , , , , on February 20, 2011 by xi'an

François Perron is visiting me for two months from Montréal and, following a discussion about the parallel implementation of MCMC algorithms—to which he also contributed with Yves Atchadé in 2005—, he remarked that a deterministic choice of permutations with the maximal contrast should do better than random or even half-random permutations. Assuming p processors or threads, with p+1 a prime number, his solution is to take element (i,j) of the permutation table as (ij) mod (n+1): here are a few examples


> ((1:10)%*%t(1:10))%%11
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]    1    2    3    4    5    6    7    8    9    10
[2,]    2    4    6    8   10    1    3    5    7     9
[3,]    3    6    9    1    4    7   10    2    5     8
[4,]    4    8    1    5    9    2    6   10    3     7
[5,]    5   10    4    9    3    8    2    7    1     6
[6,]    6    1    7    2    8    3    9    4   10     5
[7,]    7    3   10    6    2    9    5    1    8     4
[8,]    8    5    2   10    7    4    1    9    6     3
[9,]    9    7    5    3    1   10    8    6    4     2
[10,]   10    9    8    7    6    5    4    3    2     1

> ((1:16)%*%t(1:16))%%17
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
[1,]    1    2    3    4    5    6    7    8    9    10    11
[2,]    2    4    6    8   10   12   14   16    1     3     5
[3,]    3    6    9   12   15    1    4    7   10    13    16
[4,]    4    8   12   16    3    7   11   15    2     6    10
[5,]    5   10   15    3    8   13    1    6   11    16     4
[6,]    6   12    1    7   13    2    8   14    3     9    15
[7,]    7   14    4   11    1    8   15    5   12     2     9
[8,]    8   16    7   15    6   14    5   13    4    12     3
[9,]    9    1   10    2   11    3   12    4   13     5    14
[10,]   10    3   13    6   16    9    2   12    5    15     8
[11,]   11    5   16   10    4   15    9    3   14     8     2
[12,]   12    7    2   14    9    4   16   11    6     1    13
[13,]   13    9    5    1   14   10    6    2   15    11     7
[14,]   14   11    8    5    2   16   13   10    7     4     1
[15,]   15   13   11    9    7    5    3    1   16    14    12
[16,]   16   15   14   13   12   11   10    9    8     7     6


which show that the scheme provides an interestingly diverse repartition of the indices. We certainly have to try this in the revision.