SICP through §1.2.6 + REPL development

All code available on GitHub. Drawings/diagrams available here.

One thing that has made going through this book super fun so far is REPL-driven development. My current workflow uses VSCode’s Calva plugin and it is truly powerful. For instance, I can evaluate parts of an expression easily and also instrument code and do step-through debugging. I also took some time to learn about Paredit, which has made working with all those S-expressions much easier! Anyway, on to the meat of the article:

One key takeaway from §1.2 is that we must not confuse Recursive Procedures with Recursive Processes. The following exercises are an application of this concept:

;; Exercise 1.11.  A function f is defined by the rule that f(n) = n if n<3 and f(n) = f(n - 1) + 2f(n - 2) + 3f(n - 3) if n> 3. Write a procedure that computes f by means of a recursive process. Write a procedure that computes f by means of an iterative process.
(defn ex111-r [n]
  (cond (< n 3) n
        (>= n 3) (+ (ex111-r (- n 1))
                    (* 2 (ex111-r (- n 2)))
                    (* 3 (ex111-r (- n 3))))
        :else (throw (new AssertionError "invalid n"))))

(ex111-r 10)

(defn ex111-i [n]
  (let [iter (fn [a b c count]
               (if (zero? count) c
                   (recur b c (+ (* 3 a) (* 2 b) c) (dec count))))]
    (iter 0 1 2 (- n 2))))

(ex111-i 10)

This one uses a trick I learnt from doing 4Clojure exercises - partition 2 1 gives overlapping pairs in a list which work well for Pascal’s triangle.

