Dan's Thoughts Thinking somewhat carefully

Exercise 1.35 of SICP

Exercise 1.35: Show that the golden ratio (section 1.2.2) is a fixed point of the transformation , and use this fact to compute by means of the fixed-point procedure.

(define tolerance 0.00001)
(define (fixed-point f first-guess)
  (define (close-enough? v1 v2)
    (< (abs (- v1 v2)) tolerance))
  (define (try guess)
    (let ((next (f guess)))
      (if (close-enough? guess next)
          next
          (try next))))
  (try first-guess))

> (fixed-point (lambda (x) (+ 1 (/ 1 x))) 1.0)
1.6180327868852458
Actual value of golden ratio: 1.61803399
Tolerance: 0.0001

Clearly is in the defined error tolerance.

Exercise 1.31 of SICP

Exercise 1.31:
a. The sum procedure is only the simplest of a vast number of similar abstractions that can be captured as higher-order procedures. Write an analogous procedure called product that returns the product of the values of a function at points over a given range. Show how to define factorial in terms of product. Also use product to compute approximations to using the formula

b. If your product procedure generates a recursive process, write one that generates an iterative process. If it generates an iterative process, write one that generates a recursive process.

Part A:

Notice that the product function looks almost the same as sum. The only things truly different are the multiplication operation and the base case being multiplicative identity. Factorial doesn't require much except the identity function and a trivial next function.

It's easier to think about Wallis Product product if it's written slightly differently:

(define (product term a next b)
  (if (> a b)
    1
    (* (term a) (product term (next a) next b))))
 
(define (factorial n)
  (define (identity x) x)
  (define (next x) (+ 1 x))
  (product identity 1 next n))
 
(define (wallis n)
  (define (next n) (+ n 1))
  (define (term-num i) (if (even? i) (+ i 2) (+ i 1)))
  (define (term-den i) (if (even? i) (+ i 1) (+ i 2)))
  (define (fraction i) (/ (term-num i) (term-den i)))
  (* 4.0 (product fraction 1 next n)))

> (wallis 5)
2.9257142857142857
> (wallis 50)
3.1719440917427058
> (wallis 500)
3.1447232866889245
> (wallis 5000)
3.1419067029355268
> (wallis 10000)
3.1417497057380523

While the results do tend to they do so very slowly.

Part B:

This is the iterative version of product. It looks similar to the iterative version of sum.

(define (wallis n)
  (define (next n) (+ n 1))
  (define (term-num i) (if (even? i) (+ i 2) (+ i 1)))
  (define (term-den i) (if (even? i) (+ i 1) (+ i 2)))
  (define (fraction i) (/ (term-num i) (term-den i)))
  (* 4.0 (product-iter fraction 1 next n)))
 
(define (product-iter term a next b)
  (define (iter a result)
    (if (> a b)
      result
      (iter (next a) (* (term a) result))))
  (iter a 1))

Exercise 1.29 of SICP

Exercise 1.29: Simpson's Rule is a more accurate method of numerical integration than the method illustrated above. Using Simpson's Rule, the integral of a function f between a and b is approximated as

where , for some even integer n, and . (Increasing n increases the accuracy of the approximation.) Define a procedure that takes as arguments f, a, b, and n and returns the value of the integral, computed using Simpson's Rule. Use your procedure to integrate cube between 0 and 1 (with n = 100 and n = 1000), and compare the results to those of the integral procedure shown above.

(define (cube x) (* x x x))
(define (simpson-integral f a b n)
  (define h (/ (- b a) n))
  (define (next k) (+ k 1))
  (define (coeff k)
    (cond ((or (= k 0) (= k n)) 1)
          ((even? k) 2)
          (else 4)))
  (define (term k) 
    (* (coeff k) (f (+ a (* k h)))))
  (* (/ h 3.0)
     (sum term 0 next n)))
 
(define (sum term a next b)
  (if (> a b)
    0
    (+ (term a)
       (sum term (next a) next b))))

There are several things I want to note here. All the auxiliary functions to represent the leading coefficients as well as h make life a lot easier. It's also important to note that the variables a and b in sum are not performing the same role as in the integral example earlier. In this case they are just counting indices from a = 0 to b = n. It might not be immediately obvious: When the function term is called from sum a (the start of the integrating interval) is bound to it from the original call to simpson-integral . Term only receives the index which it promptly evaluates.

Exercise 1.28 of SICP

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 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 (miller-rabin-test n)
  (define (random n) (random-integer n))
  (define (expmod base exp m)
    (define (square-mod x) (remainder (* x x) m))
    (define (square-signal-root x)
      (if (and 
            (not (or (= 1 x) (= x (- m 1))))
            (= 1 (square-mod x)))
        0
        (square-mod x)))
    (cond ((= exp 0) 1)
          ((even? exp)
           (square-signal-root (expmod base (/ exp 2) m)))
          (else
            (remainder (* base (expmod base (- exp 1) m))
                       m))))  
  (define (try-it a)
    (= (expmod 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)))

