Finding Perrin Pseudo Primes_Part 1

   
Recent changes
Table of contents
Links to this page
FRONT PAGE / INDEX

Subscribe!
@ColinTheMathmo

My latest posts can be found here:
Previous blog posts:
Additionally, some earlier writings:

Finding Perrin Pseudo-primes:
Part 1

Some 20 years ago I was chatting with a friend of mine ( Andrew Lipson ) and he said something like the following:

This isn't an exact
quotation, but it's
close enough.

  • "A chap at work who does maths, but is not actually a mathematician, asked me the following.

    Consider the sequence $k(n)$ with
    • $k(1)=0,\;k(2)=2,\;k(3)=3,$ and
    • $k(n)=k(n-2)+k(n-3).$
    Why is it true that $n$ divides $k(n)$ if and only if $n$ is prime?"

My immediate response was "Well, it's not true." and Andrew's reply was swift and to the point:

OK then: Find a
counter-example.

So I started on that. I figured Andrew had already done the basic work of looking for a small counter-example, but I also figured that I would need to get warmed up and start to get a feel for the problem.

Before that, I need to put this into some historical context. This was when, if you wanted to do arithmetic on more than 32 bit numbers, you had to write your own arbitrary arithmetic package. It was before the WWW as we know it, although not pre-internet per se. So some of what follows will seem unnecessary. If you think that, I invite you to find what is common for the first 1000 counter-examples, except for a small number of them, and tell me which ones do not follow the pattern.

And so ...

n k(n) Prime Divides Confirmed
2
2
True
True
Yes
3
3
True
True
Yes
4
2
False
False
Yes
5
5
True
True
Yes
6
5
False
False
Yes
7
7
True
True
Yes
8
10
False
False
Yes
9
12
False
False
Yes
10
17
False
False
Yes
11
22
True
True
Yes
12
29
False
False
Yes
13
39
True
True
Yes
14
51
False
False
Yes
15
68
False
False
Yes
16
90
False
False
Yes
17
119
True
True
Yes
18
158
False
False
Yes
19
209
True
True
Yes
20
277
False
False
Yes
21
367
False
False
Yes
22
486
False
False
Yes
23
644
True
True
Yes
24
853
False
False
Yes
25
1130
False
False
Yes

and so on ...
OK, so working it by hand we can see that it's true up to 25, but we can also see that k(n) is starting to grow quite rapidly. Some messing about by hand shows that it's going to be exponential, just as the Fibonacci sequence is exponential. Not a surprise, really, since this is very Fibonacci-like, with each term defined by adding two previous terms.

We can work out the ratio between successive terms and see that it starts to converge quite rapidly to about 1.324736... In fact we can show that:

$$k(n)=\alpha^n+\beta^n+\gamma^n$$

where $\alpha,$ $\beta,$ and $\gamma$ are solutions to the equation:

$$x^3=x+1$$

(Usually we need constants in front of the powers of the roots, but in this case the constants all turn out to be 1).

Skip this
if you
want to.
That's not magic, that's called the "characteristic equation" and you can see that the coefficients are related to the definition.

No? Can't see that? Well the definition is:

$$k(n)=a.k(n-1)+b.k(n-2)+c.k(n-3)$$

where a=0, b=1, and c=1. Then the characteristic equation is $x^3=a.x^2+b.x+c.$

Anyway ...

The roots to that cubic are:

  • $\alpha\approx{}1.324736\ldots$
  • $\beta\approx{}-0.66236-0.56228{}i$
  • $\gamma\approx{}-0.66236+0.56228{}i$

These two last have absolute value less than 1, so as n gets bigger, $\beta^n$ and $\gamma^n$ approach 0. Since $k(n)=\alpha^n+\beta^n+\gamma^n,$ with $\beta^n$ and $\gamma^n$ disappearing, so $k(n)\rightarrow\alpha^n$ as n gets larger.

Now the number of digits in a number m (base 10) is given by log(m) (base 10), so that means the number of digits in k(n) is log(k(n)) which is roughly $\log(\alpha^n),$ which in turn is $n.\log(\alpha).$

Now $\log(\alpha)$ is about 0.122..., so by the time n is 100, k(n) will have 12 digits (actually 13, because we should round up), and it will only get worse.

So we can't go far with this by hand. But that's OK, because what we want to know is this:

  • When n is prime,
    • does n divide k(n) ?

The sequence is called
the Perrin Sequence,
so I call this the
"Perrin Test."

 
#!/usr/bin/python

from time import clock

for limit_power in [0, 1, 2, 3, 4]:

    n = 2
    limit = 4**limit_power
    t_start = clock()

    while clock()-t_start < limit:

        a,b,c = 3,0,2
        for i in xrange(n):
            a,b,c = b,c,(a+b) % n

        n += 1

    print limit, n, float(n)/2**limit_power
If we compute k(n) modulo n then it will divide if the answer is 0, and because we're working modulo n the numbers in the intermediate calculation don't grow too quickly. The only problem then is that after each n, we have to start all over again for the next one. That means that the time taken grows quadratically, and that will kill us fairly quickly - we won't get that far.

So here we have some Python code to compute k(n) mod n for as long as allowed by some time limit, and then we print how far we got. We're not actually doing any of the searching for a counter-example in this code, since we're not actually seeing whether the numbers are predicted to be prime, nor are we seeing if they actually are prime. This is just showing that it's going to be a slow process if we use this method.

Secs
n
n/sqrt(t)
1
3133
3133.0
4
6229
3114.5
16
12562
3140.5
64
25179
3147.4
256
50642
3165.1

The code multiplies the time limit by 4 each time, and we can see that we get twice as far. That means that to get twice as far takes four times as long, so, as expected, the time taken grows quadratically. The timings on my machine are shown at left.

 
#!/usr/bin/python

from time import clock
from sys import argv

def is_prime_q(n):

    if n in [2, 3, 5, 7]: return True
    if n < 11           : return False

    if n % 2 == 0       : return False
    if pow(2,n-1,n) != 1: return False
    if n % 3 == 0       : return False
    if pow(3,n-1,n) != 1: return False
    if n % 5 == 0       : return False
    if pow(5,n-1,n) != 1: return False

    d = 7
    while d*d <= n:
        if (n % d) == 0: return False
        d += 2

    return True

def is_perrin_q(n):

    if n in [2, 3]: return True

    a0, a1, a2 = 3, 0, 2

    i = n
    while i > 0:

        a0, a1, a2 = a1, a2, (a0+a1) % n
        i -= 1

    return a0 == 0

# ======================================

n = 2
t = 1

while clock() < int(argv[1]):

    is_prime = is_prime_q(n)
    is_perrin = is_perrin_q(n)
    if is_prime != is_perrin:
        print n, is_prime, is_perrin

    if t < clock():
        print 'Tested to %d in %d secs' % (n,t)
        t += 1

    n += 1
As you can see, it takes 4 minutes to get to 50 thousand, so that's not getting much done. We can derive a formula for how far we get in a given time:

$$n\;\approx\;3150\sqrt(t)$$

So here we have the code to start the testing process. We have a simple test for primality, first testing for being a small prime, then using Fermat's Little Theorem, and finally performing trial division up to the square root of the number being tested.

Then we have the routine to test the Perrin condition. We start with the first few terms in the Perrin sequence, then advance up the sequence, working modulo n to avoid having the numbers get stupidly large. In the end, we return True if the relevant term in the sequence is zero, indicating that n divides k(n).

The main loop starts at n=2 and then tests to see if the number is prime, and whether it passes the Perrin test. If these answers don't match, then the number is printed, and is a counter-example.

Specific measurements give:

  • Roughly 104 per second when n=100k
  • Roughly 42 per second when n=250k
  • Roughly 10 per second when n=1000k

We also print every second how far we've got so far, so we can measure how fast the code runs, and then see how we can improve it. Running it for 1000 seconds and fitting an equation, we get an approximate formula for its speed:

$$n\;\approx\;4575\sqrt{t}$$

This is all very rough and ready, and will fall apart as we go significantly further, but it's enough to get an idea of what's going on. We can see that getting to a million will take about 13 hours, so we really, really want to speed this up.

So what do we do? How can we get further, faster?

In Part 2 we'll look at algorithmic improvements.


<<<< Prev <<<<
The Unwise Update
:
>>>> Next >>>>
FindingPerrinPseudoPrimes Part2 ...


You should follow me on twitter @ColinTheMathmo

Comments

I've decided no longer to include comments directly via the Disqus (or any other) system. Instead, I'd be more than delighted to get emails from people who wish to make comments or engage in discussion. Comments will then be integrated into the page as and when they are appropriate.

If the number of emails/comments gets too large to handle then I might return to a semi-automated system. That's looking increasingly unlikely.


Contents

 

Links on this page

 
Site hosted by Colin and Rachel Wright:
  • Maths, Design, Juggling, Computing,
  • Embroidery, Proof-reading,
  • and other clever stuff.

Suggest a change ( <-- What does this mean?) / Send me email
Front Page / All pages by date / Site overview / Top of page

Universally Browser Friendly     Quotation from
Tim Berners-Lee
    Valid HTML 3.2!