;; Exercise 1.12. Write a procedure that computes elements of Pascal's triangle by means of a recursive process.
(defn pascal [n]
  (let [pp (fn [i prev]
             (cond (< n 0) (throw (new AssertionError "invalid n"))
                   (= n 1) [1]
                   (= n 2) [1 1]
                   (= i n) prev
                   :else (let [pairs (partition 2 1 prev)
                               new-inner (map #(apply + %) pairs)]
                           (recur (inc i) (concat [1] new-inner [1])))))]
    (pp 1 [])))

(assert (= (pascal 5) [1 4 6 4 1]))

Ah yes; invariants! Designing a loop with the right invariants can often assure you of the correctness of the loop. Exercises 1.16 and 1.18 are a great example of this idea.

;; Exercise 1.16  Design a procedure that evolves an iterative exponentiation process that uses successive squaring and uses a logarithmic number of steps, as does fast-expt.
(defn fast-expt [b n]
  (cond (zero? n) 1
        (even? n) (* (fast-expt b (/ n 2)) (fast-expt b (/ n 2)))
        :else (* b (fast-expt b (- n 1)))))
(fast-expt 2 5)

(defn fast-expt-iter [b n]
  ;; a﹒b^n is invariant
  (let [f (fn [b n a]
            (cond
              (= n 0) a
              (even? n) (recur (* b b) (/ n 2) a)
              :else (recur b (dec n) (* b a))))]
    (f b n 1)))
(fast-expt-iter 2 5)
;; Exercise 1.18.  Using the results of exercises 1.16 and 1.17, devise a procedure that generates an iterative process for multiplying two integers in terms of adding, doubling, and halving and uses a logarithmic number of steps.40

;; After thinking of the invariant, it comes naturally.
(defn fast-mult-iter [a b n]
  ;; a﹒b + n is invariant
  (let [double (fn [a] (* a 2))
        halve (fn [a] (/ a 2))]
    (cond
      (zero? b) n
      (even? b) (recur (double a) (halve b) n)
      :else (recur a (dec b) (+ n a)))))
(fast-mult-iter 11 3 0)

A logarithmic algorithm for computing Fibonacci numbers, super cool! My derivation of the transform is here: https://notability.com/n/~S4C78hC359sd2EDdHWkW.

;; Exercise 1.19
(defn fib-iter [a b p q count]
  (cond (= count 0) b
        (even? count) (recur
                       a
                       b
                       (+ (* q q) (* p p))
                       (+ (* q q) (* 2 p q))
                       (/ count 2))
        :else (recur (+ (* b q) (* a q) (* a p))
                     (+ (* b p) (* a q))
                     p
                     q
                     (dec count))))
(fib-iter 1 0 0 1 10)

For the following exercises, we use Clojure’s built-in time function, which prints out the runtime of the expression as a side-effect.

;; Exercise 1.21.  Use the smallest-divisor procedure to find the smallest divisor of each of the following numbers: 199, 1999, 19999.
(defn smallest-divisor-slow [n]
  (first (filter #(zero? (mod n %))
                 (range 2 (inc n)))))

(defn smallest-divisor [n]
  (first (filter #(zero? (mod n %))
                 (concat '(2) (range 3 (inc n) 2)))))

(smallest-divisor 199) ;; 199
(smallest-divisor 1999) ;; 1999
(smallest-divisor 19999) ;; 7

;; Ex 1.23
(defn prime-slow? [n] (= n (smallest-divisor-slow n)))
(defn prime? [n] (= n (smallest-divisor n)))
(time (dotimes [n 10] (prime? 1000037))) ;; 555 ms
(time (dotimes [n 10] (prime-slow? 1000037))) ;; 1073 ms
;; 2x speedup, as expected!

Here we apply another trick to measure the runtime of the code - we use doall since Clojure’s sequences are lazy and we want to force evaluation.

;; Ex 1.21
(defn fermat-test [n]
  (let [expmod (fn expmod [base exp m]
                 (cond (zero? exp) 1
                       (even? exp) (let [x (expmod base (/ exp 2) m)] (mod (* x x) m))
                       :else (mod (* base (expmod base (dec exp) m)) m)))
        try-it  (fn [a] (= (expmod a n n) a))]
    (try-it (+ 1 (rand-int (dec n))))))
;; Uses the mod property 
(fermat-test 19999)

(defn fast-prime? [n times]
  (cond (zero? times) true
        (fermat-test n) (fast-prime? n (- times 1))
        :else false))

;; Ex 1.22
(defn search-for-primes [lower-bound num-primes-wanted]
  (take num-primes-wanted
        (filter #(fast-prime? % 8)
                (filter odd? (range lower-bound (java.lang.Integer/MAX_VALUE))))))

(defn search-for-primes-slow [lower-bound num-primes-wanted]
  (take num-primes-wanted
        (filter prime?
                (filter odd? (range lower-bound (java.lang.Integer/MAX_VALUE))))))

;; This is approximately O(sqrt(n)) growth.
(time (doall (search-for-primes-slow 1000 3))) ;;/ (1009 1013 1019) 0.6 msecs"
(time (doall (search-for-primes-slow 10000 3))) ;; (10007 10009 10037) 4 msecs
(time (doall (search-for-primes-slow 100000 3))) ;; (100003 100019 100043) 30 msecs"
(time (doall (search-for-primes-slow 1000000 3))) ;; (1000003 1000033 1000037) 224 ms

;; Ex 1.24
;; This is approximately O(log(n)) growth and demonstrates the power of logarithmic runtimes - a multiplication of the input's cardinality results in an additive runtime increase.
(trace/trace (time (doall (search-for-primes 1000 3)))) ;; (1009 1013 1019) 0.24 msecs"
(time (doall (search-for-primes 10000 3))) ;; (10007 10009 10037) 0.50 msecs
(time (doall (search-for-primes 100000 3))) ;; (100003 100019 100043) 0.57 msecs"
(time (doall (search-for-primes 1000000 3))) ;; (1000003 1000033 1000037) 0.6 msecs"
;; Exercise 1.25.  Alyssa P. Hacker complains that we went to a lot of extra work in writing expmod. After all, she says, since we already know how to compute exponentials, we could have simply written: 
;; (define (expmod base exp m)
;;     (remainder (fast-expt base exp) m))

;; Theoretically (in terms of math), it would work. Computationally it fails since the numbers (without the intermediate mods) are too large and either overflow or take too long.
;; Exercise 1.26
;;  Previously, the process was linear recursive and you halved the search space on each level (the /2). Now, every time you halve the seach space, you make 2 recursive calls at the same level. So the runtime is now O(N) instead of O(log N).s
;; Exercise 1.28
;; Proof that if there is nontrivial square root of n it is not prime
;; https://crypto.stanford.edu/pbc/notes/numbertheory/poly.html
(defn miller-rabin [n]
  (let [expmod (fn expmod [base exp m]
                 (cond (zero? exp) 1
                       (even? exp) (let [x (expmod base (/ exp 2) m)
                                         sq (mod (* x x) m)]
                                     (if (and (not= x 1)
                                              (not= x (dec m))
                                              (= sq 1)) 0
                                         sq))
                       :else (mod (* base (expmod base (dec exp) m)) m)))
        try-it  (fn [a] (= (expmod a (dec n) n) 1))]
    (try-it 3)
    ;; (every? true? (map try-it (range 2 n)))
    ))