Perrin Pseudoprimes
May 29, 2015
We begin by writing a function that computes the Perrin numbers sequentially and tests each using a Miller-Rabin pseudoprime test:
(define (perrin-pseudo-primes n)
(let loop ((p-2 3) (p-1 0) (p 2) (i 3))
(when (< i n)
(let ((p+1 (+ p-2 p-1)))
(when (and (zero? (modulo p+1 i))
(not (prime? i)))
(display i) (newline))
(loop p-1 p p+1 (+ i 1))))))
> (time (perrin-pseudo-primes 1000000))
271441
904631
(time (perrin-pseudo-primes 1000000))
3077 collections
108124 ms elapsed cpu time, including 358 ms collecting
108680 ms elapsed real time, including 346 ms collecting
26597533392 bytes allocated, including 26600269744 bytes reclaimed
> (time (perrin-pseudo-primes 2000000))
271441
904631
(time (perrin-pseudo-primes 2000000))
11675 collections
415446 ms elapsed cpu time, including 1683 ms collecting
424037 ms elapsed real time, including 1609 ms collecting
103799938080 bytes allocated, including 103791257680 bytes reclaimed
This appears to take quadratic time, as doubling the input caused the elapsed time to quadruple. That makes sense: even though the underlying algorithm is linear (it performs a linear number of additions), the numbers of the Perrin sequence grow exponentially and the bit-cost of performing the various arithmetic operations begins to be a factor.
Next we look at the matrix methods for computing Perrin numbers, both the standard method that computes the entire Perrin number and the modular method that computes the n th Perrin number mod n:
(define (perrin n)
(if (= n 0) 3 (if (= n 1) 0 (if (= n 2) 2
(vector-ref (vector-ref
(matrix-multiply
(matrix-power perrin-matrix (- n 2))
perrin-init)
2) 0)))))
> (perrin 200)
2658793989922287946990250
(define (perrin-mod n)
(if (= n 0) 3 (if (= n 1) 0 (if (= n 2) 2
(vector-ref (vector-ref
(matrix-multiply-mod
(matrix-power-mod perrin-matrix (- n 2) n)
perrin-init n)
2) 0)))))
> (perrin-mod 200)
50
> (modulo (perrin 200) 200)
50
The two methods are identical except for the use of modular arithmetic. Then we write a function that identifies Perrin primes and use it to compute the Perrin pseudoprimes:
(define (perrin-prime? n) (zero? (perrin-mod n)))
> (perrin-prime? 271441)
#t
> (prime? 271441)
#f
(define (perrin-pseudo-primes n)
(do ((i 2 (+ i 1))) ((< n i))
(when (and (perrin-prime? i)
(not (prime? i)))
(display i) (newline))))
> (time (perrin-pseudo-primes 1000000))
271441
904631
(time (perrin-pseudo-primes 1000000))
28595 collections
618965 ms elapsed cpu time, including 1912 ms collecting
631582 ms elapsed real time, including 2122 ms collecting
240949881088 bytes allocated, including 240947324672 bytes reclaimed
That’s much longer than the additive method, but the order of growth is linear rather than quadratic, so it will soon catch up. You can run the program at http://ideone.com/SOZgf8, where you will also see the various matrix functions and the Miller-Rabin primality tester.
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).