SICP Exercise 1.28 - Miller-Rabin - "at least half the numbers will reveal a nontrivial square root of 1 modulo n"

Question

SICP Exercise 1.28

https://mitpress.mit.edu/sites/default/files/sicp/full-text/book/book-Z-H-11.html#%_thm_1.28

One variant of the Fermat test that cannot be fooled is called the Miller-Rabin test (Miller 1976; Rabin 1980). This starts from an alternate form of Fermat's Little Theorem, which states that if n is a prime number and a is any positive integer less than n, then a raised to the (n - 1)st power is congruent to 1 modulo n. To test the primality of a number n by the Miller-Rabin test, we pick a random number a < n and raise a to the (n - 1)st power modulo n using the expmod procedure. However, whenever we perform the squaring step in expmod, we check to see if we have discovered a ``nontrivial square root of 1 modulo n,'' that is, a number not equal to 1 or n - 1 whose square is equal to 1 modulo n. It is possible to prove that if such a nontrivial square root of 1 exists, then n is not prime. It is also possible to prove that if n is an odd number that is not prime, then, for at least half the numbers a < n, computing a^(n-1) in this way will reveal a nontrivial square root of 1 modulo n. (This is why the Miller-Rabin test cannot be fooled.) Modify the expmod procedure to signal if it discovers a nontrivial square root of 1, and use this to implement the Miller-Rabin test with a procedure analogous to fermat-test. Check your procedure by testing various known primes and non-primes. Hint: One convenient way to make expmod signal is to have it return 0.

I've written my own solution and its results are consistent with the solutions provided here:

http://community.schemewiki.org/?sicp-ex-1.28

15 is an odd number that is not prime, so for at least half the numbers a from 1 to 14, I expect computing expmod(a, 14, 15) will reveal a nontrivial square root of 1 modulo n, which is signified by expmod returning 0.

However, these are the results I get:

(expmod 1 14 15)
> 1
(expmod 2 14 15)
> 4
(expmod 3 14 15)
> 9
(expmod 4 14 15)
> 0
(expmod 5 14 15)
> 10
(expmod 6 14 15)
> 6
(expmod 7 14 15)
> 4
(expmod 8 14 15)
> 4
(expmod 9 14 15)
> 6
(expmod 10 14 15)
> 10
(expmod 11 14 15)
> 0
(expmod 12 14 15)
> 9
(expmod 13 14 15)
> 4
(expmod 14 14 15)
> 1

As can be seen, only 2 of these results are 0, which is way short of at least 7 as expected.

Am I misunderstanding the statement? Am I being a complete idiot? Is the code wrong? Is SICP wrong? Many thanks.

Edit 1: it was requested that I supply the exact code I'm using. Here it is, although I'm essentially just copying the solution I linked to, and aliasing remainder as mod because that's what my interpreter calls it.

 (define (square x) (* x x))

 (define remainder mod)

 (define (miller-rabin-expmod base exp m) 
   (define (squaremod-with-check x) 
     (define (check-nontrivial-sqrt1 x square) 
       (if (and (= square 1) 
                (not (= x 1)) 
                (not (= x (- m 1)))) 
           0 
           square)) 
     (check-nontrivial-sqrt1 x (remainder (square x) m))) 
   (cond ((= exp 0) 1) 
         ((even? exp) (squaremod-with-check 
                       (miller-rabin-expmod base (/ exp 2) m))) 
         (else 
          (remainder (* base (miller-rabin-expmod base (- exp 1) m)) 
                     m))))

(define expmod miller-rabin-expmod)

(print (expmod 1 14 15))
(print (expmod 2 14 15))
(print (expmod 3 14 15))
(print (expmod 4 14 15))
(print (expmod 5 14 15))
(print (expmod 6 14 15))
(print (expmod 7 14 15))
(print (expmod 8 14 15))
(print (expmod 9 14 15))
(print (expmod 10 14 15))
(print (expmod 11 14 15))
(print (expmod 12 14 15))
(print (expmod 13 14 15))
(print (expmod 14 14 15))

Edit 2: I have now also manually calculated the steps of expmod(a, 14, 15) (which always recurses via exp = 14, exp = 7, exp = 6, exp = 3, exp = 2, exp = 1, exp = 0), for all values of a from 1 to 14, and I'm certain that only a = 4 and a = 11 encounter a nontrivial square root of 1. So I'm inclined to think that SICP is either wrong about this, or is not expressing itself clearly.

Answer 1

SICP is wrong because it uses the wrong definition of a Miller–Rabin witness (see Keith Conrad, The Miller–Rabin Test). In particular, the following line is wrong:

