## Two Powering Predicates

### August 6, 2010

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 whethernis of the formpwith^{k}pprime, and if it is outputs the primep.1. [Initialize] Set

a← 1.2. [Compute GCD] Set

a=a+ 1 and computeb←amod^{n}nusing one of the powering algorithms in the ring Z/nZ, and then compute the GCDp← (b−a,n).3. [Finished?] If

p= 1,nis not a prime power, and terminate the algorithm. Otherwise, using a compositeness test, test whetherpis composite. If it is, go to Step 2.4. [Final test] (Here

pis almost certainly prime.) Using a primality test prove thatpis prime. If it is not (an exceedingly rare occurrence), go to Step 2. Otherwise, by dividingnbyprepeatedly, check whethernis a power ofpor not. If it is not,nis not a prime power and terminate. If it is, outputpand 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 replacingnbyp. This occurs rarely however and so it is not necessary to bother.Second, to check whether

nis a power ofp, instead of repeatedly dividing byp, 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*^{1} 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.

Something went wrong during editing of the exercise, and the code given for the

`square?`

function was incorrect. It has now been fixed.[…] method, then multiply to determine if the original number is a square. But that’s slow. In a previous exercise, we used a method devised by Henri Cohen to calculate the quadratic residues of n to various […]