Chao's Random Thoughts

Keep looking, don't settle

Scala Cake Pattern

The summary of my thoughts on the Scala Cake Pattern.

SICP Section 1.3

Procedures as Arguments

Why 'high order procedures' is useful?

Often the same programming pattern will be used with a number of different procedures. To express such patterns as concepts, we will need to construct procedures that can accept procedures as arguments or return procedures as values. Procedures that manipulate procedures are called high-order procedures, they can serve as powerful abstraction mechanisms, vastly increaseing the expressive power of our language.

Lambda

  • Using lambda to create anonymous functions:
    • (define (plus4 x) (+ x 4)) is equivalent to (define plus4 (lambda (x) (+ x 4)))
  • Using lambda(let) to create local variables
Sometimes we can use internal definitions to get the same effect as with let.

And in which cases we couldn't get the same effect?

Abstractions and first-class procedures

As programmers, we should be alert to oppotunities to identify the underlying abstractions in our programs and to build upon then and generalize them to create more powerful abstractions. This is not to say one should always write programs in the most abstract way possible; expert programmers know how to choose the level of abstraction appropriate to their task. But it's important to be able to think in terms of these abstractions, so that we can be ready to apply them in new contexts.

Exercises

1.29 Answer:

This is not very hard, just simply translate the simpson rule definition to code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

(define (simpson-rule f a b n)
  ; still use 'define' since we haven't learned 'let' till now
  (define h (/ (- b a) n))
  ; the procedure to compute 'yk', the helper method used by 'term'
  (define (y k) (f (+ a (* k h))))
  ; get the 'next and 'term for the sum procedure
  (define (next k) (+ k 1))
  (define (term k)
    (* (cond ((or (= k 1) (= k n)) 1)
             ((even? k) 2)
             ((odd? k) 4))
       (y k)))
  (/ (* h (sum term 0 next n)) 3))

(define (sum term a next b)
  (if (> a b)
    0
    (+ (term a)
       (sum term (next a) next b))))

(define (integral f a b dx)
  (define (add-dx x) (+ x dx))
  (* (sum f (+ a (/ dx 2.0)) add-dx b) dx))

(define (cube x) (* x x x))

(simpson-rule cube 0 1 100.0)
(simpson-rule cube 0 1 1000.0)
(integral cube 0 1 0.01)
(integral cube 0 1 0.001)

;; the output:
; 0.24999998999999987
; 0.24999999999900027
; 0.24998750000000042
; 0.249999875000001

1.30 Answer:

Put both the recursive and iterative versions here to compare:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

(define (sum term a next b)
  (if (> a b)
    0
    (+ (term a)
       (sum term (next a) next b))))

(define (sum-iter term a next b)
  (define (iter a result)
    (if (> a b)
        result
        (iter (next a) (+ (term a) result))))
  (iter a 0))


(sum (lambda (x) x) 1 (lambda (x) (+ x 1)) 10)
(sum-iter (lambda (x) x) 1 (lambda (x) (+ x 1)) 10)

(sum (lambda (x) (* x x)) 1 (lambda (x) (+ x 1)) 5)
(sum-iter (lambda (x) (* x x)) 1 (lambda (x) (+ x 1)) 5)

They produce the same results as expected.

1.31 Answer:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

(define (product term a next b)
  (if (> a b)
    1
    (* (term a)
       (product term (next a) next b))))

(define (product-iter term a next b)
  (define (iter cur result)
    (if (> cur b)
      result
      (iter (next cur) (* (term cur) result))))
  (iter a 1))

(define (factorial n)
  (define (identity x) x)
  (product identity 1 inc n))

(define (inc x) (+ x 1))


(product (lambda (x) x) 1 (lambda (x) (+ x 1)) 6)
(product-iter (lambda (x) x) 1 (lambda (x) (+ x 1)) 6)
(factorial 6)
(factorial 10)
(product (lambda (x) (* x x)) 1 (lambda (x) (+ x 1)) 5)
(product-iter (lambda (x) (* x x)) 1 (lambda (x) (+ x 1)) 5)

;; compute the approximations to PI using the formula:
;; PI   2*4*4*6*6*8...
;; -- = --------------
;; 4    3*3*5*5*7*7...

;; here, term is (2n/(2n + 1) * 2(n+1) /(2n + 1)), n = 1, 2, 3, ...
(/
  (product-iter
    (lambda (n)
      (* (* 2 n) (* 2 (+ n 1.0))))
    1
    inc
    50)
  (product-iter
    (lambda (n)
      (* (+ (* 2 n) 1.0) (+ (* 2 n) 1.0)))
    1
    inc
    50)
)

There might be better ways to translate the formula. The solution above the is not veru accurate and if the upper bound is too large (say 100), it will yield to 'nan+'.

1.32 Answer:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

(define (accumulate combiner null-value term a next b)
  (if (> a b)
    null-value
    (combiner (term a)
              (accumulate combiner null-value term (next a) next b))))

(define (accumulate-iter combiner null-value term a next b)
  (define (iter cur result)
    (if (> cur b)
      result
      (iter (next cur) (combiner (term cur) result))))
  (iter a null-value))

(define (sum term a next b)
  ; (accumulate + 0 term a next b))
  (accumulate-iter + 0 term a next b))

(define (product term a next b)
  ; (accumulate * 1 term a next b))
  (accumulate-iter * 1 term a next b))


(sum (lambda (x) x) 1 (lambda (x) (+ x 1)) 6)
(product (lambda (x) x) 1 (lambda (x) (+ x 1)) 6)
(sum (lambda (x) (* x x)) 1 (lambda (x) (+ x 1)) 5)
(product (lambda (x) (* x x)) 1 (lambda (x) (+ x 1)) 5)

1.33 Answer:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42

(define (filtered-accumulate combiner null-value term a next b filter)
  (if (> a b)
    null-value
    (if (filter a)
      (combiner (term a)
                (filtered-accumulate combiner null-value term (next a) next b filter))
      (filtered-accumulate combiner null-value term (next a) next b filter))))


(define (prime? n)
  (= n (smallest-divisor n)))

(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 (square x) (* x x))

(define (divides? a b)
  (= (remainder b a) 0))

(define (inc x) (+ x 1))

(define (identity x) x)

; return a procedure Integer -> Boolean
; that check if the given integer is relative prime with n
(define (relative-prime? n)
  (define (f a) (= (gcd n a) 1))
  f)

(define (gcd a b)
  (if (= b 0)
    a
    (gcd b (remainder a b))))

(filtered-accumulate + 0 square 2 inc 10 prime?)
(filtered-accumulate * 1 identity 1 inc 10 (relative-prime? 10))

1.34 Answer:

Apply (f f) with get (f 2), while tring to apply (f 2), error will be issued since a procedure is expected but an integer parameter 2 is given.

1.35 Answer:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

(define tolerance 0.00001)

(define (fixed-point f first-guess)
  (define (close-enough x y)
    (< (abs (- x y)) tolerance))
  (define (try guess)
    (let ((next (f guess)))
      (if (close-enough next guess)
        next
        (try next))))
  (try first-guess))


; since the golden-ratio x satisfies that x^2 = x + 1
(define (golden-ratio)
  (fixed-point (lambda (x) (+ 1 (/ 1 x))) 1.0))


(golden-ratio)

1.36 Answer:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70

(define tolerance 0.00001)

(define (fixed-point f first-guess)
  (define (close-enough x y)
    (< (abs (- x y)) tolerance))
  (define (try guess)
    (display "current guess is: " )
    (display guess)
    (newline)
    (let ((next (f guess)))
      (if (close-enough next guess)
        next
        (try next))))
  (try first-guess))


; without average damping
(fixed-point (lambda (x) (/ (log 1000) (log x))) 2.0)

; output below:
; current guess is: 2.0
; current guess is: 9.965784284662087
; current guess is: 3.004472209841214
; current guess is: 6.279195757507157
; current guess is: 3.759850702401539
; current guess is: 5.215843784925895
; current guess is: 4.182207192401397
; current guess is: 4.8277650983445906
; current guess is: 4.387593384662677
; current guess is: 4.671250085763899
; current guess is: 4.481403616895052
; current guess is: 4.6053657460929
; current guess is: 4.5230849678718865
; current guess is: 4.577114682047341
; current guess is: 4.541382480151454
; current guess is: 4.564903245230833
; current guess is: 4.549372679303342
; current guess is: 4.559606491913287
; current guess is: 4.552853875788271
; current guess is: 4.557305529748263
; current guess is: 4.554369064436181
; current guess is: 4.556305311532999
; current guess is: 4.555028263573554
; current guess is: 4.555870396702851
; current guess is: 4.555315001192079
; current guess is: 4.5556812635433275
; current guess is: 4.555439715736846
; current guess is: 4.555599009998291
; current guess is: 4.555493957531389
; current guess is: 4.555563237292884
; current guess is: 4.555517548417651
; current guess is: 4.555547679306398
; current guess is: 4.555527808516254
; current guess is: 4.555540912917957
; 4.555532270803653

; with average damping
(fixed-point (lambda (x) (/ (+ x (/ (log 1000) (log x))) 2)) 2.0)

; output below:
; current guess is: 2.0
; current guess is: 5.9828921423310435
; current guess is: 4.922168721308343
; current guess is: 4.628224318195455
; current guess is: 4.568346513136242
; current guess is: 4.5577305909237005
; current guess is: 4.555909809045131
; current guess is: 4.555599411610624
; current guess is: 4.5555465521473675
; 4.555537551999825

1.37 Answer:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

; k-term finite continued fraction
(define (cont-frac n d k)
  (define (frac i)
    (if (= i k)
      (/ (n i) (d i))
      (/ (n i) (+ (d i) (frac (+ i 1))))))
  (frac 1))

(define (cont-frac-iter n d k)
  (define (iter i result)
    (if (= i 0)
      result
      (iter (- i 1) (/ (n i) (+ (d i) result)))))
  (iter k 0))


; k = 10 produce 0.6179775280898876
(cont-frac (lambda (i) 1.0)
           (lambda (i) 1.0)
           10)

(cont-frac-iter (lambda (i) 1.0)
                (lambda (i) 1.0)
                10)

; k = 11 produce 0.6180555555555556
(cont-frac (lambda (i) 1.0)
           (lambda (i) 1.0)
           11)

(cont-frac-iter (lambda (i) 1.0)
                (lambda (i) 1.0)
                11)


; TODO: implement a procedure to retry different k until the first
; one that produce the value accurate to 4 decimal.

1.38 Answer:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

(define (cont-frac n d k)
  (define (iter i result)
    (if (= i 0)
      result
      (iter (- i 1) (/ (n i) (+ (d i) result)))))
  (iter k 0))


; if    i % 3 == 2, (d i) = 2 * ((i + 1) / 3)
; else  (d i) = 1
(+ (cont-frac (lambda (i) 1.0)
              (lambda (i) (let ((r (remainder i 3)))
                            (if (= r 2)
                              (* 2.0 (/ (+ i 1) 3))
                              1.0)))
              10)
   2)

1.39 Answer:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

; if i = 1, (n i) = -x else (n i) = -x^2
; (d i) = 2 * i - 1
(define (tan-cf x k)
  (cont-frac (lambda (i) (if (= i 1)
                           x
                           (- 0(* x x))))
             (lambda (i) (- (* 2 i) 1))
             k))

(define (cont-frac n d k)
  (define (iter i result)
    (if (= i 0)
      result
      (iter (- i 1) (/ (n i) (+ (d i) result)))))
  (iter k 0))


; tan (pi / 8)
(tan-cf 0.4 10)

; tan (pi / 4)
(tan-cf 0.7853981633974483 10)

1.40 Answer:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50

(define (average-damp f)
  (define (average a b)
    (/ (+ a b) 2))
  (lambda (x) (average x (f x))))

(define (square x) (* x x))

(define tolerance 0.00001)

(define (fixed-point f first-guess)
  (define (close-enough x y)
    (< (abs (- x y)) tolerance))
  (define (try guess)
    (let ((next (f guess)))
      (if (close-enough next guess)
        next
        (try next))))
  (try first-guess))

;; re-write the square-root procedure
(define (sqrt x)
  (fixed-point (average-damp (lambda (y) (/ x y))) 1.0))

;; cube root, similar with above
(define (cube-root x)
  (fixed-point (average-damp (lambda (y) (/ x (square y)))) 1.0))


(define dx 0.00001)

(define (deriv g)
  (lambda (x) (/ (- (g (+ x dx)) (g x)) dx)))

(define (cube x)
  (* x x x))

(define (newton-transform g)
  (lambda (x)
    (- x (/ (g x) ((deriv g) x)))))

(define (newton-method g guess)
  (fixed-point (newton-transform g) guess))

(define (cubic a b c)
  (lambda (x)
    (+ (* x x x) (* a (square x)) (* b x) c)))

; zero of x^3 + x^2 + x + 1
; should be -1
(newton-method (cubic 1 1 1) 1.0)

1.41 Answer:

1
2
3
4
5
6
7
8
9
10
11
12
13
14

;; Note: it's a little confusing at first, 'double' is not apply the function
;; 'f' and double the result, it's 'DOUBLE' the application of the function
;; 'f'.
(define (double f)
  (lambda (x) (f (f x))))

(define (inc x)
  (+ x 1))

((double inc) 1)
(((double double) inc) 1)

; this will get 21
(((double (double double)) inc) 5)

1.42 Answer:

1
2
3
4
5
6
7
8
9

(define (compose f g)
  (lambda (x) (f (g x))))

(define (square x) (* x x))

(define (inc x) (+ x 1))


((compose square inc) 6)

1.43 Answer:

1
2
3
4
5
6
7
8
9
10
11
12
13
14

(define (repeated f n)
  (if (= n 1)
    f
    (compose f (repeated f (- n 1)))))

(define (compose f g)
  (lambda (x) (f (g x))))

(define (square x) (* x x))

(define (inc x) (+ x 1))


((repeated square 2) 5)

1.44 Answer:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

(define (n-fold-smooth f n)
  (repeated (smooth f) n))

(define (smooth f)
  (let ((dx 0.0001))
    (lambda (x)
      (/ (+ (f (- x dx)) (f x) (f (+ x dx))) 3))))

(define (repeated f n)
  (if (= n 1)
    f
    (compose f (repeated f (- n 1)))))

(define (compose f g)
  (lambda (x) (f (g x))))

((n-fold-smooth (lambda (x) (/ 1.0 x)) 3) 6)
((lambda (x) (/ 1.0 x)) 6)

1.45 Answer:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41

(define tolerance 0.00001)

;; from the experiments, the nth root requires at least log2(n) average damps
(define (nth-root x n)
  (fixed-point
    ((repeated average-damp (floor (/ (log n) (log 2))))
    (lambda (y) (/ x (expt y (- n 1)))))
    1.0))

(define (fixed-point f first-guess)
  (define (close-enough? x y)
    (< (abs (- (/ x y) 1)) tolerance))
  (define (try guess)
    (let ((next (f guess)))
      ; (display next)
      ; (newline)
      (if (close-enough? next guess)
        next
        (try next))))
  (try first-guess))

(define (average-damp f)
  (define (average x y)
    (/ (+ x y) 2))
  (lambda (x) (average x (f x))))

(define (repeated f n)
  (if (= n 1)
    f
    (compose f (repeated f (- n 1)))))

(define (compose f g)
  (lambda (x) (f (g x))))


(nth-root 16 4)
(nth-root 32 5)
(nth-root 64 6)
(nth-root 128 7)
(nth-root 256 8)
(nth-root 10000000000000000000000000000000000000000000000000000000000000000 64)

1.46 Answer:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

(define tolerance 0.00001)

(define (iterative-improve enough? improve)
  (define (try guess)
    (let ((next-guess (improve guess)))
      (if (enough? guess next-guess)
        guess
        (try next-guess))))
  try)

(define (fixed-point f first-guess)
  ((iterative-improve close-enough? f) first-guess))

(define (sqrt x)
  ((iterative-improve
      close-enough?
      (lambda (y) (average y (/ x y))))
      1.0))

(define (average x y)
  (/ (+ x y) 2))

(define (close-enough? x y)
  (< (abs (/ (- x y) y)) tolerance))


(sqrt 4)
(sqrt 6)
(sqrt 9)

(fixed-point cos 1.0)

SICP Section 1.2

Recursive Process V.S. Iterative Process

  • Recursive process

    The substitution model of a recursive process reveals a shape of expansion followed by constraction. The expansion occurs as the process builds up a chain of defered operations. The contraction occurs as the operations are actually performed.

  • Iterative process

    An iterative process is one whose state can be summarized by a fixed number of state variables, together with a fixed rule that describes how the state variables should be updated as the process moves from state to state and an (optional) end test that specifies the conditions under which the process should terminate.

In the iterative case, the program variables provide a complete description of the state of the process at any point, while in the recursive case, there is some additional "hidden" information, maintained by the interpreter and not contained in the program variables.

In contrasting iteration and recursion, we must be careful not to confuse the notion of a recursive process with the notion of a recursive procedure. When we describe a procedure as recursive, we are referring to the syntactic fact that the procedure definition refers (either directly or indirectly) to the procedure itself. But when we describe a process as following a pattern that is, say, linearly recursive, we are speaking about how the process evolves, not about the syntax of how a procedure is written. It may seem disturbing that we refer to a recursive procedure as generating an iterative process.

Lame's Theorem

If Euclid's Algorithm requires k steps to compute the GCD of some pair, then the smaller number in the pair must be greater than or equal to the kth Fibonacci number.

Fermat's Little Theorem

If n is a prime number and a is any positive integer less than n, then a raised to the nth power is congruent to a modulo n.

Exercises

1.9 Answer: The first process is recursive and the second one is iterative.

1.10 Answer: from the definition of (A x y), we have:

(A 1 10) = A(0 (A 1 9)) = 2 * (A 1 9) = ... = 2^9 * (A 1 1) = 2^10

(A 2 4) = 65536

(A 3 3) = 65536

(f n) computes 2n

(g n) computes 2^n

(h n) computes 2^2^2... (n times)

1.11 Answer: This is similar with the example of computing Fibonacci number.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

(define (f-recursive n)
  (cond ((= n 0) 0)
        ((= n 1) 1)
        ((= n 2) 2)
        (else (+ (f-recursive (- n 1))
                 (* (f-recursive (- n 2)) 2)
                 (* (f-recursive (- n 3)) 3)))))


(define (f-iterative n)
  (f-iter 0 1 2 n))

(define (f-iter a b c count)
  (cond ((= count 0) a)
        ((= count 1) b)
        ((= count 2) c)
        (else (f-iter b c (+ c (* b 2) (* a 3)) (- count 1)))))

1.12

1
2
3
4
5
6

; both row and col indexes start at 0
(define (pascal-triangle row col)
  (cond ((= col 0) 1)
        ((= row col) 1)
        (else (+ (pascal-triangle (- row 1) (- col 1))
                 (pascal-triangle (- row 1) col)))))

1.15 Answer:

a. from the the procedure definition, suppose p will be called t times, we have:

a / (3^t) <= 0.1, which leads to:

log(10a) <= t, so t = ceiling (log(10a)) = 5, (the base of the log is 3)

b. both the space and number of steps depend on how many times p is called, so it's O(t) = O(log(a)).

1.16

1
2
3
4
5
6
7
8
9
10
11
12

(define (fast-expt b n)
  (fast-expt-iter b n 1))

(define (fast-expt-iter b n product)
  (cond ((= n 0) product)
        ((even? n) (fast-expt-iter (square b) (/ n 2) product))
        (else (fast-expt-iter b (- n 1) (* b product)))))

(define (square x) (* x x))

(define (even? n)
  (= (remainder n 2) 0))

1.17

1
2
3
4
5
6
7

(define (mul a b)
  (cond ((= b 0) 0)
        ((even? b) (mul (double a) (halve b)))
        (else (+ a (mul a (- b 1))))))

(define (double x) (* x 2))
(define (halve x) (/ x 2))

1.18

1
2
3
4
5
6
7
8
9
10
11

(define (mul a b)
  (mul-iter a b 0))

(define (mul-iter a b sum)
  (cond ((= b 0) sum)
        ((even? b) (mul-iter (double a) (halve b) sum))
        (else (mul-iter a (- b 1) (+ a sum)))))


(define (double x) (* x 2))
(define (halve x) (/ x 2))

1.19

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

; use the matrix multiplication to represent the transformation

;                       0 1
; (fib(n-1), fib(n)) *        = (fib(n), fib(n-1) + fib(n)) = (fib(n), fib(n+1))
;                       1 1

; a <- a + b
; b <- a

; the state transformation above is just a special case of the below one when
; p = 0 and q = 1

; a <- bq + aq + ap
; b <- bp + aq

(define (fib n)
  (fib-iter 1 0 0 1 n))

(define (fib-iter a b p q count)
  (cond ((= count 0) b)
        ((even? count)
         (fib-iter
           a
           b
           (+ (square p) (square q))
           (+ (square q) (* p q 2))
           (/ count 2)))
  (else (fib-iter
          (+ (* b q) (* a q) (* a p))
          (+ (* b p) (* a q))
          p
          q
          (- count 1)))))

(define (square x) (* x x))

1.21

1
2
3
4
5
6
7
8

; the following expressions produce
; 199
; 1999
; 7

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

1.25 Answer:

Since Scheme has built-in support for arbitrary precision arithmetic, the procedure will produce the same result as the original expmod, however, it will be very inefficient since the huge number arithmetic will take much longer time than the numbers can be represented by a single computer word.

The original expmod used the successive squaring, the numbers to be processed will never be larger than m^2.

1.26 Answer:

With the calling of square, the original problem can be reduced to a sub problem with half of the size at each of the step when even? test is true. So T(n) = T(n/2) = O(logn)

However, if the explicit multiplication used instead, the recursive call of expmod will be evaluated twice, it not only just compute the sub problems two time, it is actually a tree recursion like the first solution for computing fibnacci sequence, so the number of expmod calls grow exponentially, which conclues that T(n) = O(2^n * logn) = O(n)

