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.

Advertisements

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

Leave a Reply

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

WordPress.com Logo

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

Google photo

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

Twitter picture

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

Facebook photo

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

Connecting to %s

%d bloggers like this: