265 lines
8.7 KiB
Scheme
265 lines
8.7 KiB
Scheme
(display "ex-1.21") (newline)
|
|
|
|
(define (smallest-divisor n)
|
|
(find-divisor n 2))
|
|
|
|
(define (find-divisor n test-divisor)
|
|
(cond ((> (square test-divisor) n) n)
|
|
((divides? test-divisor n) test-divisor)
|
|
(else (find-divisor n (+ test-divisor 1)))))
|
|
|
|
(define (divides? a b) (= (remainder b a) 0))
|
|
|
|
(define (prime? n) (= n (smallest-divisor n)))
|
|
|
|
(display "(smallest-divisor 33) = ")
|
|
(display (smallest-divisor 33))
|
|
(newline)
|
|
(display "(prime? 37) = ")
|
|
(display (prime? 37))
|
|
(newline)
|
|
|
|
;a^n = a (mod n) ; fermat's theorem
|
|
(define (expmod base exp m)
|
|
(cond ((= exp 0) 1)
|
|
((even? exp)
|
|
(remainder (square (expmod base (/ exp 2) m))
|
|
m))
|
|
(else
|
|
(remainder (* base (expmod base (- exp 1) m))
|
|
m))))
|
|
|
|
(display "(expmod 11 23 3) = ")
|
|
(display (expmod 11 23 3))
|
|
(newline)
|
|
|
|
(define (fermat-test n)
|
|
(define (try-it a)
|
|
(= (expmod a n n) a))
|
|
(try-it (+ 1 (random (- n 1)))))
|
|
|
|
(define (fast-prime? n times)
|
|
(cond ((= times 0) #t)
|
|
((fermat-test n) (fast-prime? n (- times 1)))
|
|
(else #f)))
|
|
|
|
(display "(fast-prime? 11) = ")
|
|
(display (fast-prime? 11 1))
|
|
|
|
(display "\n\nex-1.22\n")
|
|
(display "(smallest-divisor 199) = ")
|
|
(display (smallest-divisor 199)) (newline)
|
|
(display "(smallest-divisor 1999) = ")
|
|
(display (smallest-divisor 1999)) (newline)
|
|
(display "(smallest-divisor 19999) = ")
|
|
(display (smallest-divisor 19999)) (newline)
|
|
|
|
(newline) (display "ex-1.23 - more tests in comments")
|
|
|
|
(define (timed-prime-test n)
|
|
(newline)
|
|
(display n)
|
|
(start-prime-test n (runtime)))
|
|
(define (start-prime-test n start-time)
|
|
(if (prime? n)
|
|
(report-prime (- (runtime) start-time))))
|
|
(define (report-prime elapsed-time)
|
|
(display " *** ")
|
|
(display elapsed-time)
|
|
)
|
|
|
|
(define (search-for-primes n a)
|
|
(cond ((= a 3) 0)
|
|
((prime? n) (timed-prime-test n) (search-for-primes (+ n 1) (+ a 1)))
|
|
((search-for-primes (+ n 1) a))))
|
|
|
|
(newline) (display "search-for-primes 1000")
|
|
(search-for-primes 1000 0)
|
|
|
|
;(newline)(newline) (display "search-for-primes 100000000")
|
|
;(search-for-primes 100000000 0)
|
|
;
|
|
;(newline)(newline) (display "search-for-primes 1000000000")
|
|
;(search-for-primes 1000000000 0)
|
|
;
|
|
;(newline)(newline) (display "search-for-primes 10000000000")
|
|
;(search-for-primes 10000000000 0)
|
|
;
|
|
;(newline)(newline) (display "search-for-primes 100000000000")
|
|
;(search-for-primes 100000000000 0)
|
|
|
|
; It is amazing how exact the time represents the expected increase. If we
|
|
; increase n by factor 10 the time increase by factor sqrt(10) = 3.3. I knew
|
|
; this, but it is cool to see it.
|
|
; search-for-primes 1000000000
|
|
; 1000000007 *** 4.0000000000000036e-2
|
|
; 1000000009 *** .04999999999999999
|
|
; 1000000021 *** 5.0000000000000044e-2
|
|
;
|
|
; search-for-primes 10000000000
|
|
; 10000000019 *** .15000000000000002
|
|
; 10000000033 *** .1399999999999999
|
|
; 10000000061 *** .1399999999999999
|
|
;
|
|
; search-for-primes 100000000000
|
|
; 100000000003 *** .4500000000000002
|
|
; 100000000019 *** .43999999999999995
|
|
; 100000000057 *** .4299999999999997
|
|
|
|
(display "ex-1.23") (newline)
|
|
|
|
(define (next n)
|
|
(if (= n 2) 3 (+ n 2)))
|
|
|
|
(define (find-divisor n test-divisor)
|
|
(cond ((> (square test-divisor) n) n)
|
|
((divides? test-divisor n) test-divisor)
|
|
(else (find-divisor n (next test-divisor)))))
|
|
|
|
;(timed-prime-test 100000000003)
|
|
;(timed-prime-test 100000000019)
|
|
;(timed-prime-test 100000000057)
|
|
|
|
; I did not really expect the time to half, but to 2/3. That is because we skip
|
|
; one number per every three numbers. For example, instead of 3,4,5 we do 3,5
|
|
; 100000000003 *** .30000000000000004
|
|
; 100000000019 *** .29
|
|
; 100000000057 *** .28
|
|
|
|
(newline) (display "ex-1.24")
|
|
|
|
(define (prime? n) (fast-prime? n 100))
|
|
|
|
(define (start-prime-test n start-time)
|
|
(if (prime? n)
|
|
(report-prime (- (runtime) start-time))))
|
|
|
|
(newline)(newline) (display "search-for-primes 1000000000000")
|
|
(search-for-primes 1000000000000 0)
|
|
|
|
(newline)(newline) (display "search-for-primes 1000000000000000000000")
|
|
(search-for-primes 1000000000000000000000000 0)
|
|
|
|
; This very nicely shows that in order to double the execution time we have to
|
|
; square the input. For example:
|
|
|
|
;search-for-primes 1000000000000
|
|
; 1.9999999999999962e-2 / 9.999999999999981e-3 = 2
|
|
|
|
; I feel like all discrepancies must be explained by the random function. If
|
|
; the factor is bigger it takes more steps.
|
|
;1000000000039 *** 9.999999999999981e-3
|
|
;1000000000061 *** 9.999999999999981e-3
|
|
;1000000000063 *** 9.999999999999981e-3
|
|
|
|
;search-for-primes 1000000000000000000000
|
|
;1000000000000000000000007 *** 1.9999999999999962e-2
|
|
;1000000000000000000000049 *** .02999999999999997
|
|
;1000000000000000000000121 *** 3.0000000000000027e-2
|
|
|
|
(newline)(newline) (display "ex-1.25") (newline)
|
|
|
|
;(define (square m)
|
|
; (display "square ")(display m)(newline)
|
|
; (* m m))
|
|
|
|
(define (expmod-naiv base exp m)
|
|
(remainder (fast-expt base exp) m))
|
|
|
|
; Terminates quickly.
|
|
(display "(expmod 10000000000 1000000 100) ; runs quickly") (newline)
|
|
(expmod 10000000000 1000000 100)
|
|
|
|
; Scheme can handle arbitrary precision arithmetic. However,
|
|
; when the numbers get really big that takes more time. When
|
|
; using expmod the number to be squared stays lower than the
|
|
; prime to be tested.
|
|
; When we do the naiv approach the numbers become really big
|
|
; and the same call will not terminate.
|
|
(display "(expmod-naiv 10000000000 1000000 100) ; doesn't terminate")
|
|
;(expmod-naiv 10000000000 1000000 100)
|
|
|
|
(newline)(newline) (display "ex-1.26 - see comment") (newline)
|
|
; With that function for every of the log(n) steps the
|
|
; function is called twice which means the complexity
|
|
; grows exponantially which gives us log(n)^2 a O(n) process.
|
|
|
|
(newline) (display "ex-1.27") (newline)
|
|
|
|
; Do fermat tests for all a in [2, p - 1].
|
|
(define (fermat-test-all n) (fermat-test-all-iter 2 n))
|
|
|
|
(define (fermat-test-all-iter x n)
|
|
(cond ((= x n) #t)
|
|
((= (expmod x n n) x) (fermat-test-all-iter (+ x 1) n))
|
|
(else (display x)(newline) #f)))
|
|
|
|
|
|
; 561, 1105, 1729, 2465, 2821, and 6601
|
|
(display (fermat-test-all 561)) (newline)
|
|
(display (fermat-test-all 1105)) (newline)
|
|
(display (fermat-test-all 1729)) (newline)
|
|
(display (fermat-test-all 2465)) (newline)
|
|
(display (fermat-test-all 2821)) (newline)
|
|
(display (fermat-test-all 6601)) (newline)
|
|
; I am really unhappy with that exercise, because it wouldn't
|
|
; worked if we used the condition a ** (p - 1) = 1 (mod p) and
|
|
; I don't know the exact reason for that.
|
|
|
|
(newline) (display "ex-1.28") (newline)
|
|
|
|
; 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 an-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.
|
|
|
|
(define (is-nontrivial-square n m)
|
|
(cond ((= n 1) n)
|
|
((= n (- m 1)) n)
|
|
((= (remainder (square n) m) 1) 0) ; non-trivial square root
|
|
(else n)))
|
|
|
|
(define (expmod-mr base exp m)
|
|
(cond ((= exp 0) 1)
|
|
((even? exp)
|
|
(remainder (square (is-nontrivial-square (expmod-mr base (/ exp 2) m) m))
|
|
m))
|
|
(else
|
|
(remainder (* base (expmod-mr base (- exp 1) m))
|
|
m))))
|
|
|
|
(define (miller-rabin-test n)
|
|
(define (try-it a)
|
|
(= (expmod-mr a (- n 1) n) 1))
|
|
(try-it (+ 1 (random (- n 1)))))
|
|
|
|
(define (fast-prime? n times)
|
|
(cond ((= times 0) #t)
|
|
((miller-rabin-test n) (fast-prime? n (- times 1)))
|
|
(else #f)))
|
|
|
|
(display (fast-prime? 17 10)) (newline)
|
|
(display (fast-prime? 561 10)) (newline)
|
|
(display (fast-prime? 1105 10)) (newline)
|
|
(display (fast-prime? 1729 10)) (newline)
|
|
(display (fast-prime? 231082301002031230 10)) (newline)
|
|
(display (fast-prime? 1000000000000000000000121 10)) (newline)
|
|
(display (fast-prime? 2 10)) (newline)
|
|
(display (fast-prime? 4 10)) (newline)
|
|
|