1.27

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

(define (fermat-test n)
  (fermat-test-iter n 2))

(define (fermat-test-iter n a)
  (cond ((= n a) true)
        ((not (= (expmod a n n) a)) false)
        (else (fermat-test-iter n (+ a 1)))))

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

(define (square x) (* x x))


; Carmichael numbers
(fermat-test 561)
(fermat-test 1105)
(fermat-test 1729)
(fermat-test 2465)
(fermat-test 2821)
(fermat-test 6601)
; these all give #t

(prime? 561)
(prime? 1105)
(prime? 1729)
(prime? 2465)
(prime? 2821)
(prime? 6601)
; these all give #f

1.28

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38

(define (miller-rabin-test n)
  (define (try-it a)
    (= (expmod-with-signal a (- n 1) n) 1))
  (try-it (+ 2 (random (- n 2)))))

(define (expmod-with-signal base exp m)
  (cond ((= exp 0) 1)
        ((even? exp)
         (square-with-signal (expmod-with-signal base (/ exp 2) m) m))
        (else
         (remainder (* base (expmod-with-signal base (- exp 1) m)) m))))

(define (square-with-signal a n)
  (if (and (not (= a 1))
           (not (= a (- n 1)))
           (= (remainder (* a a) n) 1))
      0
      (remainder (* a a) n)))


(miller-rabin-test 561)
(miller-rabin-test 1105)
(miller-rabin-test 1729)
(miller-rabin-test 2465)
(miller-rabin-test 2821)
(miller-rabin-test 6601)
; these will all give #f with quite a chance, but with enough runs, it will
; have some false positive

(newline)

(miller-rabin-test 7)
(miller-rabin-test 23)
(miller-rabin-test 103)
(miller-rabin-test 1009)
(miller-rabin-test 1019)
(miller-rabin-test 1000033)
; these will always all give #t

SICP Section 1.1

Computational processes are abstract beings that inhabit computers. As they evolve, processes manipulate other abstract things called data. The evolution of a process is directed by pattern of rules called a program.

Why use Lisp for the book?

The most significant unique feature of Lisp is the fact that Lisp descriptions of processes, called procedures, can themselves be represented and manipulated as Lisp data. The importance of this is that there are powerful program-design techniques that rely on the ability to blur the traditional distinction between "passive" data and "active" processes.

The three mechanisms of combining simple ideas to form more complex ideas in any powerful programming languages:

  1. primitive expressions
  2. means of combination
  3. means of abstraction

Substitution model

  • Applicative order

    Evaluate the operator and operands first and then applies the resulting procedure to the resulting arguments. "Fully expand and then reduce"

  • Normal order

    Don't evaluate the operands until their values are needed, instead, substitute operand expressions for parameters until it obtained an expression involving only primitive operators, and would then perform the evaluation. "Evaluate the arguments and then apply"

Example: Square Roots by Newton's Method

The contrast between function and procedure is a reflection of the general distinction between describing properties of things and describing how to do things, or, as it is sometimes refered to, the distinction between declarative knowledge and imperative knowledge. In mathematics we are usually concerned with declarative (what is) descriptions, whereas in computer science we are usually concerned with imperative (how to) descriptions.

Exercises

1.3 Define a procedure that takes three numbers as arguments and return the sum of squares of the two larger numbers.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

(define (square x) (* x x))

(define (sum-of-squares a b) (+ (square a) (square b)))

(define (min a b)
  (if (< a b)
     a
     b))

(define (max a b)
  (if (< a b)
     b
     a))

(define (sum-of-squares-of-2-largest a b c)
  (sum-of-squares (max a b) (max c (min a b))))

1.4 Observe that our model of evaluation allows for combinations whose operators are compound expressions. Use this observation to describe the behavior of the following procedure:

1
2

(define (a-plus-abs-b a b)
  ((if (> b 0) + -) a b))

if b is greater than 0, the operator that will apply to a and b is +, or else it will be -, so the a-plus-abs-b will always result in a plus abs of b.

1.5 Answer: the interpreter that uses applicative order evaluation will hang due to the infinite recursive call of 'p', while an interpreter that uses normal order evaluation will get 0.

Note: p is a function, (p) is a call to function p.

1.6 Answer: the interpreter will hang due to the infinite recursive call to sqrt-iter. Since List uses applicative order evaluation, in the definition of new-if, the else-clause will always be evaluated no matter the result of the predicate, thus lead to infinite recursive call to sqrt-iter. That's why if needs to be a special form, the predicate expression is evaluated first, and the result determines whether to evaluate the consequent or the alternative expression.

1.7 Answer: For small values, the absolute tolerance 0.001 is too large, so the results become inaccurate. For example, (sqrt 0.001) gives 0.04124542607499115 on my machine. (ubuntu 12.10 x86_64); And for large values, due to the precision limitation of float-point representation, the guess couldn't be refined to a value that can be represented within the tolerance. In such cases, the program can endlessly alternate between two guesses that are more than 0.001 away from the true square root.

so, instead of using absolute tolerance, we changed to use the relative tolerance of two continuous guess values. This can be demonstrated with the below updated good-enough? procedure:

1
2

(define (good-enough? guess x)
  (< (abs (- (improve guess x) guess)) 0.001))

1.8

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

(define (cbrt-iter guess x)
  (if (good-enough? guess x)
    guess
    (cbrt-iter (improve guess x) x)))

(define (improve guess x)
  (/ (+ (/ x (square guess)) (* 2 guess)) 3))

(define (good-enough? guess x)
  (< (abs (- (improve guess x) guess)) 0.001))

(define (square x)
  (* x x))

