Fast Perrin Test

Recent changes
Table of contents
Links to this page


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

Fast Perrin Test - 2015/05/19

For context you probably want to read these first:

So we've got scaffolding to look for Perrin Pseudo-Primes (PPPs), assuming any exist (which they do) but as we run the existing code we find that it's spending pretty much all its time in the test as to whether n divides k(n).

n = 2
for n in xrange(2,1+10*5):

    is_prime = is_prime_q(n)
    is_perrin = is_perrin_q(n)

    if is_prime != is_perrin:
        print n
At right is the main code loop for finding a PPP. We start with n=2. We ask if it's prime, we ask if it passes the Perrin test, and if those two things aren't the same, we print n as a PPP. Using the obvious Perrin test this takes about 60 seconds to get up to $n=10^5.$ Worse, though, pretty much all of the time, over 99%, is spent in the Perrin test routine.

We have to make that faster.

$K_0 = \begin{bmatrix} k_0 \\ k_1 \\ k_2 \\ \end{bmatrix}, K_1 = \begin{bmatrix} k_1 \\ k_2 \\ k_3 \\ \end{bmatrix}, K_n = \begin{bmatrix} k_n \\ k_{n+1} \\ k_{n+2} \\ \end{bmatrix} $
So here's a really clever idea. We start, not necessarily obviously, by writing three consecutive terms from the Perrin sequence as a vector. Here at left are some examples. The question is, given one of these vectors, can we somehow transform it into another?

As soon as we ask that question we think of matrices, because matrices represent transformations from vectors to vectors, so we want a matrix $T$ such that $K_1=T.K_0.$ Such a matrix is easy to construct.


$$ \begin{bmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 1 & 0 \\ \end{bmatrix} \begin{bmatrix} k_0 \\ k_1 \\ k_2 \\ \end{bmatrix} = \begin{bmatrix} k_1 \\ k_2 \\ k_0+k_1 \\ \end{bmatrix} = \begin{bmatrix} k_1 \\ k_2 \\ k_3 \\ \end{bmatrix} $$

Applying $T$ lots of times steps us through the vectors, so we get:

$$ K_n = T^nK_0 $$

def multiply_matrices(A,B,n):
    A0,A1,A2,A3,A4,A5,A6,A7,A8 = A
    B0,B1,B2,B3,B4,B5,B6,B7,B8 = B
    return [
        ( A0*B0 + A1*B3 + A2*B6 ) % n,
        ( A0*B1 + A1*B4 + A2*B7 ) % n,
        ( A0*B2 + A1*B5 + A2*B8 ) % n,
        ( A3*B0 + A4*B3 + A5*B6 ) % n,
        ( A3*B1 + A4*B4 + A5*B7 ) % n,
        ( A3*B2 + A4*B5 + A5*B8 ) % n,
        ( A6*B0 + A7*B3 + A8*B6 ) % n,
        ( A6*B1 + A7*B4 + A8*B7 ) % n,
        ( A6*B2 + A7*B5 + A8*B8 ) % n,

def compute_perrin_number(n):

    power = n
    C = [ 1, 0, 0, 0, 1, 0, 0, 0, 1 ]
    A = [ 0, 1, 0, 0, 0, 1, 1, 1, 0 ]

    while power>0:
        if (power & 1)==1:
            C = multiply_matrices(C,A,n)
        A = multiply_matrices(A,A,n)
        power = power/2

    return ( C[0]*3 + C[2]*2 ) % n

def is_perrin_q(n):

# Pre-filtering code omitted ...

    return compute_perrin_number(n) == 0
So if we want to compute the $n^{th}$ Perrin term we can simply compute $K_n=T^nK_0$ and the first element of $K_n$ is the term we want.

Then we can use Russian Peasant Multiplication to compute that quickly and efficiently. The code shown here does that for us. There are as yet a few micro-optimisations we can make, but we'll leave that until after we've done the main algorithmic improvements, because time taken now could be wasted as the code changes. For similar reasons I'm not using any libraries such as sage, numpy, or scipy to do the sums. At this stage it's worth keeping it all "on the surface" and visible.

So what about running this?

Remember, it was taking about 60 seconds to get to $10^5$. How long does this version take?

1.19 seconds

Well that's not bad. (!)

We have been running the code for 100 seconds to see how far it gets. Previously it got to about n=109159. What about now?

n = 6594026

Yup, we get about 60 times as far.

But what's more, we've got our first Perrin Pseudo-Prime. And in fact we've got two!

271441, 904631

Interestingly, in both cases they are predicted to be prime by the Perrin test whereas they are in fact composite:

  • $271441 = 521^2$
  • $904631 = 7\times13\times9941$

If it's always true that a PPP is a composite that's incorrectly predicted to be prime (and it is true - the clue is in the name) then there are more optimisations we can make. But for now we can look at the timing and see that in the 100 seconds we got as far as $n=6594026,$ with 22.61 seconds in the prime testing routine, and 66.93 seconds in the Perrin testing routine.

Next, we see if we can take advantage of the mathematical observation about the failures, and whether we can further reuse calculations we've already done.

<<<< Prev <<<<
Russian Peasant Multiplication
>>>> Next >>>>
Sieve Of Eratosthenes In Python ...

You should follow me on twitter @ColinTheMathmo


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.



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!