Searching for Perrin pseudoprimes

A week ago I wrote about Perrin numbers, numbers Pn defined by a recurrence relation similar to Fibonacci numbers. If n is prime, Pn mod n = 0, and the converse is nearly always true. That is, if Pn mod n = 0, n is usually prime. The exceptions are called Perrin pseudoprimes.

Matt McIrvin wrote an excellent post explaining how to compute Perrin pseudoprimes. Here I’m just going to elaborate on a couple points in his post.

Matt’s first point is that if you want to search for Perrin pseudoprimes, the most direct approach won’t get you very far. The obvious thing to do is compute Pn and then see whether it has remainder 0 when you divide by n. The problem is that Pn grows exponentially with n. In fact, Pn is approximately ρn where ρ = 1.3247… is the plastic number. This means that Pn has about n log10 ρ digits. So searching for pseudoprimes less than one billion would require working with numbers with over 100,000,000 digits. This can be done, but it’s slow and unnecessary.

Since the goal is to compute Pn mod n rather than Pnper se, we can carry out all calculations mod n and avoid extended precision arithmetic as long as n itself can fit in an ordinary precision integer. If we want to find pseudoprimes less than one billion, we calculate Pn mod n for each n up to N = 109. This only requires ordinary arithmetic.

However, this approach takes O(N2) time unless we’re clever. We have to compute Pn mod n separately for each n, and the most direct approach takes n steps. This leads to Matt’s second point: use matrix multiplication (mod n) to calculate Pn mod n. This requires calculating the nth power of a 3×3 matrix, which can be done in O(log n) time using fast exponentiation. This makes the search for pseudoprimes less than N require O(N log N) rather than O(N2) time. This is enough to make the search for pseudoprimes less than a billion practical.