## 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.

In Python2

See 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:

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).