Wrong claim. It is also possible to prove that if n is an odd number that is not prime, then, for at least half the numbers a < n, computing a^{n - 1} in this way will reveal a nontrivial square root of 1 modulo n.

You can verify this to be false, an example is when n = 9.

Following Theorem 2.9 in the above reference, the correct statement should be as such:

Correct claim. Let n > 1 be an odd number that is not prime. Then we can write n - 1 as (2^e)k such that e &geq; 1 and k is odd. (For example, if n = 21, we can write 21 - 1 = 20 = (2^2)·5, so that e = 2 ≥ 1 and k = 5 is odd.) It is possible to prove that for at least half the numbers a < n, we will have a^k ≢ 1 and a^k ≢ n - 1 and a^{2k} ≢ n - 1 and ... and a^{2^{e-1}k} ≢ n - 1 modulo n.

So, for n = 21, we can prove that for at least half the numbers a < 21, we will have a^5 ≢ 1 and a^5 ≢ 20 and a^10 ≢ 20. We get the following table (with all values given modulo 21):

+----+-----+-------+
| a  | a^5 |  a^10 | MILLER–RABIN WITNESS? 
+----+-----+-------+
|  1 |   1 |     1 | NO, a^5 ≡ 1
|  2 |  11 |    16 | YES
|  3 |  12 |    18 | YES
|  4 |  16 |     4 | YES
|  5 |  17 |    16 | YES
|  6 |   6 |    15 | YES
|  7 |   7 |     7 | YES
|  8 |   1 |     1 | NO, a^5 ≡ 1
|  9 |  18 |     9 | YES
| 10 |  19 |     4 | YES
| 11 |   2 |     4 | YES
| 12 |   3 |     9 | YES
| 13 |  13 |     1 | YES
| 14 |  14 |     7 | YES
| 15 |  15 |    15 | YES
| 16 |   4 |    16 | YES
| 17 |   5 |     4 | YES
| 18 |   9 |    18 | YES
| 19 |  10 |    16 | YES
| 20 |  20 |     1 | NO, a^5 ≡ 20
+----+-----+-------+

And sure enough, more than half of the a < 21 (in fact, more than 75%) satisfy all three congruences a^5 ≢ 1, a^5 ≢ 20 and a^10 ≢ 20. (We call them Miller–Rabin witnesses; because they witness the fact that n is not prime. In general, many primality tests rely on some property that holds for all prime numbers—if you show that such a property fails for some number, then that number cannot be prime. The more witnesses, the better the primality test works.)

EDIT. Just as an example to illustrate primality, here is the table for n = 13. Naturally there cannot be any Miller–Rabin witnesses for 13 as it is prime; there is no non-primality to witness. Since n = 13, we have n - 1 = 12 = (2^2)·3, so that e = 2 ≥ 1 and k = 3 is odd. Thus, as argued in page 1 of Keith Conrad's expository paper, all a < 13 will satisfy at least one of the three congruences a^3 ≡ 1, a^3 ≡ 12, a^6 ≡ 12. And sure enough:

+----+-----+-------+
| a  | a^3 |   a^6 | MILLER–RABIN WITNESS? 
+----+-----+-------+
|  1 |   1 |     1 | NO, a^3 ≡ 1
|  2 |   8 |    12 | NO, a^6 ≡ 12
|  3 |   1 |     1 | NO, a^3 ≡ 1
|  4 |  12 |     1 | NO, a^3 ≡ 12
|  5 |   8 |    12 | NO, a^6 ≡ 12
|  6 |   8 |    12 | NO, a^6 ≡ 12
|  7 |   5 |    12 | NO, a^6 ≡ 12
|  8 |   5 |    12 | NO, a^6 ≡ 12
|  9 |   1 |     1 | NO, a^3 ≡ 1
| 10 |  12 |     1 | NO, a^3 ≡ 12
| 11 |   5 |    12 | NO, a^6 ≡ 12
| 12 |  12 |     1 | NO, a^3 ≡ 12
+----+-----+-------+

Answer 2

I found a paper that covers this test and a proof of the specific result that over half of the values between 2 and n-2 will result in a nontrivial square root of 1. (Theorem 4.1)

• https://kconrad.math.uconn.edu/blurbs/ugradnumthy/millerrabin.pdf

I made this code to double check it

(define (print x) (display x) (newline))

