Subscribe!

• Colins Blog

• Russian Peasant Multiplication
• FindingPerrinPseudoPrimes Part2
• FindingPerrinPseudoPrimes Part1
• The Unwise Update
• Miles Per Gallon
• Tracking An Item On Hacker News
• Hacker News User Ages
• Poking The Dusty Corners
• There Is No Time For This
• Publically Sharing Links
• Learning Times Tables
• Graceful Degradation
• Diagramming Maths Topics
• On The Rack
• Square Root By Long Division
• Beyond The Boundary
• Fill In The Gaps
• Software Checklist
• NASA Space Crews
• The Birthday Paradox
• The Trapezium Conundrum
• Revisiting The Ant
• The Ant And The Rubber Band
• Irrationals Exist
• Multiple Choice Probability Puzzle
• Random Eratosthenes
• Wrapping Up Square Dissection
• Dissecting A Square Part 2
• Dissecting A Circle
• Dissecting A Square
• An Oddity In Tennis
• Decision Tree For Tennis
• Decision Trees In Games
• A Matter Of Convention
• Do You Nourish Or Tarnish
• Binary Search Reconsidered
• Two Equals Four
• The Lost Property Office
• The Forgiving User Interface
• Setting Up RSS
• Withdrawing From Hacker News

• Random Writings.
• Colins Blog 2010
• Colins Blog 2009
• Colins Blog 2008
• Colins Blog 2007
• Colins Blog Before 2007

Fast Perrin Test - 2015/05/19

For context you probably want to read these first:

• FindingPerrinPseudoPrimes Part1
• FindingPerrinPseudoPrimes Part2
• Russian Peasant Multiplication

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 stdout.flush()

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.

Here:

$$ \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?

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?

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!

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$

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

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.