1 665c255d 2023-08-04 jrmu ;; calculate base^exp modulo m
2 665c255d 2023-08-04 jrmu (define (expmod base exp m)
3 665c255d 2023-08-04 jrmu (cond ((= exp 0) 1)
5 665c255d 2023-08-04 jrmu (remainder (square (expmod base (/ exp 2) m)) m))
6 665c255d 2023-08-04 jrmu (else (remainder (* base (expmod base (- exp 1) m)) m))))
8 665c255d 2023-08-04 jrmu ;; tests if integer n passes fermat's little theorem
9 665c255d 2023-08-04 jrmu (define (fermat-prime? n)
10 665c255d 2023-08-04 jrmu (define (fermat-test a)
11 665c255d 2023-08-04 jrmu (cond ((= a n) #t)
12 665c255d 2023-08-04 jrmu ((not (= (expmod a n n) a)) #f)
13 665c255d 2023-08-04 jrmu (else (fermat-test (+ a 1)))))
14 665c255d 2023-08-04 jrmu (fermat-test 1))
16 665c255d 2023-08-04 jrmu (define (list-primes upper)
17 665c255d 2023-08-04 jrmu (define (test i)
18 665c255d 2023-08-04 jrmu (cond ((= i upper) (display "Finished"))
19 665c255d 2023-08-04 jrmu ((fermat-prime? i) (begin (display i)
20 665c255d 2023-08-04 jrmu (newline))))
21 665c255d 2023-08-04 jrmu (test (+ i 1)))
24 665c255d 2023-08-04 jrmu (define (test-case actual expected)
25 665c255d 2023-08-04 jrmu (load-option 'format)
27 665c255d 2023-08-04 jrmu (format #t "Actual: ~A Expected: ~A" actual expected))
29 665c255d 2023-08-04 jrmu ;; (test-case (fermat-prime? 2) #t)
30 665c255d 2023-08-04 jrmu ;; (test-case (fermat-prime? 3) #t)
31 665c255d 2023-08-04 jrmu ;; (test-case (fermat-prime? 4) #f)
32 665c255d 2023-08-04 jrmu ;; (test-case (fermat-prime? 5) #t)
33 665c255d 2023-08-04 jrmu ;; (test-case (fermat-prime? 6) #f)
34 665c255d 2023-08-04 jrmu ;; (test-case (fermat-prime? 7) #t)
35 665c255d 2023-08-04 jrmu ;; (test-case (fermat-prime? 8) #f)
36 665c255d 2023-08-04 jrmu ;; (test-case (fermat-prime? 9) #f)
39 665c255d 2023-08-04 jrmu ;; (list-primes 10000)
41 665c255d 2023-08-04 jrmu ;; Carmichael Numbers
42 665c255d 2023-08-04 jrmu (test-case (fermat-prime? 561) #f)
43 665c255d 2023-08-04 jrmu (test-case (fermat-prime? 1105) #f)
44 665c255d 2023-08-04 jrmu (test-case (fermat-prime? 1729) #f)
45 665c255d 2023-08-04 jrmu (test-case (fermat-prime? 2465) #f)
46 665c255d 2023-08-04 jrmu (test-case (fermat-prime? 2821) #f)
47 665c255d 2023-08-04 jrmu (test-case (fermat-prime? 6601) #f)
50 665c255d 2023-08-04 jrmu ;; Exercise 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.
52 665c255d 2023-08-04 jrmu (expmod base exp m)
53 665c255d 2023-08-04 jrmu (= (expmod a (n-1) n) 1)
55 665c255d 2023-08-04 jrmu ;; calculate base^exp modulo m but signal if
56 665c255d 2023-08-04 jrmu (define (expmod base exp m)
57 665c255d 2023-08-04 jrmu (cond ((= exp 0) 1)
58 665c255d 2023-08-04 jrmu ((even? exp)
59 665c255d 2023-08-04 jrmu (remainder (square (expmod base (/ exp 2) m)) m))
60 665c255d 2023-08-04 jrmu (else (remainder (* base (expmod base (- exp 1) m)) m))))