(define (cbrt x)
  (cbrt-iter 1.0 x))

The only difference between cbrt and sqrt is the improve procedure.

Coursera Online Functional Programming Course - a Retrospective

Here comes the summary of the Coursera course Fucntional Programming Principles in Scala It's a little bit late, it has been two months since I finished the course. This is the first Coursera course I followed from the very beginning to the end and accomplished all the programming assignments with full scores, which helped me to get the certificate with distinction.

/images/posts/coursera_scala_certificate.jpg

First the excellent parts of the course:

  • Martin Odersky is a pretty good teacher, and all the lectures are well designed. I enjoyed watching the lectures a lot and most of which I have read more than once.
  • The programming assignments are well structured, with detailed instructions and step by step guide. In each of the assignments, the whole task is split into several steps with approprivate level of abstractions.
  • There are also many useful supporting materials like scala style guide, sbt tutorial, instructions on tools setup for the course, which made it easier to concentrate on the course content. For me, I don't have scala develop environment setup before the course, and by following the tools setup section in only took me less than half an hour to get everything ready for trying the example code and the assignments.

As for the less-good things, I would say:

  • Both the lectures and the assignments are quite easy for an experienced programmers (No functional programming background needed), I was expecting more challenging stuff.
  • There are no official solutions distributed. (The course will be offered again some time later) However, I still think for those students who passed the course should be qualified to get the solutions so that they can compare those with their own to see where they can still improve.
  • It's a pity that this is only the first half of an advanced undergraduate course that Martin taught on campus. I am interested in the other half.

After taking this cousre, I got a deeper understanding of the functional programming basics and it made me feel more comfortable while picking up SICP again (after 3 years). Now, I am convinced that I am able to go through SICP and finish most of the exercises. Also, I had a firmer grasp of Scala even though I didn't write more than 100 lines of Scala before; I understand more about the scala syntax, idioms and even the motivations behind some of the language structures. E.g., call by name/value, lazy evaluation, currying, pattern matching and so on. I will publish my detailed notes on the lectures and assignments to this blog later.

To conclude, I really enjoyed taking the course and many thanks to Martin , the TAs, and also the coursera staff for offering such a wonderful course.

Python Abc

Introduction

This module provides the infrastructure for defining abstract base classes (ABCs) in Python. The ABCs define a minimal set of methods that establish the characteristic behavior of the type. For more details about this, see PEP 3119.

Highlights

The module provides a metaclass used to create ABCs. An ABC can be subclassed directly. The class also has a 'register' method to register unrelated concrete classes (including built-in classes) and unrelated ABCs as 'virtual subclasses'

Also there are two decorators abstractmethod and abstractproperty, which will set the function object's attribute '__isabstractmethod__' to True. Only when all of the abstract methods and abstract properties are overriden, can a class that has a metaclass derived from ABCMeta be instantiated.

Code comments

ABCMeta.__new__

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

def __new__(mcls, name, bases, namespace):
    cls = super(ABCMeta, mcls).__new__(mcls, name, bases, namespace)
    # Compute set of abstract method names
    abstracts = set(name
                 for name, value in namespace.items()
                 if getattr(value, "__isabstractmethod__", False))
    for base in bases:
        for name in getattr(base, "__abstractmethods__", set()):
            value = getattr(cls, name, None)
            if getattr(value, "__isabstractmethod__", False):
                abstracts.add(name)
    cls.__abstractmethods__ = frozenset(abstracts)

    # chaoc: caches are the typical usages of weak references

    # Set up inheritance registry
    cls._abc_registry = WeakSet()
    cls._abc_cache = WeakSet()
    cls._abc_negative_cache = WeakSet()
    cls._abc_negative_cache_version = ABCMeta._abc_invalidation_counter
    return cls
  • It first creates a 'type' object cls. (super(ABCMeta, mcls) is 'type')
  • Iterate through all the attributes (including all the attributes inherited from all the bases), if any of them have '__isabstractmethod__' set to true, add it to cls's __abstractmethods__.
  • Initialize the attributes '_abc_registry', '_abc_cache', '_abc_negative_cache' and '_abc_negative_cache_version', which are used to speed up the check in __instancecheck__ and __subclasscheck__.

ABCMeta.register(cls, subclass)

See the added comments in line

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

def register(cls, subclass):
    """Register a virtual subclass of an ABC."""
    # chaoc: sanity check to make sure that subclass is class type
    # either a type object (for new style classes) or the type is
    # the same as ClassType(for old style classes)
    if not isinstance(subclass, (type, types.ClassType)):
        raise TypeError("Can only register classes")
    if issubclass(subclass, cls):
        return  # Already a subclass
    # Subtle: test for cycles *after* testing for "already a subclass";
    # this means we allow X.register(X) and interpret it as a no-op.
    if issubclass(cls, subclass):
        # This would create a cycle, which is bad for the algorithm below
        raise RuntimeError("Refusing to create an inheritance cycle")
    cls._abc_registry.add(subclass)
    ABCMeta._abc_invalidation_counter += 1  # Invalidate negative cache

ABCMeta.__instancecheck__

See the added comments in line

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

def __instancecheck__(cls, instance):
    """Override for isinstance(instance, cls)."""
    # Inline the cache checking when it's simple.
    subclass = getattr(instance, '__class__', None)
    if subclass is not None and subclass in cls._abc_cache:
        return True
    subtype = type(instance)
    # Old-style instances
    if subtype is _InstanceType:
        subtype = subclass
    # chaoc: subtype will also be subclass for old style classes
    # as assigned in the above step
    if subtype is subclass or subclass is None:
        # chaoc: check if the negative cache is safe to use or not
        if (cls._abc_negative_cache_version ==
            ABCMeta._abc_invalidation_counter and
            subtype in cls._abc_negative_cache):
            return False
        # Fall back to the subclass check.
        return cls.__subclasscheck__(subtype)
    return (cls.__subclasscheck__(subclass) or
            cls.__subclasscheck__(subtype))