Carmichael Numbers: 561, 1105, 1729, 2465, 2821, 6601
Primes: 2, 3, 10169, 31337, 1000249, 382176531747930913347461352433
Non-primes: 6, 27, 49, 1024

> (fast-prime? 561 10)
#f
> (fast-prime? 1105 10)
#f
> (fast-prime? 1729 10)
#f
> (fast-prime? 2465 10)
#f
> (fast-prime? 2821 10)
#f
> (fast-prime? 6601 10)
#f
> (fast-prime? 2 10)
#t
> (fast-prime? 3 10)
#t
> (fast-prime? 10169 10)
#t
> (fast-prime? 31337 10)
#t
> (fast-prime? 1000249 10)
#t
> (fast-prime? 382176531747930913347461352433 10)
#t
> (fast-prime? 6 10)
#f
> (fast-prime? 27 10)
#f
> (fast-prime? 49 10)
#f
> (fast-prime? 1024 10)
#f

Exercise 1.22 of SICP

Exercise 1.22: Most Lisp implementations include a primitive called runtime (my implementation of scheme doesn't. I define a work around) that returns an integer that specifies the amount of time the system has been running (measured, for example, in microseconds). The following timed-prime-test procedure, when called with an integer n, prints n and checks to see if n is prime. If n is prime, the procedure prints three asterisks followed by the amount of time used in performing the test.

(define (runtime) (time->seconds (current-time)))
(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))

Using this procedure, write a procedure search-for-primes that checks the primality of consecutive odd integers in a specified range. Use your procedure to find the three smallest primes larger than 1000; larger than 10,000; larger than 100,000; larger than 1,000,000. Note the time needed to test each prime. Since the testing algorithm has order of growth of , you should expect that testing for primes around 10,000 should take about times as long as testing for primes around 1000. Do your timing data bear this out? How well do the data for 100,000 and 1,000,000 support the prediction? Is your result compatible with the notion that programs on your machine run in time proportional to the number of steps required for the computation?

> (timed-prime-test 68720001023)
68720001023 *** .4850001335144043

> (timed-prime-test 70368760954879)
70368760954879 *** 16.844000101089478

Since this exercise is a bit dated testing numbers on the order 10⁶ for primality doesn't give my system any noticeable workout. I test my ranges starting from 10¹⁰.
> (search-for-primes (expt 10 10) (- (expt 10 11) 1))

10000000019 *** .14100003242492676
10000000033 *** .1399998664855957
10000000061 *** .14000010490417483

According to the theoretical calculations the numbers on the order of 10¹¹ should take seconds.
Error: 17%

> (search-for-primes (expt 10 11) (- (expt 10 12) 1))

100000000003 *** .5310001373291016
100000000019 *** .5320000648498535
100000000057 *** .54699993133544923

seconds.

> (search-for-primes (expt 10 12) (- (expt 10 13) 1))

1000000000039 *** 1.7660000324249268
1000000000061 *** 1.75
1000000000063 *** 1.76500010490417483

seconds.
Error: 2%

> (search-for-primes (expt 10 13) (- (expt 10 14) 1))

10000000000037 *** 5.672000169754028
10000000000051 *** 5.672000169754028
10000000000099 *** 5.65599989891052253

It seems that is indeed a good predictor of the time it takes to check for primality of the algorithm.

(define (search-for-primes beg end)
  (define (search beg end num-primes)
    (cond ((= num-primes 3) num-primes)
          ((> beg end) num-primes)
          ((even? beg) (search (+ beg 1) end num-primes))
          ((prime? beg) (timed-prime-test beg) 
                        (search (+ beg 2) end (+ 1 num-primes)))
          (else
            (search (+ beg 2) end num-primes))))
  (search beg end 0))
 
(define (runtime) (time->seconds (current-time)))
(define (square x) (* x x))
(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)))
 
(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))

Exercise 1.13 of SICP

Exercise 1.13: Prove that Fib(n) is the closest integer to n/5, where . Hint: Let Use induction and the definition of the Fibonacci numbers (see section 1.2.2) to prove that .

F_n is a linear second-order recurrence with constant coefficients. The characteristic equation for it is:

This means that the closed form solution is of the form
The constants and are determined by the initial conditions bellow.






The closed expression has the following form:

Now that the derivation of the closed form for F_n is done here's the proof by induction.





Perform the induction step:
Assume:


Because

The proof that is straightforward.

for example and

Clearly for n>30 is a good approximation.