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

Pages: 1 2

### 5 Responses to “Perrin Pseudoprimes”

1. Rutger said

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 = 78501
```
2. danaj said

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:

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

3. Mike said

@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”.

4. programmingpraxis said

@Mike: You are correct. Fixed. Thank you.

5. danaj said

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