ABCMeta.__subclasscheck__

The code and comment in this function is very clear and straightforward.

Just make sure the different cases needed to check:

  1. check the subclass hook
  2. check if it's a direct subclass through __mro__
  3. check if it's a subclass of a registered class (issubclass is called to do recursive check)
  4. check if it's a subclass of a subclass (issubclass is called to do recursive check)

In this post, we only talk about the defitions of ABCMeta. We will see the typical usages in the collections module.

Python Bisect

Introduction

This module provides support for maintaining a list in sorted order without having to sort the list after each insertion. It uses a basic bisection algorithm similar with the classic binary search.

Code comments

The above binsearch function will return the first of the elements that equals to the target if there are more than one. Sometimes it's required to find the left most one or the right most one. We can achieve this by using the bisect() functions. Let's examine the pre and post conditions, and also loop invariants.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

def bisect_right(a, x, lo=0, hi=None):

    # Let l be the original lo passed in, and h be hi. In order to differentiate.

    # @requires 0 <= l && l < h && h <= len(a)
    # @requires is_sorted(a[l:h])
    # @ensures all e in a[l:result] have e <= x, and all e in a[result:h] have e > x

    while lo < hi:
        # @loop_invariant l <= lo && lo <= hi && hi <= h
        # @loop_invariant lo == l || a[lo - 1] <= x
        # @loop_invariant hi == h || a[hi] > x
        mid = (lo + hi) // 2
        # @assert lower <= mid && mid < higher
        if x < a[mid]:
            hi = mid
        else:
            lo = mid + 1
    return lo

After exiting the loop, we have lo == hi. Now we have to distinguish some cases:

  • If lo == l, from the third conjunct, we know that a[hi] > x, and since lo == hi, we have a[lo] > x. In this case, all e in a[l:h] > x
  • If hi == h, from the second conjunct, we know that a[lo - 1] < x, in this case, all e in a[l:h] <= x
  • if lo != l and hi != h, from the second and third conjuncts, we know that a[lo - 1] <= x, a[hi] > x. The post condition still holds.

We can do the same analysis on the function bisect_left().

Note: In the pre condition I explicitly written down that lo < hi is required, the code will directly return lo when hi is less than or equal to lo, but that's simply meaningless, if this is the case, the sorted order of 'a' after we insert element before the index returned by the function.

About the function prototype

In the doc of the module, it's bisect.bisect_right(a, x, lo=0, hi=len(a)) while in the source code, it's def bisect_left(a, x, lo=0, hi=None)

bisect.bisect_right(a, x, lo=0, hi=len(a)) is not valid python code, you will get error like this: NameError: name 'a' is not defined. This is because default values are computed and bound at the function definition time rather than when you call the function. This means that you can't have a default which is dependent on something that is not known until the function is called.

Pythonic stuff

Override function definitions in Python

1
2
3
4
5

# Overwrite above definitions with a fast C implementation
try:
    from _bisect import *
except ImportError:
    pass

Gotcha -- Mutable default arguments

See here for more details.

Support LaTex in Octopress

It took me some time to finally get latex math formulas working in Octopress. If you googled ‘Octopress latex’, you can get quite a few online resources about how to support latex in octopress, with various levels of complexity. In this post, I will write down how I achieve this.

The initial attempt

As I installed the jekyll-rst plugin to use rst to write my posts, I thought it should be easy to write latex math because docutils has native support for it since version 0.8 (A :math: role and also a .. math: directive introduced for that). However, after I tried to use these in a octopress post, I found that the post will be rendered to empty. For example, if I insert the following rst code into my post, the whole post becomes empty; but after removing this line, everything is fine.

1

The area of a circle is :math:`A_\text{c} = (\pi/4) d^2`.

I also verified that the exact same code can be successfully converted to valid html using the ‘rst2html.py’ script on my system, so I guess maybe something is wrong in ‘RbST’. I found that in RbST, it has its own copies of rst2html and rst2latex tools under /gems/RbST-0.1.3/lib/rst2parts ,

1
2
3
4

chuchao@chuchao:~/.rvm/gems/ruby-1.9.2-p320/gems/RbST-0.1.3/lib/rst2parts$ ls
__init__.py  rst2html.py  rst2latex.py  transform.py  transform.pyc
chuchao@chuchao:~/.rvm/gems/ruby-1.9.2-p320/gems/RbST-0.1.3/lib/rst2parts$ ls ..
rbst.rb  rst2parts

which will be used in rbst.rb. I have even tried to change rbst.rb to use the rst2html.py installed on my system, but this also didn’t get any luck.

1
2
3
4
5
6
7
8
9

class RbST

  @@executable_path = File.expand_path(File.join(File.dirname(__FILE__), "rst2parts"))
  @@executables = {
    # :html  => File.join(@@executable_path, "rst2html.py"),
    :html  => "/usr/local/bin/rst2html.py",
    :latex => File.join(@@executable_path, "rst2latex.py")
  }
...

Finally, I gave up on this and opened an issue on this for the jekyll-rst plugin. Hope the author can fix this.

Switch back to markdown and using kramdown

After a google search about this issue, seems that the simplest solution is to use kramdown, which has built-in support for latex.

First, install kramdown: gem install karmdown

Then, add mathjax configs into tag, in octopress, just add the below code into /source/_includes/custom/head.html:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

