Programming Praxis

The first if identifies obvious composites, and the second if does trial division to a million, which will be faster than anything else we can do. Then f-loop gathers factors until f exceeds the cube root of n, using both trial division and Pollard’s rho method, a-loop finds a suitable witness, and the last two clauses of the cond make the determination.

(define (pocklington? n)
  (if (not (pseudoprime? n)) #f
    (if (< n #e1e12) (= (car (td n #e1e6)) n)
      (let* ((n-1 (- n 1))
             (root2 (iroot 2 n)) (root3 (iroot 3 n))
             (fs (td n-1 1000)) (f (apply * fs)))
        (let f-loop ((f f) (fs fs))
          (if (< f root3)
              (let ((p (rho (/ n-1 f) #e1e6)))
                (if p (f-loop (* f p) (cons p fs)) 'unknown))
              (let a-loop ((a 2))
                (if (not (= (expm a n-1 n) 1)) #f
                  (let loop ((fs fs))
                    (let ((g (gcd (- (expm a (/ n-1 (car fs)) n) 1) n)))
                      (cond ((< 1 g n) #f)
                            ((= g n) (a-loop (+ a 1)))
                            ((pair? (cdr fs)) (loop (cdr fs)))
                            ((< root2 f) #t)
                            (else (let* ((ds (digits n f))
                                         (c2 (car ds)) (c1 (cadr ds))
                                         (s (- (* c1 c1) (* 4 c2))))
                                    (not (square? s)))))))))))))))

Here are some examples. Note the test of 3825123056546413051, which is a strong pseudoprime to bases 2, 3 and 5, but is correctly identified as composite by the Pocklington test. By the way, the symbol unknown is considered to be true by Scheme, which you will want to consider as you are checking return values from the function:

> (pocklington? (mersenne 89))
#t
> (pocklington? (mersenne 521))
#t
> (pocklington? (mersenne 607))
unknown
> (pseudoprime? 3825123056546413051)
#t
> (pocklington? 3825123056546413051)
#f

We stole much code from the Standard Prelude and previous exercises, too much to make a list here, but you can see it all assembled at http://programmingpraxis.codepad.org/UB4OxZv1.

Pages: 1 2