A week ago I wrote about Perrin numbers, numbers *P*_{n} defined by a recurrence relation similar to Fibonacci numbers. If *n* is prime, *P*_{n} mod *n* = 0, and the converse is nearly always true. That is, if *P*_{n} 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 *P*_{n} and then see whether it has remainder 0 when you divide by *n*. The problem is that *P*_{n} grows exponentially with n. In fact, *P*_{n} is approximately ρ^{n} where ρ = 1.3247… is the plastic number. This means that *P*_{n} has about *n* log_{10} ρ 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 *P*_{n} mod *n* rather than *P*_{n} *per 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 *P*_{n} mod *n* for each *n* up to *N* = 10^{9}. This only requires ordinary arithmetic.

However, this approach takes O(*N*^{2}) time unless we’re clever. We have to compute *P*_{n} 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 *P*_{n} mod *n*. This requires calculating the *n*th 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(*N*^{2}) time. This is enough to make the search for pseudoprimes less than a billion practical.