[seqfan] Re: Proving a number to be prime using the Tower of Hanoi.

[Rick Shepherd]

Hi Daniel, Robert, and SeqFans,

Interesting observations, Daniel.  Here are some of mine from a few days
ago.
Robert, you posted about 341, etc., before I had time to.

I haven't worked with continued fractions much, so hadn't thought of them as
necessarily being useful for probable-prime tests.  In this case (at
least), perhaps
they're even better as a sort of "probable-squarefree" confirmation test,
albeit possibly "expensive."
(Of course there are squarefree numbers that do not produce continued
fractions of this form.)

The continued fraction for A000295(n+1) / A000295(n) has the form [2; (2^n
- n - 2) / n, n]
if and only if n > 2 and n is a term of A015919.

(I'll be correcting the "formula" in A015919 to indicate that A006935 is
also a subset.)

In particular, for n > 1, A000295(prime(n) + 1) / A000295(prime(n)) = [2;
A164740(n), prime(n)],
where prime(n) is the n-th prime.

Examining 3 <= n <= 500000, the continued fractions having this form occur
exactly
for all 41537 primes and all 172 (squarefree) Fermat pseudoprimes to base 2
(elements of
A001567) and both even pseudoprimes (161038 and 215326 in A006935) in this
range.

At least for "smaller numbers," n is highly likely to be prime
(respectively, odd, or squarefree) if its
corresponding continued fraction is of this form.  For 3 <= n <= 500000,
99.58% of these cases are prime;
99.995% are odd, and 100% are squarefree.

See also A158358, "Pseudoprimes to base 2 that are *not* squarefree," the
first term of which
is 1194649 (actually a square, 1093^2, where 1093 is one of only two known
Wieferich primes (A001220)).
These pseudoprimes are relatively rare.

Rick

On Sat, Oct 18, 2014 at 2:32 AM, Robert Munafo <mrob27 at gmail.com> wrote:

> On Tue, 14 Oct 2014 11:17:23, Daniel Joyce <hlauk.h.bogart at gmail.com>
> wrote:
> > The tower of Hanoi and the sum of its pieces in all the moves it takes to
> move to a different peg is related to the primes and a sequence in OEIS.
> > ...
>
> Did you try 341? :-)
>
> Looking at A000295, note the direct formula a(n)=2^n - n - 1, and the
> recurrence formula a(1)=0, a(n)=2*a(n-1)+n-1, n>1.
>
> Now consider a given "possible prime" p. Daniel has suggested looking at
> the continued fraction (c.f.) expansion of the fraction a(p+1)/a(p). For
> example, is p is 5, we are considering the c.f. expansion of a(6)/a(5),
> which is 57/26.
>
> For any given p, a(p) will be 2^p-p-1 (example: 32-5-1 = 26), and the next
> term a(p+1) is 2*a(p) + (p+1)-1, or more simply 2*a(p)+p. (example:
> 2*26+6-1 = 52+5 = 57)
>
> In general, we want to find the terms of the c.f. expansion of
> [2*(2^p-p-1)+p]/(2^p-p-1). You can see right away that this is 2 +
> p/(2^p-p-1). So the first term of the c.f. is 2. Continuing, we have to
> turn it into
>
>   2 + 1/[(2^p-p-1)/p]
>
> and consider what happens with the next term.
>
> If p is, in fact, prime, then by Fermat's Little Theorem, 2^p-2 is a
> multiple of p. If we had "(2^p - p - 2)/p" it would be an integer and our
> c.f. would be done. Instead we are too big by 1, so there will be 1/p left
> over. In other words, (2^p - p - 1)/p = k + 1/p for some integer k. That
> 1/p results in one more c.f. term and then the c.f. terminates: 2 + 1/(k +
> 1/(p + 0)).
>
> If p is composite, then *usually* 2^p-2 is not a multiple of p, and the
> c.f. ends up having more terms. But not always... there are composites for
> which 2^p-2 *is* a multiple of p. These are called "*Fermat pseudoprimes*"
> to base 2.
>
> I'll use Maxima, because it has a nice c.f. function and you don't have to
> set a precision limit. Here is what happens when p=5 (n is numerator, d is
> denominator)
>
> (%i8) p:5;
> (%o8)                                  5
> (%i9) d:2^p-p-1;
> (%o9)                                 26
> (%i10) n:2*d+p;
> (%o10)                                57
> (%i11) cf(n/d);
> (%o11)                             [2, 5, 5]
>
>
> So far, so good. Now we'll let p be the composite number 9, and then the
> known prime 37:
>
> (%i1) p:9;
> (%o1)                                  9
> (%i2) d:2^p-p-1;
> (%o2)                                 502
> (%i3) n:2*d+p;
> (%o3)                                1013
> (%i4) cf(n/d);
> (%o4)                          [2, 55, 1, 3, 2]
>
> (%i12) p:37;
> (%o12)                                37
> (%i13) d:2^p-p-1;
> (%o13)                           137438953434
> (%i14) n:2*d+p;
> (%o14)                           274877906905
> (%i15) cf(n/d);
> (%o15)                        [2, 3714566309, 37]
>
>
> Everything seems fine so far. But here's what happens with 341:
>
> (%i16) p:341;
> (%o16)                                341
> (%i17) d:2^p-p-1;
> (%o17)
> 44794894843556084211148845611368885562432909944692990697999782019275837\
> 42360321890761754986543214231210
> (%i18) n:2*d+p;
> (%o18)
> 89589789687112168422297691222737771124865819889385981395999564038551674\
> 84720643781523509973086428462761
> (%i19) cf(n/d);
> (%o19) [2,
> 1313633279869679888889995472474160866933516420665483598181811789421\
> 5788100763407304286671514789484549, 341]
> (%i20) factor(341);
> (%o20)                               11 31
>
>
> In line %o19 we see that the c.f. expansion terminates after just 3 terms
> (the middle one is big), but in line %o20 we see that 341 = 11 * 31.
>
> On the plus side, if you actually manage to complete a Tower of Hanoi with
> 341 disks, you will have out-lived the Universe, which is quite an
> accomplishment, primes or no primes.
>
> --
>   Robert Munafo  --  mrob.com
>   Follow me at: gplus.to/mrob - fb.com/mrob27 - twitter.com/mrob_27 -
> mrob27.wordpress.com - youtube.com/user/mrob143 - rilybot.blogspot.com
>
> _______________________________________________
>
> Seqfan Mailing list - http://list.seqfan.eu/
>