# sicp-ex-1.28

Maggyero

QUESTION

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.

``` (import (scheme small))
(import (srfi 27))

(define (miller-rabin-test n rounds)
(define (try-it rounds)
(define a (+ 1 (random-integer (- n 1))))
(cond
((= rounds 0) #t)
((= (expmod a n n) a) (try-it (- rounds 1)))
(else #f)))
(define (expmod base e mod)
(cond
((= e 0) 1)
((even? e) (squaremod (expmod base (/ e 2) mod) mod))
(else (remainder (* base (expmod base (- e 1) mod)) mod))))
(define (squaremod x mod)
(define y (remainder (square x) mod))
(if (and (= y 1) (not (= x 1)) (not (= x (- mod 1))))
0
y))
(and (> n 1) (try-it rounds)))

(display (miller-rabin-test 2 10)) (newline)
(display (miller-rabin-test 3 10)) (newline)
(display (miller-rabin-test 5 10)) (newline)
(display (miller-rabin-test 7 10)) (newline)
(display (miller-rabin-test 0 10)) (newline)
(display (miller-rabin-test 1 10)) (newline)
(display (miller-rabin-test 4 10)) (newline)
(display (miller-rabin-test 6 10)) (newline)
(display (miller-rabin-test 8 10)) (newline)
(display (miller-rabin-test 9 10)) (newline)
(display (miller-rabin-test 561 10)) (newline)
(display (miller-rabin-test 1105 10)) (newline)
(display (miller-rabin-test 1729 10)) (newline)
(display (miller-rabin-test 2465 10)) (newline)
(display (miller-rabin-test 2821 10)) (newline)
(display (miller-rabin-test 6601 10)) (newline)
```

This exercise is based on sicp-ex-1.27.

```
;; ex 1.28

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

(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 (miller-rabin-test n)
(define (try-it a)
(define (check-it x)
(and (not (= x 0)) (= x 1)))
(check-it (miller-rabin-expmod a (- n 1) n)))
(try-it (+ 1 (random (- n 1)))))

(define (fast-prime? n times)
(cond ((= times 0) true)
((miller-rabin-test n) (fast-prime? n (- times 1)))
(else false)))

(define (prime? n)
; Perform the test how many times?
; Use 100 as an arbitrary value.
(fast-prime? n 100))

(define (report-prime n expected)
(define (report-result n result expected)
(newline)
(display n)
(display ": ")
(display result)
(display ": ")
(display (if (eq? result expected) "OK" "FOOLED")))
(report-result n (prime? n) expected))

(report-prime 2 true)
(report-prime 7 true)
(report-prime 13 true)
(report-prime 15 false)
(report-prime 37 true)
(report-prime 39 false)

(report-prime 561 false)  ; Carmichael number
(report-prime 1105 false) ; Carmichael number
(report-prime 1729 false) ; Carmichael number
(report-prime 2465 false) ; Carmichael number
(report-prime 2821 false) ; Carmichael number
(report-prime 6601 false) ; Carmichael number

```

Another solution that avoids nested functions. Similar to http://wiki.drewhess.com/wiki/SICP_exercise_1.28. However, this version does not intermingle the non-trivial-sqrt property check with returning the correct result in expmod. Instead, the non-trivial-sqrt property check is a dedicated function which returns true or false, as it should be. Only works for n>0. Uses the 'let' special form not introduced at this point of the book.

```
(define (miller-rabin n)
(miller-rabin-test (- n 1) n))

(define (miller-rabin-test a n)
(cond ((= a 0) true)
; expmod is congruent to 1 modulo n
((= (expmod a (- n 1) n) 1) (miller-rabin-test (- a 1) n))
(else false)))

(define (expmod base exp m)
(cond ((= exp 0) 1)
((even? exp)
(let ((x (expmod base (/ exp 2) m)))
(if (non-trivial-sqrt? x m) 0 (remainder (square x) m))))
(else
(remainder (* base (expmod base (- exp 1) m))
m))))

(define (non-trivial-sqrt? n m)
(cond ((= n 1) false)
((= n (- m 1)) false)
; book reads: whose square is equal to 1 modulo n
; however, what was meant is square is congruent 1 modulo n
(else (= (remainder (square n) m) 1))))

```

Another solution in GNU Guile:

``` ;;;; Under Creative Commons Attribution-ShareAlike 4.0
;;;; International. See

(define-module (net ricketyspace sicp one twentyeight)
#:use-module (srfi srfi-1)
#:export (miller-rabin-test prime? run-tests))

(define (sqmod x m)
"Return x^2 if `x^2 mod m` is not equal to `1 mod m`
and x != m - 1 and x != 1; 0 otherwise."
(let ((square (* x x)))
(cond ((and  (= (remainder square m) 1) ; 1 mod m = 1
(not (= x (1- m)))
(not (= x 1)))
0)
(else square))))

(define (expmod base exp m)
(cond ((= exp 0) 1)
((even? exp)
(remainder (sqmod (expmod base (/ exp 2) m) m)
m))
(else
(remainder (* base (expmod base (- exp 1) m))
m))))

(define (miller-rabin-test n)
(define (pass? a)
(= (expmod a (1- n) n) 1))
(fold (lambda (a p) (and (pass? a) p)) #t (iota (1- n) 1)))

(define (prime? n)
(if (miller-rabin-test n) #t #f))

;;; Tests

(define (carmichael-numbers-pass?)
"Return #t if the sample carmichael numbers are detected as non-prime."
(let ((numbers '(561 1105 1729 2465 2821 6601)))
(cons "carmichael-numbers-pass?"
(fold (lambda (n p) (and (not (prime? n)) p)) #t numbers))))

(define (prime-numbers-pass?)
"Return #t if the sample prime numbers are detected as prime"
(let ((numbers '(311 641 829 599 809 127 419 13 431 883)))
(cons "prime-numbers-pass?"
(fold (lambda (n p) (and (prime? n) p)) #t numbers))))

(define (even-numbers-pass?)
"Return #t if the sample even numbers are detected as non-prime"
(let ((numbers '(302 640 828 594 804 128 414 12 436 888)))
(cons "prime-numbers-pass?"
(fold (lambda (n p) (and (not (prime? n)) p)) #t numbers))))

(define (run-tests)
(map (lambda (test) (test)) (list carmichael-numbers-pass?
prime-numbers-pass?
even-numbers-pass?)))

;;; Guile REPL
;;;
;;; scheme@(guile-user)> ,use (net ricketyspace sicp one twentyeight)
;;; scheme@(guile-user)> (run-tests)
;;; \$18 = (("carmichael-numbers-pass?" . #t)
;;;        ("prime-numbers-pass?" . #t)
;;;        ("prime-numbers-pass?" . #t))

```

tiendo1011

When n = 9, the test doesn't return 1 for any a <= n - 1 Which makes me wonder the validity of the following 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.

So I dig deeper and find this answer on StackOverflow: https://stackoverflow.com/a/59834347/4632735 Here is my attempt trying to convert the answer to scheme code

``` (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)))))
```