Programming Praxis

We begin with the algorithm to compute integer square roots by Newton’s method. The Standard Prelude already had an implementation of that function, but Cohen’s version is both prettier and faster, and generates less garbage:

(define (isqrt n)
  (if (not (and (positive? n) (integer? n)))
      (error 'isqrt "must be positive integer")
      (let loop ((x n))
        (let ((y (quotient (+ x (quotient n x)) 2)))
          (if (< y x) (loop y) x)))))

We will replace the existing implementation of the function in the Standard Prelude with the one shown above.

Here is our version of Cohen’s square? function. The pre-computations are done once when the function is loaded, not when it is called, and are visible only to the square? function. As Cohen notes, less than one percent of numbers survive the four table lookups, so the function is quite fast, an order of magnitude faster than the naive calculation of the square of the integer square root:

(define square?
  (let ((q11 (make-vector 11 #f))
        (q63 (make-vector 63 #f))
        (q64 (make-vector 64 #f))
        (q65 (make-vector 65 #f)))
    (do ((k 0 (+ k 1))) ((< 5 k))
      (vector-set! q11 (modulo (* k k) 11) #t))
    (do ((k 0 (+ k 1))) ((< 31 k))
      (vector-set! q63 (modulo (* k k) 63) #t))
    (do ((k 0 (+ k 1))) ((< 31 k))
      (vector-set! q64 (modulo (* k k) 64) #t))
    (do ((k 0 (+ k 1))) ((< 32 k))
      (vector-set! q65 (modulo (* k k) 65) #t))
    (lambda (n)
      (if (not (integer? n)) (error 'square? "must be integer")
        (if (< n 1) #f
          (if (not (vector-ref q64 (modulo n 64))) #f
            (let ((r (modulo n 45045)))
              (if (not (vector-ref q63 (modulo r 63))) #f
                (if (not (vector-ref q65 (modulo r 65))) #f
                  (if (not (vector-ref q11 (modulo r 11))) #f
                    (let ((q (isqrt n)))
                      (if (= (* q q) n) q #f))))))))))))

The implementation of the prime-power checker involves an interesting story. The version of Cohen’s book that I checked out of the library gave this algorithm, which differs from the one above:

Given a positive integer n, this algorithm tests whether n is of the form p^k with p prime, and if it is outputs the prime p.

1. [Initialize] Set a ← 1.

2. [Compute GCD] Set a = a + 1 and compute b ← aⁿ mod n using one of the powering algorithms in the ring Z/nZ, and then compute the GCD p ← (b−a, n).

3. [Finished?] If p = 1, n is not a prime power, and terminate the algorithm. Otherwise, using a compositeness test, test whether p is composite. If it is, go to Step 2.

4. [Final test] (Here p is almost certainly prime.) Using a primality test prove that p is prime. If it is not (an exceedingly rare occurrence), go to Step 2. Otherwise, by dividing n by p repeatedly, check whether n is a power of p or not. If it is not, n is not a prime power and terminate. If it is, output p and terminate the algorithm.

We have been a little sloppy in this algorithm. First, when one finds that p > 1 is composite in Step 3, we could use this information for the next trial. In fact, we could apply the algorithm recursively by replacing n by p. This occurs rarely however and so it is not necessary to bother.

Second, to check whether n is a power of p, instead of repeatedly dividing by p, one can use a binary search analogous to the binary powering algorithm. We leave this as an exercise for the reader.

I implemented that algorithm as shown below, fixing the second bit of sloppiness with the ilog function from the Standard Prelude:

(define (prime-power? n)
  (if (even? n) (if (= (expt 2 (ilog 2 n)) n) 2 #f)
    (let loop ((a 2))
      (let* ((b (expm a n n)) (p (gcd (- b a) n)))
        (cond ((= p 1) #f)
              ((prime? p) (if (= (expt p (ilog p n)) n) p #f))
              (else (loop (+ a 1))))))))

A few quick tests showed that it seemed to be working, so I launched a loop to compute the prime powers less than a thousand. It didn’t take long before I was stuck at 45 = 3 × 3 × 5. A bit of hand calculation quickly showed the problem, and a bit of head-scratching suggested that fixing the first bit of sloppiness was not optional, so I wrote a second version of the function:

(define (prime-power? n)
  (if (even? n) (if (= (expt 2 (ilog 2 n)) n) 2 #f)
    (let loop ((a 2) (n n))
      (let* ((b (expm a n n)) (p (gcd (- b a) n)))
        (cond ((= p 1) #f)
              ((prime? p) (if (= (expt p (ilog p n)) n) p #f))
              (else (loop (+ a 1) p)))))))

That fixed the problem at 45, and I happily re-launched my loop to compute the prime powers less than a thousand. It stopped at 561. It doesn’t take much knowledge of number theory to recognize 561 as the smallest Carmichael number, which defeats Fermat’s Little Theorem and thus breaks the algorithm. Again, the fix was simple; since all Carmichael numbers have three factors, they can’t be prime powers, and Carmichael numbers are easily recognized when a = b. Here’s the next version:

(define (prime-power? n)
  (if (even? n) (if (= (expt 2 (ilog 2 n)) n) 2 #f)
    (let loop ((a 2) (n n))
      (let* ((b (expm a n n)) (p (gcd (- b a) n)))
        (cond ((or (= p 1) (= a b)) #f)
              ((prime? p) (if (= (expt p (ilog p n)) n) p #f))
              (else (loop (+ a 1) p)))))))

That fixed all the Carmichael numbers, but when I compared my list to Sloane’s A0009061, I realized I was missing all the prime numbers p, because p¹ is certainly a prime power. So I wrote one more version:

(define (prime-power? n)
  (cond ((not (integer? n)) (error 'prime-power? "must be integer"))
        ((< n 1) #f) ((= n 1) 1) ((prime? n) n)
        ((even? n) (if (= (expt 2 (ilog 2 n)) n) 2 #f))
        (else (let loop ((a 2) (n n))
                (let* ((b (expm a n n)) (p (gcd (- b a) n)))
                  (cond ((or (= p 1) (= a b)) #f)
                        ((prime? p) (if (= (expt p (ilog p n)) n) p #f))
                        (else (loop (+ a 1) p))))))))

That version also forces (prime-power? 1) to be 1, which is consistent with Sloane, although some mathematicians disagree with that definition.

After all that fun, I went looking for the original problem. What happened is that Cohen published the book and my library bought a copy. But readers quickly pointed out the problem with the prime power test, and a second printing of the book fixed the error. You can find the errata for the book at ftp://megrez.math.u-bordeaux.fr/pub/cohenbook.

The lesson for the programmer is clear: Even if you trust the source of borrowed code, test it yourself, because once you incorporate it into your program you become responsible for it.

The prime-power? function uses ilog and expm from the Standard Prelude and prime? from the exercise on the Baillie-Wagstaff primality checker. You can run the code, including examples, at http://programmingpraxis.codepad.org/hcMUSvP6.

Pages: 1 2