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.
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 ≥ 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
+----+-----+-------+