Programming Praxis

There is a lot of code common to the two versions of the Lucas pseudoprimality test, so we will abstract much of the code to auxiliary functions. Here is the code to implement Selfridge’s selection of P and Q:

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

The biggest part of the code computes the Lucas chain. We compute and return a value qkd which is used only in the strong version of the test so that we can share the code between the two tests. It’s imperative to compute the various intermediate values in the correct order; I made more than one mistake with out-of-order computations during the development of this code. For those not familiar with Scheme, the let* version of let computes its results left to right, with each computed value available to subsequent computations:

(define (chain n u v u2 v2 d q m)
  (let loop ((u u) (v v) (u2 u2) (v2 v2) (q q) (qkd q) (m m))
    (if (zero? m) (values u v qkd)
      (let* ((u2 (modulo (* u2 v2) n))
             (v2 (modulo (- (* v2 v2) (* 2 q)) n))
             (q (modulo (* q q) n)))
        (if (odd? m)
            (let* ((t1 (* u2 v)) (t2 (* u v2))
                   (t3 (* v2 v)) (t4 (* u2 u d))
                   (u (+ t1 t2)) (v (+ t3 t4))
                   (qkd (modulo (* qkd q) n)))
              (loop (modulo (quotient (if (odd? u) (+ u n) u) 2) n)
                    (modulo (quotient (if (odd? v) (+ v n) v) 2) n)
                    u2 v2 q qkd (quotient m 2)))
            (loop u v u2 v2 q qkd (quotient m 2)))))))

With this, the standard version of the test is quite simple. The m value in the chain is (n + 1) / 2, and the only return value needed from chain is u:

(define (standard-lucas-pseudoprime? n)
  (call-with-values
    (lambda () (selfridge n))
    (lambda (d p q)
      (if (zero? p) (= n d)
        (call-with-values
          (lambda () (chain n 0 2 1 p d q (/ (+ n 1) 2)))
          (lambda (u v qkd) (zero? u)))))))

For the strong version of the test, we need to split n + 1 into its constituent powers of two:

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

Most of the strong test is reused from the standard test. The m value in the chain is ⌊t / 2⌋. The primary difference comes after the chain is computed:

(define (strong-lucas-pseudoprime? n)
  (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 () (chain n 1 p 1 p d q (quotient t 2)))
              (lambda (u v qkd)
                (if (or (zero? u) (zero? v)) #t
                    (let loop ((r 1) (v v) (qkd qkd))
                      (if (= r s) #f
                        (let* ((v (modulo (- (* v v) (* 2 qkd)) n))
                             (qkd (modulo (* qkd qkd) n)))
                        (if (zero? v) #t (loop (+ r 1) v qkd))))))))))))))

With all that, the prime? function is straight forward. We choose the primes less than a hundred for the trial division step of the algorithm, thougn many implementations prefer a limit of 1000; since ps is bound inside the define but outside the lambda, it is computed once, at compile time, and saved from one call of the function to the next. Here’s prime?:

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

We give several examples below. The first four are small calculations useful for testing; you can augment this implementation to print the results of intermediate calculations in strategic places, using them to debug your code. Then 323 and 5777 are, respectively, the first standard Lucas pseudoprime and the first strong Lucas pseudoprime, as you can see from the displayed results. The number that identified the error with the implementation in the previous exercise, 3825123056546413051, is tested next, so you can see that the function now works properly. The final test, 2⁸⁹ − 1, is a moderate-sized prime:

> (do ((ns (list 79 83 89 111 323 5777 3825123056546413051 (- (expt 2 89) 1)) (cdr ns))) ((null? ns))
    (let ((n (car ns)))
      (display n) (display " ")
      (display (strong-pseudoprime? n 2)) (display " ")
      (display (strong-pseudoprime? n 3)) (display " ")
      (display (standard-lucas-pseudoprime? n)) (display " ")
      (display (strong-lucas-pseudoprime? n)) (display " ")
      (display (prime? n)) (newline)))
83 #t #t #t #t #t
89 #t #t #t #t #t
111 #f #f #f #f #f
323 #f #f #t #f #f
5777 #f #f #t #t #f
3825123056546413051 #t #t #f #f #f
618970019642690137449562111 #t #t #t #t #t

We reused much code from several previous exercises, which you can see assembled at http://programmingpraxis.codepad.org/bHofkLUf. A Python version is also available.

Pages: 1 2