Perrin Pseudoprimes
May 29, 2015
The Perrin Numbers Pn are defined as P0 = 3, P1 = 0, P2 = 2, and Pn+3 = Pn+1 + Pn for n > 2. If Pn (mod n) ≡ 0, then n is most likely prime. Perrin’s formula always reports that a prime number is prime, but sometimes reports a composite number is prime, though seldom: there are only two pseudoprimes, 271441 and 904631, less than a million.
The Perrin sequence A001608 is computed by sequential addition. An individual member of the Perrin sequence can be computed by matrix exponentiation followed by matrix multiplication:
![P_n = \left[ \begin{array}{ccc} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 1 & 0 \end{array} \right] ^n \left[ \begin{array}{c} 3 \\ 0 \\ 2 \end{array} \right]](https://s0.wp.com/latex.php?latex=P_n+%3D+%5Cleft%5B+%5Cbegin%7Barray%7D%7Bccc%7D+0+%26+1+%26+0+%5C%5C+0+%26+0+%26+1+%5C%5C+1+%26+1+%26+0+%5Cend%7Barray%7D+%5Cright%5D+%5En+%5Cleft%5B+%5Cbegin%7Barray%7D%7Bc%7D+3+%5C%5C+0+%5C%5C+2+%5Cend%7Barray%7D+%5Cright%5D+&bg=ffffff&fg=555555&s=0&c=20201002)
The terms of the Perrin sequence grow exponentially at a rate of 1.32471795, which is known at the plastic number.
The Perrin pseudoprimality test can be implemented using matrix exponentiation followed by matrix multiplication with all operations performed modulo n. Sloane gives a list of Perrin pseudoprimes at A013998.
Your task is to write a function to determine if a number is a Perrin pseudoprime and to find all Perrin pseudoprimes less than a million. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.
Programming Praxis 
In Python2
# Cached version def cached_perrin(up_until): cache = {0:3, 1:0, 2:2} max_cache_size = 100000 #approx 1 GB memory, perhaps better to just use itertools lru_function def p(n): if n < max_cache_size and n in cache: return cache[n] else: v = p(n-2) + p(n-3) if n < max_cache_size: cache[n] = v return v count = 0 for n in range(1, up_until): if p(n) % n == 0: count +=1 return count # for perrin(n) we only need to keep track of the last 3 values def iterative_perrin(up_until): last_3 = [3,0,2] count = 2 # 2 primes in the first 3 numbers for n in range(3, up_until): v = last_3[(n-3)%3] + last_3[(n-2)%3] last_3[n%3] = v if v % n == 0: count +=1 return count up_until = 1000000 print "Total =", cached_perrin(up_until) print "Total =", iterative_perrin(up_until) Second method is much faster, both print Total = 78501See also Colin Wright’s blog: http://www.solipsys.co.uk/new/FastPerrinTest.html?ColinsBlog
In Perl, but cheating since the work is done in C:
use ntheory ":all"; forcomposites { say if is_perrin_pseudoprime($_) } 1e6;Runs in 0.7s on my Macbook with the latest code.
It does modular exponentiation of the generating matrix. If the input is less than half a half-native-word, this can be done easily, otherwise care must be taken with proper mulmod and addmod code. The trace is compared to zero.
Based on some of Colin’s notes on his blog, I recently added a pre-test that, for numbers divisible by small primes, checks the divisibility of precalculated small Perrin sequences. This is only about 60 bytes of data to check divisibility for 2,3,5,7,11,13. I’m hoping his future entries give more ideas for speeding this up. As we’d expect, the BPSW test is faster and is known to have *no* counterexamples for any 64-bit inputs.
@Praxis, In the first sentence “The Perrin Numbers Pn are defined as P0 = 3, P1 = 2, P2 = 2”, I think it should be “P1 = 0”.
@Mike: You are correct. Fixed. Thank you.
@Rutger, cached_perrin() indicates 1 is prime. iterative_perrin() is double counting 3. These are both easily fixed. Since the first Perrin pseudoprime is 271441, we should expect a count to 10,000 to yield 1229 and a count to 100,000 to yield 9592. The count to 1M should get 78500 (because of the two pseudoprimes).