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]

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.


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

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

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: