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

**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.

Solution

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
+----+-----+-------+
```

- Reverse the digits of a number using a list
- Exercise 12.10 from the book Scheme and the art of programming
- How to create Alist from List using syntax-rules in Scheme?
- Anonymous lambdas directly referring to themselves
- Why can't a Scheme macro with the name "if" be defined?
- How to implement "if" in normal order in Scheme?
- Map-reduce functional outline
- How to use let-values in Gambit Scheme?
- Execute command line from Scheme (Guile)
- How do I define a sub environment in scheme?
- Sources for learning about Scheme Macros: define-syntax and syntax-rules
- What is "Call By Name"?
- DrRacket's call/cc
- Racket: Conditional Statements Not Matching String Equality in Gambling Simulation"
- What's wrong with this coin change python implementation?
- Encode "ä", "ö", "ü" and "ß" in scheme (get german text from file)
- Using AND with the apply function in Scheme
- Combination of list-ref and index-of in racket
- How to undefine a variable in Scheme?
- Running Scheme from the command line
- Function that splits word
- How does the yin-yang puzzle work?
- Dr Racket Recursing Without Returning Inital Parent Node Within Function
- Chez Scheme FFI Procedure Doesn't Work After Change to Apple Silicon
- How do we convert this Scheme (Lisp) function into C#
- What is #~ in scheme?
- Racket : Run scheme file in terminal and show result
- Recursive numeric equality in Scheme
- Under any Scheme standard, is (let (x y z) x) valid code?
- Iterate through the letters of the alphabet in Racket