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