<!-- mathjax config similar to math.stackexchange -->

<script type="text/x-mathjax-config">
  MathJax.Hub.Config({
    tex2jax: {
      inlineMath: [ ['$','$'], ["\\(","\\)"] ],
      processEscapes: true
    }
  });
</script>

<script type="text/x-mathjax-config">
    MathJax.Hub.Config({
      tex2jax: {
        skipTags: ['script', 'noscript', 'style', 'textarea', 'pre', 'code']
      }
    });
</script>

<script type="text/x-mathjax-config">
    MathJax.Hub.Queue(function() {
        var all = MathJax.Hub.getAllJax(), i;
        for(i=0; i < all.length; i += 1) {
            all[i].SourceElement().parentNode.className += ' has-jax';
        }
    });
</script>

<script type="text/javascript"
   src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
</script>

See here for more details. After this, we are ready to test latex math in our post. For example:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

$$
\begin{align*}
  & \phi(x,y) = \phi \left(\sum_{i=1}^n x_ie_i, \sum_{j=1}^n y_je_j \right)
  = \sum_{i=1}^n \sum_{j=1}^n x_i y_j \phi(e_i, e_j) = \\
  & (x_1, \ldots, x_n) \left( \begin{array}{ccc}
      \phi(e_1, e_1) & \cdots & \phi(e_1, e_n) \\
      \vdots & \ddots & \vdots \\
      \phi(e_n, e_1) & \cdots & \phi(e_n, e_n)
    \end{array} \right)
  \left( \begin{array}{c}
      y_1 \\
      \vdots \\
      y_n
    \end{array} \right)
\end{align*}
$$

will render as

$$ \begin{align} & \phi(x,y) = \phi \left(\sum{i=1}^n x_ie_i, \sum{j=1}^n y_je_j \right) = \sum{i=1}^n \sum{j=1}^n x_i y_j \phi(e_i, e_j) = \ & (x_1, \ldots, x_n) \left( \begin{array}{ccc} \phi(e_1, e_1) & \cdots & \phi(e_1, e_n) \ \vdots & \ddots & \vdots \ \phi(e_n, e_1) & \cdots & \phi(e_n, e_n) \end{array} \right) \left( \begin{array}{c} y_1 \ \vdots \ y_n \end{array} \right) \end{align} $$

And for inline latex code, just use $\exp(-\frac{x^2}{2})$, which will give $\exp(-\frac{x2}{2})$.

Python Heapq

Introduction

The module provides an implementation of heap queue algorithm, also known as priority queue algorithm.

Highlights

  • Zero-based indexing is used, so the children's index of node with index k are (2*k + 1) and (2*k + 2) respectively.
  • Internally a 'min heap' is maintained rather than 'max heap', which is more generally used in algorithm textbooks.
  • Three general functions based on heaps are also provided:
1
2
3
4
5
6
7
8
9
10

# Merge multiple sorted inputs into a single sorted output
# (for example, merge timestamped entries from multiple log files)
heapq.merge(*iterables)

# The following two functions are effectively like
# sorted(iterable, key=key, reverse=True)[:n] and
# sorted(iterable, key=key)[:n],
# but they perform best with smaller values of n
heapq.nlargest(n, iterable[, key])
heapq.nsmallest(n, iterable[, key]))

Pythonic stuff

Rich comparison methods

1
2
3

def cmp_lt(x, y):
    # use __lt__ if available; otherwise, try __le__
    return x < y if hasattr(x, '__lt__') else (not y <= x)

In python, there are 6 so called "Rich Comparison" methods, x < y calls x.__lt__(y) and others are similar (__le__ and <=; __gt__ and >=; __eq__ and ==; __ne__ and <>). Arguments to rich comparison methods are never coerced. see coercion

Code comments

Why heapify(x) is O(n)?

This is not obvious at first by seeing the code, given that there is a 'while' loop in _siftup and also a while loop in _siftdown(called in _siftup). Let's look into it further:

  1. in the while loop of _siftup, it takes O(L) time for nodes that L levels above leaves.
  2. and in the while loop of _siftdown called in _siftup, it takes at most L steps, so _siftdown is O(L).
  3. since we have n/4 nodes in level 1, n/8 nodes in level 2, and finally one root node, which is lg(n) levels above leaf, so the total amount in the while loop of heapify is:
n/4 * c + n/8 * c + n/16 * 3c + ... + 1 * lg(n) * c, and let n/4 = 2^k, after simplification, we get:
c * 2^k(1/2^0 + 2/2^1 + 3/2^2 + ... + (k+1)/2^k), as the limit of (k+1)/2^k is 0 when k is infinite, so the term in the brackets bound to a constant, from this we can conclude that heapify is O(2^k), which is O(n).

Why it continues to find the smaller child until a leaf is hit in _siftup?

As explained in the comment by the module author, this is a ad hoc to reduce the comparisons on the following operations on the heap.

Introduction on Reading Source of Python Standard Libraries Series

I am going to start the series of posts on reading the source of python standard modules. I will go with the pure python modules first, and maybe later I can continue with C implementations of the modules. Let's see how far I could go.

What will be included

  • A brief introduction of the module. (It should be very short, people can go to the standard library doc for more information.)
  • Special highlights about the important APIs, implementation details.
  • Python features/idioms/tricks/gotchas that worth the whistle, especially those I was not familiar with
  • Detail explanations about the tricky part of the code if any

What will not be included

Also, alone the way, I may start another series on some specific 'advanced topics' in python, like descriptor, decorator, method resolution order(mro) and so on. Mainly about why they are introduced into python, how they are used and the typical use cases. This is inspired by the blogs about python history

This post will also be used to track my progress.