(define (assert p) (unless p (error 'assert-failed)))

(define (power-of-two-split m)
  ;; write m = 2^e k
  (let loop ((e 0) (k m))
    (if (even? k)
        (loop (+ e 1) (quotient k 2))
        (cons e k))))

(define (exp/mod a k n)
  ;; a^k (mod n)
  (cond ((= k 0) 1)
        ((= k 1) (modulo a n))
        (else (modulo (* a (exp/mod a (- k 1) n)) n))))

(define (miller-rabin a n)
  (assert (odd? n))
  (assert (= 3 (modulo n 4))) ;; only handles e=1 case, need to use power-of-two-split for full test
  (let ((k (quotient (- n 1) 2)))
    (exp/mod a k n)))

(define (test n)
  (for ((i (in-range 2 (- n 2))))
    (let ((m (miller-rabin i n)))
      (print `(,i -> ,m squared ,(exp/mod m 2 n))))))

(test 15)

it prints the following result

(2 -> 8 squared 4)
(3 -> 12 squared 9)
(4 -> 4 squared 1)
(5 -> 5 squared 10)
(6 -> 6 squared 6)
(7 -> 13 squared 4)
(8 -> 2 squared 4)
(9 -> 9 squared 6)
(10 -> 10 squared 10)
(11 -> 11 squared 1)
(12 -> 3 squared 9)

So checking against the formal definition of a miller-rabin witness, they are actually all wintesses:

Definition 2.3. For odd n > 1, write n − 1 = 2^ek with k odd and pick a ∈ {1, . . . , n − 1}. We say a is a Miller–Rabin witness for n if all of the congruences in are false: * a^k = 1 mod n and a^{(2^i)k} = −1 mod n for all i ∈ {0, . . . , e − 1}.

and you can see none of the values in the 'm' column are 1 and none of the values in the squared column are 14. So they are all witnesses, so >50% are.

The code that's doing "check-nontrivial-sqrt1" isn't relevant in this particular case of n=3 (mod 4) because e=1 in that case.

---

Update:

I just realized the reason why we have many witnesses but we don't always find a square root from them all:

the idea behind Miller–Rabin witnesses is to find an unexpected square root of 1 mod n. This is not always what we actually find, since the premise a n−1 ≡ 1 mod n for prime n might not actually be true for composite n.

here's a table of a^(n-1) (mod n) for n=15

(2 -> 4)
(3 -> 9)
(4 -> 1)
(5 -> 10)
(6 -> 6)
(7 -> 4)
(8 -> 4)
(9 -> 6)
(10 -> 10)
(11 -> 1)
(12 -> 9)

as you can see, only two times does the congruence a^(n-1) = 1 (mod n) actually hold.

Answer 3

From @Memes's answer, I've gone ahead and added scheme code for it:

(define (display-all . vs)
  (for-each display vs))

(define (find-e-k n)
  (define (find-e-k-iter possible-k possible-e)
    (if (= (remainder possible-k 2) 0)
        (find-e-k-iter (/ possible-k 2) (+ possible-e 1))
        (values possible-e possible-k)))
  (find-e-k-iter (- n 1) 0))

; first-witness-case-test: (a ^ k) mod n # 1
(define (first-witness-case-test a k n)
  (not (= (expmod a k n) 1)))

; second-witness-case-test: all a ^ ((2 ^ i) * k) (with i = {0..e-1}) mod n # (n - 1)
(define (second-witness-case-test a e k n)
  (define (second-witness-case-test-iter a i k n)
    (cond ((= i -1) true)
          (else (let ()
                 (define witness (not (= (expmod a (* (fast-expt 2 i) k) n) (- n 1))))
                 (if witness
                 (second-witness-case-test-iter a (- i 1) k n)
                 false)))))
  (second-witness-case-test-iter a (- e 1) k n))

(define (miller-rabin-test n)
  (define (try-it a e k)
    (if (and (first-witness-case-test a k n) (second-witness-case-test a e k n))
        (display-all "is not prime, with a = " a "\n")
        (if (< a (- n 1))
            (try-it (+ a 1) e k)
            (display "is prime\n"))))
  (cond ((< n 2) (display "not prime"))
        ((= (remainder n 2) 0) (display "not prime\n"))
        (else (let ()
               (define-values (e k) (find-e-k n))
               (try-it 1 e k)))))

Answer 4

One variant of the Fermat test that cannot be fooled is called the Miller-Rabin test (Miller 1976; Rabin 1980). This starts from an alternate form of Fermat's Little Theorem, which states that if n is a prime number and a is any positive integer less than n, then a raised to the (n - 1)st power is congruent to 1 modulo n.

• n is prime, a < n ==> a^(n-1) = 1 (mod n)

To test the primality of a number n by the Miller-Rabin test, we pick a random number a < n and raise a to the (n - 1)st power modulo n using the expmod procedure.

So pick a randomly and check if a^(n-1) = 1 (mod n). if it's not then you know n is not prime.

However, whenever we perform the squaring step in expmod, we check to see if we have discovered a ``nontrivial square root of 1 modulo n,'' that is, a number not equal to 1 or n - 1 whose square is equal to 1 modulo n.

This is talking about an extra check being added inside the expmod function. This might be the thing you overlooked.

It is possible to prove that if such a nontrivial square root of 1 exists, then n is not prime.

Let's go over this in detail. A nontrivial square root of 1 would be a number x such that x^2 = 1 (mod n). and x is not 1 or -1.

Why would one of these indicate n not being prime?

We know that x^2 - 1 = (x - 1)(x + 1) and (working modulo n) both x-1 and x+1 are not zero, but their product is zero. This implies we have a composite modulus and you can factor it by taking the GCD of those two values.

It is also possible to prove that if n is an odd number that is not prime, then, for at least half the numbers a < n, computing a^(n-1) in this way will reveal a nontrivial square root of 1 modulo n. (This is why the Miller-Rabin test cannot be fooled.)

This again is talking about adding an internal test to the squaring branch of the expmod function.

Modify the expmod procedure to signal if it discovers a nontrivial square root of 1, and use this to implement the Miller-Rabin test with a procedure analogous to fermat-test. Check your procedure by testing various known primes and non-primes. Hint: One convenient way to make expmod signal is to have it return 0.

Hope that helps! Ask if you need any more guidance.

Answer 5

From [Cormen et al. 2009], p. 955:

Theorem 31.31 (Fermat’s theorem)
If p is prime, then a^{p − 1} ≡ 1 (mod p) for all a ∈ Z_p*.

From [Cormen et al. 2009], p. 956:

Theorem 31.34
If p is an odd prime and e ≥ 1, then the equation x² ≡ 1 (mod p^e) has only two solutions, namely x = 1 and x = −1.

[…]

A number x is a nontrivial square root of 1, modulo n, if it satisfies the equation x² ≡ 1 (mod n) but x is equivalent to neither of the two “trivial” square roots: 1 or −1, modulo n. For example, 6 is a nontrivial square root of 1, modulo 35. We shall use the following corollary to Theorem 31.34 in the correctness proof in Section 31.8 for the Miller-Rabin primality-testing procedure.

Corollary 31.35
If there exists a nontrivial square root of 1, modulo n, then n is composite.

From [Cormen et al. 2009], p. 971:

Theorem 31.38
If n is an odd composite number, then the number of witnesses to the compositeness of n is at least (n − 1)/2.

Incorrect claim in SICP:

It is also possible to prove that if n is an odd number that is not prime, then, for at least half the numbers a < n, computing a^(n-1) in this way will reveal a nontrivial square root of 1 modulo n.

Corrected claim (equivalent to theorem 31.38 in [Cormen et al. 2009]):

It is also possible to prove that if n is an odd number that is not prime, then, for at least half the numbers a < n, computing a^(n-1) in this way will reveal a nontrivial square root of 1 modulo n or a^(n-1) will not be equal to 1 module n.

SICP is wrong because a is a witness to the compositeness of n does not mean that there exists a^k a nontrivial square root of 1 modulo n, but it means that there exists a^k a nontrivial square root of 1 modulo n or aⁿ⁻¹ ≢ 1 (mod n). The reason for the second condition is that if there exists a^k a nontrivial square root of 1 modulo n then aⁿ⁻¹ ≡ 1 (mod n), and by contrapositive of theorem 31.31 in [Cormen et al. 2009], if aⁿ⁻¹ ≢ 1 (mod n) then n is composite.

Example. — For n = 9, we expect at least (n − 1)/2 = 4 witnesses to the compositeness of n. The expmod procedure reveals 6 witnesses:

(display (expmod 1 8 9)) (newline)  ; 1
(display (expmod 2 8 9)) (newline)  ; 4 (witness: a^(n-1) ≢ 1 (mod n))
(display (expmod 3 8 9)) (newline)  ; 0 (witness: nontrivial square root)
(display (expmod 4 8 9)) (newline)  ; 7 (witness: a^(n-1) ≢ 1 (mod n))
(display (expmod 5 8 9)) (newline)  ; 7 (witness: a^(n-1) ≢ 1 (mod n))
(display (expmod 6 8 9)) (newline)  ; 0 (witness: nontrivial square root)
(display (expmod 7 8 9)) (newline)  ; 4 (witness: a^(n-1) ≢ 1 (mod n))
(display (expmod 8 8 9)) (newline)  ; 1

Reference

Cormen (Thomas), Leiserson (Charles), Rivest (Ronald), Stein (Clifford), Introduction to Algorithms, 3rd ed., MIT Press, Cambridge, Massachusetts, 2009.