<< Previous exercise (1.27) | sicp-solutions | Next exercise (1.29) >>
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 ;;;; <https://creativecommons.org/licenses/by-sa/4.0/>. (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))
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)))))
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.
ANSWER