Lucas Pseudoprimes

February 4, 2014

We begin with Selfridge’s procedure for choosing P and Q, which returns P = Q = 0 if d is a factor of n:

(define (selfridge n)
  (let loop ((d-abs 5) (sign 1))
    (let ((d (* d-abs sign)))
      (cond ((< 1 (gcd d n)) (values d 0 0))
            ((= (jacobi d n) -1) (values d 1 (/ (- 1 d) 4)))
            (else (loop (+ d-abs 2) (- sign)))))))

Then the standard Lucas test computes the U sequence:

(define (standard-lucas-pseudoprime? n)
  ; assumes odd n > 1 not a square
  (call-with-values
    (lambda () (selfridge n))
    (lambda (d p q)
      (if (zero? p) (= n d)
        (call-with-values
          (lambda () (lucas p q n (+ n 1)))
          (lambda (u v k) (zero? u)))))))

The strong Lucas test uses an auxiliary function to remove the factors of 2 from n + 1:

(define (powers-of-two n)
  (let loop ((s 0) (n n))
    (if (odd? n) (values s n)
      (loop (+ s 1) (/ n 2)))))

(define (strong-lucas-pseudoprime? n)
  ; assumes odd positive integer not a square
  (call-with-values
    (lambda () (selfridge n))
    (lambda (d p q)
      (if (zero? p) (= n d)
        (call-with-values
          (lambda () (powers-of-two (+ n 1)))
          (lambda (s t)
            (call-with-values
              (lambda () (lucas p q n t))
              (lambda (u v k)
                (if (or (zero? u) (zero? v)) #t
                  (let loop ((r 1) (v v) (k k))
                    (if (= r s) #f
                      (let* ((v (modulo (- (* v v) (* 2 k)) n))
                             (k (modulo (* k k) n)))
                        (if (zero? v) #t (loop (+ r 1) v k))))))))))))))

Then the Baillie-Wagstaff pseudoprimality test works like this, where we have added a strong pseudoprime test to base 3:

(define prime?
  (let ((ps (primes 100)))
    (lambda (n)
      (if (not (integer? n)) (error 'prime? "must be integer"))
      (if (or (< n 2) (square? n)) #f
        (let loop ((ps ps))
          (if (pair? ps)
              (if (zero? (modulo n (car ps))) (= n (car ps)) (loop (cdr ps)))
              (and (strong-pseudoprime? n 2)
                   (strong-pseudoprime? n 3)
                   (strong-lucas-pseudoprime? n))))))))

Since we got it wrong twice before, it is important to test our new functions. At oeis.org, sequence A217120 has the standard Lucas pseudoprimes and sequence A217255 has the strong Lucas pseudoprimes:

> (let ((primes48000 (primes 48000)))
    (list-of n ; standard lucas pseudoprimes A217120
      (n range 3 48000 2)
      (not (square? n))
      (standard-lucas-pseudoprime? n)
      (not (member n primes48000))))
  (323 377 1159 1829 3827 5459 5777 9071 9179 10877 11419 11663
   13919 14839 16109 16211 18407 18971 19043 22499 23407 24569
   25199 25877 26069 27323 32759 34943 35207 39059 39203 39689
   40309 44099 46979 47879)

> (let ((primes325000 (primes 325000)))
    (list-of n ; strong lucas pseudoprimes A217255
      (n range 3 325000 2)
      (not (square? n))
      (strong-lucas-pseudoprime? n)
      (not (member n primes325000))))
  (5459 5777 10877 16109 18971 22499 24569 40309 58519 75077
   97439 113573 115639 130139 155819 158399 161027 162133
   176399 176471 189419 192509 197801 224369 230691 243629
   253259 268349 288919 313499 324899)

To test the prime? function, we tested all numbers in the ranges x to x + 100000 for x ∈ {103, 106, 109, 1012}, comparing to a list of prime numbers generated by a segmented sieve. We also test all 10000 base-2 and base-3 pseudoprimes from A072276, making sure they are reported as composite:

(do ((xs (list #e1e3 #e1e6 #e1e9 #e1e12) (cdr xs))) ((null? xs))
  (let loop ((x (car xs)) (ps (primes (car xs) (+ (car xs) 100000))))
    (when (pair? ps)
      (when (and (= x (car ps)) (not (prime? x)))
        (display "prime reported as composite ") (display x) (newline))
      (when (and (prime? x) (not (= x (car ps))))
        (display "composite reported as prime ") (display x) (newline))
      (loop (+ x 1) (if (= x (car ps)) (cdr ps) ps)))))

(with-input-from-file "23primes" (lambda ()
  (do ((x (read) (read))) ((eof-object? x))
    (when (prime? x) (display x) (newline)))))

We reused much code from previous exercises, including the strong pseudoprime test from the Miller-Rabin primality checker. Though we don’t show them here, we did similar tests comparing strong pseudoprimes base 2 to A001262 and strong pseudoprimes base 3 to A020229. None of the tests showed any errors. You can see and run the complete program at http://programmingpraxis.codepad.org/nb2VLgXb.

About these ads

Pages: 1 2

One Response to “Lucas Pseudoprimes”

  1. danaj said

    I believe your test description came from Nicely’s site rather than the source material. Baillie and Wagstaff, “Lucas Pseudoprimes”, October 1980. Page 1401. Pomerance, Selfridge, Wagstaff, “The Pseudoprimes to 25 x 10^9″, July 1980. Page 1024-1025. Pomerance, “Are there counter-examples to the Baillie-PSW Primality Test?”, 1984.

    Baillie and Wagstaff page 1401: 1. Trial divide to some convenient limit (1000 is given as an example, but this is is just a practical step, not a requirement of the test). 2. if n is not a strong probable prime to base 2, then n is composite. 3. Choose P & Q by method A or B (A is the Selfridge method). 4. If n is not a strong Lucas probable prime with parameters P and Q, then n is a composite.

    Step 1 is left out in PSW80 and P84 — it’s an obvious practical step, but not necessary. Step 2 is is “strong probable prime to base 2″ in all papers. Some later authors weaken this (e.g. Chen and Greene 2003) but there is no reason to prefer the standard test to the strong test. Step 3 is selection of parameters, which P84 limits to just the Selfridge method, but the other two allow both methods (though I’ve not seen anyone use method B). Step 4 is the strong Lucas test in BW80. It is just the standard version in PSW80 (which doesn’t mention strong Lucas pseudoprimes at all) and P84.

    Lots of software take different steps and still call it a BPSW test. Step 1 is just a pretest so given software can do whatever it thinks makes things faster. If one is creating a primality test now, you would want to use at leas the strong Lucas test — it runs as fast or faster than the standard test and the pseudoprimes are a subset. Pari, for example, uses a variant of the extra strong test. I prefer the full extra strong test using the Baillie parameters (A217719) since it has fewer pseudoprimes, runs faster than the standard or strong test, and maintains the same residue class behavior.

    (caveat: link to my page) A comparison of some of the tests and the number of pseudoprimes can be seen here: http://www.sti15.com/nt/pseudoprimes.html
    An incomplete list of various software that that performs the BPSW test and the specific tests they use can be found on the Wikipedia BPSW talk page.

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 )

Twitter picture

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

Facebook photo

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

Google+ photo

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

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 630 other followers

%d bloggers like this: