math « Dan's Thoughts

25Jul/090

Exercise 2.9 of SICP

Exercise 2.9: The width of an interval is half of the difference between its upper and lower bounds. The width is a measure of the uncertainty of the number specified by the interval. For some arithmetic operations the width of the result of combining two intervals is a function only of the widths of the argument intervals, whereas for others the width of the combination is not a function of the widths of the argument intervals. Show that the width of the sum (or difference) of two intervals is a function only of the widths of the intervals being added (or subtracted). Give examples to show that this is not true for multiplication or division.

$(a-\epsilon_1,a+\epsilon_1) \pm (b-\epsilon_2,b+\epsilon_2)$
$w_a=\frac{1}{2}(a+\epsilon_1-(a-\epsilon_1))=\epsilon_1$
$w_b=\frac{1}{2}(b+\epsilon_2-(b-\epsilon_2))=\epsilon_2$

Addition:

$(a-\epsilon_1,a+\epsilon_1)+(b-\epsilon_2,b+\epsilon_2)$
$\left((a+b)-(\epsilon_1+\epsilon_1),(a+b)+(\epsilon_1+\epsilon_1))\right)$
$w_{a+b}=\epsilon_1+\epsilon_2=w_a+w_b$

Subtraction:

$(a-\epsilon_1,a+\epsilon_1)-(b-\epsilon_2,b+\epsilon_2)$
$\left((a-b)-(\epsilon_1+\epsilon_2),(a-b)+(\epsilon_1+\epsilon_2)\right)$
$w_{a-b}=\epsilon_1+\epsilon_2= w_a+w_b$

Multiplication:

$(a-\epsilon_1,a+\epsilon_1) \cdot (b-\epsilon_2,b+\epsilon_2)$

$(a+\epsilon_1)\cdot(b+\epsilon_2)=ab+a\epsilon_2+b\epsilon_1+\epsilon_1\epsilon_2$
$(a+\epsilon_1)\cdot(b-\epsilon_2)=ab-a\epsilon_2+b\epsilon_1-\epsilon_1\epsilon_2$
$(a-\epsilon_1)\cdot(b+\epsilon_2)=ab+a\epsilon_2-b\epsilon_1-\epsilon_1\epsilon_2$
$(a-\epsilon_1)\cdot(b-\epsilon_2)=ab-a\epsilon_2-b\epsilon_1+\epsilon_1\epsilon_2$

If any of the four intervals expressions are used as interval endpoints their width formula will not give the tidy results of addition/subtraction.

Division:

$(a-\epsilon_1,a+\epsilon_1) \cdot (\frac{1}{b+\epsilon_2},\frac{1}{b-\epsilon_2})$

$(a-\epsilon_1)\cdot\frac{1}{b+\epsilon_2}=\frac{a-\epsilon_1}{b+\epsilon_2}$
$(a+\epsilon_1)\cdot\frac{1}{b+\epsilon_2}=\frac{a+\epsilon_1}{b+\epsilon_2}$
$(a-\epsilon_1)\cdot\frac{1}{b-\epsilon_2}=\frac{a-\epsilon_1}{b-\epsilon_2}$
$(a+\epsilon_1)\cdot\frac{1}{b-\epsilon_2}=\frac{a+\epsilon_1}{b-\epsilon_2}$

Examples:

w₁=0.68
w₂=0.235
w_total=0.915
> (print (add-interval (make-interval 6.12 7.48) (make-interval 4.465 4.935)))
[10.585,12.415]
w_actual=0.915

w₁=0.68
w₂=0.235
w_total=0.915
> (print (sub-interval (make-interval 6.12 7.48) (make-interval 4.465 4.935)))
[1.1850000000000005,3.0150000000000006]
w_actual=0.915 (ignoring the floating point number errors)

w₁=0.68
w₂=0.235
w_{naive_total}=0.15980
Of course this doesn't make sense since accuracy cannot be magically gained by multiplying two interval together. More on this in the division example.
> (print (mul-interval (make-interval 6.12 7.48) (make-interval 4.465 4.935)))
[1.1850000000000005,3.0150000000000006]
w_actual=4.794

w₁=0.68
w₂=0.235
w_{naive_total}=2.89361702127659574468
> (print (div-interval (make-interval 6.12 7.48) (make-interval 4.465 4.935)))
[1.2401215805471126,1.6752519596864504]
w_actual=.2175651895696689
How is it that w_actual is smaller than each w₁, w₂? w_{naive_total} is just plain wrong and should be discarded.

The answer is that w_actual isn't smaller. It's actually larger. The percentage of error that w_actual represents, for it's interval, is actually greater than either w₁=10% and w₂=5%. w_actual = 14.9%.
$6.8 \pm 0.68$
$0.68 = \frac{0.68}{6.8}=10\%$
$4.7 \pm 0.235$
$0.235 = \frac{0.235}{4.7}=5\%$
$0.2175651895696689 = \frac{.2175651895696689}{\left(\frac{1}{2}\cdot(1.2401215805471126+1.6752519596864504)\right)}=14.9\%$

As expected the error is actually higher even though the aggregate number for w_actual looks smaller.

Filed under: SICP, lisp, math No Comments

24Jul/090

Exercise 2.8 of SICP

Exercise 2.8: Using reasoning analogous to Alyssa's, describe how the difference of two intervals may be computed. Define a corresponding subtraction procedure, called sub-interval.

In the resistor example, the resistor is measured at 6.8 ohm with a tolerance of 10%. This means the actual value is $6.8 \pm 0.68$ , which makes the interval [6.12 , 7.48]. The other resistor is $4.7 \pm 0.235$ with the interval of [4.465 , 4.935]. Naively subtracting the two intervals would give, [1.655 , 2.545]. Which doesn't seem wrong but it is. It is wrong because we are biasing the tolerance in one direction for both measurements. We assume that the first value is $6.8-0.68$ and second $4.7-0.235$ . There is no basis to assume that the following case isn't just as valid: $6.8-0.68$ and $4.7+0.235$ .
It's even clearer with the interval of $3 \pm 1$ [2 , 4] and $2 \pm 1$ [1 , 3]. Naively subtracting these terms yields an interval of 0 length, [1,1] or $1 \pm 0$ which is nonsense.
In order to make sure no tolerance bias exists that makes calculations seem more accurate than is possible, the worst case measurement error must be assumed. The intervals will become as wide as possible to not give a false sense of accuracy.

(define (sub-interval x y)
  (make-interval (- (lower-bound x) (upper-bound y))
                 (- (upper-bound x) (lower-bound y))))
(define (make-interval a b) (cons a b))
 
(define (lower-bound inter) (car inter))
(define (upper-bound inter) (cdr inter))
 
(define (add-interval x y)
  (make-interval (+ (lower-bound x) (lower-bound y))
                 (+ (upper-bound x) (upper-bound y))))
(define (add-interval x y)
  (make-interval (+ (lower-bound x) (lower-bound y))
                 (+ (upper-bound x) (upper-bound y))))
(define (div-interval x y)
  (mul-interval x 
                (make-interval (/ 1.0 (upper-bound y))
                               (/ 1.0 (lower-bound y)))))
(define (mul-interval x y)
  (let ((p1 (* (lower-bound x) (lower-bound y)))
        (p2 (* (lower-bound x) (upper-bound y)))
        (p3 (* (upper-bound x) (lower-bound y)))
        (p4 (* (upper-bound x) (upper-bound y))))
    (make-interval (min p1 p2 p3 p4)
                   (max p1 p2 p3 p4))))
 
(define (print int)
  (newline)
  (display "[")
  (display (lower-bound int))
  (display ",")
  (display (upper-bound int))
  (display "]")
  (newline))

> (print (sub-interval (make-interval 6.12 7.48) (make-interval 4.465 4.935)))
[1.1850000000000005,3.0150000000000006]

Filed under: SICP, lisp, math No Comments

19Jul/090

Exercise 2.3 of SICP

Exercise 2.3: Implement a representation for rectangles in a plane. (Hint: You may want to make use of exercise 2.2.) In terms of your constructors and selectors, create procedures that compute the perimeter and the area of a given rectangle. Now implement a different representation for rectangles. Can you design your system with suitable abstraction barriers, so that the same perimeter and area procedures will work using either representation?

I have decided to use a diagonal segment for representation of make-rectangle and a point with both dimensional magnitudes for make-rectangle2.

(define (make-segment p1 p2)
  (cons p1 p2))
(define (start-segment s)
  (car s))
(define (end-segment s)
  (cdr s))
(define (average a b) (/ (+ a b) 2.0))
(define (make-point x y)
    (cons x y))
(define (x-point p)
  (car p))
(define (y-point p)
  (cdr p))
(define (midpoint-segment s)
  (make-point
    (average
        (x-point (start-segment s))
        (x-point (end-segment s)))
    (average
        (y-point (start-segment s))
        (y-point (end-segment s)))))
 
(define (dimension direction r)
  (- (direction (end-segment r))
     (direction (start-segment r))))
 
(define (width r) (dimension y-point r))
(define (length r) (dimension x-point r))
 
(define (make-rectangle s) s)
 
(define (make-rectangle2 p length width)
  (make-segment p
                (make-point
                  (+ (x-point p) length)
                  (+ (y-point p) width))))
 
(define (area r)
  (* (length r)
     (width r)))
 
(define (perimeter r)
  (* 2 (+ (length r)
          (width r))))

(define R
  (make-rectangle (make-segment
                    (make-point 0 0)
                    (make-point 4 4))))

> (area R)
16
> (perimeter R)
16

(define R2 (make-rectangle2
             (make-point 0 0)
             4 4))

> (area R2)
16
> (perimeter R2)
16

(define R
  (make-rectangle (make-segment
                    (make-point 2 1)
                    (make-point 5 4))))

> (area R)
9
> (perimeter R)
12

(define R2 (make-rectangle2
             (make-point 2 1)
             3 3))

> (area R)
9
> (perimeter R)
12

Filed under: SICP, lisp, math No Comments

15Jul/090

Exercise 1.45 of SICP

Exercise 1.45: We saw in section 1.3.3 that attempting to compute square roots by naively finding a fixed point of $y \mapsto \frac{x}{y}$ does not converge, and that this can be fixed by average damping. The same method works for finding cube roots as fixed points of the average-damped $y \mapsto \frac{x}{y^2}$ . Unfortunately, the process does not work for fourth roots -- a single average damp is not enough to make a fixed-point search for $y \mapsto \frac{x}{y^3}$ converge. On the other hand, if we average damp twice (i.e., use the average damp of the average damp of $y \mapsto \frac{x}{y^3}$ ) the fixed-point search does converge. Do some experiments to determine how many average damps are required to compute nth roots as a fixed-point search based upon repeated average damping of $y \mapsto \frac{x}{y^{n-1}}$ . Use this to implement a simple procedure for computing n^th roots using fixed-point, average-damp, and the repeated procedure of exercise 1.43. Assume that any arithmetic operations you need are available as primitives.

Doing experiments on the average number of damps to get the fixed-point procedure to converge seems to be:

Degree of root	Power of 2	Number of average damps
2	2¹	1
4	2²	2
8	2³	3
16	2⁴	4
32	2⁵	5

Using this basis it makes sense to write a procedure discrete-log that computes the number of powers of 2 in a given number. It's a bit like a very primitive implementation log₂ function hence the name. Discrete-log is $\Theta(\log_2 n)$ because it halves the number of computations every iteration.

(define (root num degree)
  (define (discrete-log n)
    (if (< n 2)
      0
      (+ 1 (discrete-log (/ n 2)))))
  (define (compose f g)
    (lambda (x) (f (g x))))
  (define (repeated-compose f n)
    (if (< n 2)
      f
      (compose f (repeated-compose f (- n 1)))))
  (define (average-damp f)
    (lambda (x) (/ (+ x (f x)) 2)))
  (let ((number-of-compositions (discrete-log degree))
        (initial-function (lambda (x) (/ num (expt x (- degree 1))))))
    (cond ((= degree 0) 1)
          ((= degree 1) num)
          (else
            (fixed-point 
              ((repeated-compose average-damp number-of-compositions) initial-function)
              1.0)))))
(define (fixed-point f guess)
  (define error 0.001)
  (define (close-enough? x y) (< (abs (- 1 (/ x y))) error))
  (define (try guess)
    (let ((new-guess (f guess)))
      (if (close-enough? guess new-guess)
        new-guess
        (try new-guess))))
  (try guess))

> (root 2 0)
1
> (root 2 1)
2
> (root 2 2)
1.4142135623746899
> (root (expt 2 50) 50)
1.9999956962054166
> (root 0.00000001 2)
1.0000000000082464e-4

This program summarizes most of Chapter 1. It has lexical scope, passing functions as parameters, functions that construct and return other functions and the square root algorithm with an improved close-enough? procedure from exercise 1.7.
The repeated-compose and discrete-log could be done iteratively instead of recursively like this.

(define (repeated-compose f n)
  (define (compose f g)
    (lambda (x) (f (g x))))
  (define (iter n func)
    (if (< n 2)
      func
      (iter (- n 1) 
            (compose f func))))
  (iter n f))

(define (discrete-log n)
  (define (iter n exponent)
    (if (< n 2)
      exponent
      (iter (/ n 2) (+ 1 exponent))))
  (iter n 0))

Filed under: SICP, Uncategorized, algorithms, lisp, math No Comments

14Jul/090

Exercise 1.44 of SICP

Exercise 1.44: The idea of smoothing a function is an important concept in signal processing. If f is a function and dx is some small number, then the smoothed version of f is the function whose value at a point x is the average of f(x - dx), f(x), and f(x + dx). Write a procedure smooth that takes as input a procedure that computes f and returns a procedure that computes the smoothed f. It is sometimes valuable to repeatedly smooth a function (that is, smooth the smoothed function, and so on) to obtained the n-fold smoothed function. Show how to generate the n-fold smoothed function of any given function using smooth and repeated from exercise 1.43.

(define (compose f g)
  (lambda (x) 
    (f (g x))))
(define (repeated fn n)
  (if (= n 1)
    fn
    (compose 
      fn (repeated fn (- n 1)))))
(define (smooth f)
  (define (average a b c)
    (/ (+ a b c) 3))
  (let ((dx 0.00001))
  (lambda (x)
    (average (f (- x dx)) (f x) (f (+ x dx))))))
 
(define (n-fold-smooth f n)
    (if (not (> n 0))
      f
      (let ((smoothed-n (repeated smooth n)))
        (smoothed-n f))))
 
(define (S x)
  (define (>= a b) (not (< a b)))
  (cond ((< x 0) 0.0)
        ((>= x 0) 1.0)))

S is a modified (the argument is not floored) Heaviside step function.
> ((n-fold-smooth S 0) 0)
1
> ((n-fold-smooth S 1) 0)
.6666666666666666
> ((n-fold-smooth S 15) 0)
.5409601581500249

Filed under: SICP, lisp, math No Comments

10Jul/090

Exercise 1.40 of SICP

Exercise 1.40: Define a procedure cubic that can be used together with the newtons-method procedure in expressions of the form

 (newtons-method (cubic a b c) 1)

to approximate zeros of the cubic x³ + ax² + bx + c.

(define (cubic a b c)
  (lambda (x)
    (+
      (* x x x)
      (* a x x)
      (* b x)
      c)))
 
(define dx 0.00001)
(define (deriv g)
  (lambda (x)
    (/ (- (g (+ x dx)) (g x))
       dx)))
 
(define (newton-transform g)
  (lambda (x)
    (- x
       (/ (g x) ((deriv g) x)))))
 
(define (newtons-method g guess)
  (fixed-point (newton-transform g) guess))
 
(define (fixed-point f first-guess)
  (define tolerance 0.00001)
  (define (close-enough? v1 v2)
    (&lt; (abs (- v1 v2)) tolerance))
  (define (try guess)
    (let ((next (f guess)))
      (if (close-enough? guess next)
          next
          (try next))))
  (try first-guess))

> (newtons-method (cubic -2 -9 18) -2)
2.000000000000876
> (newtons-method (cubic 8 3 6) 1)
-7.711875650348891

Note that Newton's Method only finds a root not all roots or even all real roots.

Filed under: SICP, lisp, math No Comments

9Jul/090

Exercise 1.39 of SICP

Exercise 1.39: A continued fraction representation of the tangent function was published in 1770 by the German mathematician J.H. Lambert:

$\tan x = \cfrac{x}{1-\cfrac{x^2}{3-\cfrac{x^2}{5-{\ddots\,}}}}$

where x is in radians. Define a procedure (tan-cf x k) that computes an approximation to the tangent function based on Lambert's formula. K specifies the number of terms to compute, as in exercise 1.37.

(define (tan-cf x k)
  (define (n i)
    (if (= i 1)
      x
      (* x x)))
  (define (d i)
    (- (* 2 i) 1))
  (define (cf i)
    (if (= i (+ k 1))
      0
      (/ (n i) 
         (- (d i) (cf (+ i 1))))))
  (if (not (> k 0))
    0
    (cf 1)))

> (tan-cf (sqrt 2) 12)
6.334119167042199
$\tan \sqrt 2 = 6.3341191670421915540568332642277$

Filed under: SICP, lisp, math No Comments

8Jul/090

Exercise 1.38 of SICP

Exercise 1.38: In 1737, the Swiss mathematician Leonhard Euler published a memoir De Fractionibus Continuis, which included a continued fraction expansion for e - 2, where e is the base of the natural logarithms. In this fraction, the N_i are all 1, and the D_i are successively 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, .... Write a program that uses your cont-frac procedure from exercise 1.37 to approximate e, based on Euler's expansion.

Since all the continued fraction works is done from the previous exercise the trickiest part is defining a function for the D_i's. Procedure d should always return 1 except for the special indices. For me it was key to notice that every 3i-1 index, where is is i=1...n, was the even number I care about. And that even number was of the form 2i.

(define n (lambda (i) 1.0))
(define (d i)
  (if (= 0 (remainder (+ i 1) 3))
    (* 2 (/ (+ i 1) 3))
    1))
 
(define (cont-frac n d k)
  (define (cf i)
    (if (= i (+ k 1))
      0
      (/ (n i) 
         (+ (d i) (cf (+ i 1))))))
  (if (not (> k 0))
    0
    (cf 1)))
 
(define (cont-frac-iter n d k)
  (define (iter k result)
    (if (= k 0)
      result
      (iter (- k 1) (/ (n k) (+ (d k) result)))))
  (iter k 0))

> (cont-frac n d 11)
.7182818352059925
> (cont-frac-iter n d 11)
.7182818352059925
$e-2 = 0.7182818284590452353$
$\epsilon = |e-2-0.7182818352059925| = 0.0000000067469472647$

Filed under: SICP, lisp, math No Comments

7Jul/090

Exercise 1.37 of SICP

Exercise 1.37:
a. An infinite continued fraction is an expression of the form

$f = \cfrac{N_1}{D_1 + \cfrac{N_2}{D_2 + \cfrac{N_3}{D_3 +\cdots}}}$

As an example, one can show that the infinite continued fraction expansion with the N_i and the D_i all equal to 1 produces $\frac{1}{\phi}$ , where $\phi$ is the golden ratio (described in section 1.2.2). One way to approximate an infinite continued fraction is to truncate the expansion after a given number of terms. Such a truncation -- a so-called k-term finite continued fraction -- has the form

$f = \cfrac{N_1}{D_1 + \cfrac{N_2}{{\ddots\,} + \cfrac{N_k}{D_k}}}$

Suppose that n and d are procedures of one argument (the term index i) that return the N_i and D_i of the terms of the continued fraction. Define a procedure cont-frac such that evaluating (cont-frac n d k) computes the value of the k-term finite continued fraction. Check your procedure by approximating $\frac{1}{\phi}$ using

(cont-frac (lambda (i) 1.0)
           (lambda (i) 1.0)
           k)

for successive values of k. How large must you make k in order to get an approximation that is accurate to 4 decimal places?

(define (cont-frac n d k)
  (define (cf i)
    (if (= i (+ k 1))
      0
      (/ (n i) 
         (+ (d i) (cf (+ i 1))))))
  (if (not (> k 0))
    0
    (cf 1)))

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

.6180555555555556
It takes 11 iterations to get to 4 decimal places within the reciprocal of golden ratio.
$\frac{1}{\phi} = 0.618033988$
It's amusing that the LaTeX looks almost identical to the above definition in scheme:

\cfrac{N_1}{D_1 + \cfrac{N_2}{D_2 + \cfrac{N_3}{D_3 +\cdots}}}

b. If your cont-frac 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.

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

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

.6180555555555556

Filed under: SICP, lisp, math No Comments

6Jul/090

Exercise 1.36 of SICP

Exercise 1.36: Modify fixed-point so that it prints the sequence of approximations it generates, using the newline and display primitives shown in exercise 1.22. Then find a solution to $x^x = 1000$ by finding a fixed point of $x \mapsto \frac{\log 1000}{\log x}$ . (Use Scheme's primitive log procedure, which computes natural logarithms.) Compare the number of steps this takes with and without average damping. (Note that you cannot start fixed-point with a guess of 1, as this would cause division by log(1) = 0.)

(define (average a b) (/ (+ a b) 2))
(define tolerance 0.00001)
(define (fixed-point f first-guess)
  (define (close-enough? v1 v2)
    (< (abs (- v1 v2)) tolerance))
  (define (try guess)
    (display guess)
    (newline)
    (let ((next (f guess)))
      (if (close-enough? guess next)
          next
          (try next))))
  (try first-guess))

> (fixed-point (lambda (x) (/ (log 1000) (log x))) 2)

2
9.965784284662087
3.004472209841214
6.279195757507157
3.759850702401539
5.215843784925895
4.182207192401397
4.8277650983445906
4.387593384662677
4.671250085763899
4.481403616895052
4.6053657460929
4.5230849678718865
4.577114682047341
4.541382480151454
4.564903245230833
4.549372679303342
4.559606491913287
4.552853875788271
4.557305529748263
4.554369064436181
4.556305311532999
4.555028263573554
4.555870396702851
4.555315001192079
4.5556812635433275
4.555439715736846
4.555599009998291
4.555493957531389
4.555563237292884
4.555517548417651
4.555547679306398
4.555527808516254
4.555540912917957
4.555532270803653

$4.555532270803653^{4.555532270803653} = 999.9913579312362$
35 Steps.

Average damping: $x = \frac{\log 1000}{\log x} = \frac{1}{2}\left( \frac{\log 1000}{\log x} + x\right)$
> (fixed-point (lambda (x) (average (/ (log 1000) (log x)) x)) 2)