**Xi'an's Og » R**, and kindly contributed to R-bloggers)

**T**his 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.

**T**hen I started thinking rather than (merely and merrily) programming and realised that the terms of the sequence were all congruential to 8 modulo 17, hence that a power of 6, 7 or 8, also had to be congruencial to 8 if it was part of the sequence. Since

there is no solution in z if, for a given power *a*, all integers between 1 and 16 are not congruential to 8 modulo 17. Here is the check:

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

so this eliminates the power 8 (as well as 4 and 12), but not the power 7… Now, the check for the power 7 tells us the value of z is congruencial to 15 modulo 17, hence that the term of the sequence is equal to 15^{17}, 32^{17} or even more. This explains why I could not see any power equal to 7 in the first million trials. We are nonetheless lucky in that the first trial works:

> 15^7 [1] 170859375 > (170859375-8)/17 [1] 10050551 > 8+10050551*17-170859375 [1] 0 > ispower(8+10050551*17) $is [1] TRUE $pow [1] 15

For the power 6, the value of z is congruencial to 6 modulo 17, hence the term of the sequence is equal to 6^{6}, 23^{6} or even more. We again are lucky in that the first trial works:

> ispower(8+2744*17) $is [1] TRUE $pow [1] 6

Filed under: R, Statistics Tagged: congruence, prime number, R

**leave a comment**for the author, please follow the link and comment on his blog:

**Xi'an's Og » R**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...