math « Dan's Thoughts

28May/090

Exercise 1.12 of SICP

Exercise 1.12: The following pattern of numbers is called Pascal's triangle.

$\begin{tabular}{rccccccccc}row=0:& & & & & 1\\\noalign{\smallskip\smallskip}row=1:& & & & 1 & & 1\\\noalign{\smallskip\smallskip}row=2:& & & 1 & & 2 & & 1\\\noalign{\smallskip\smallskip}row=3:& & 1 & & 3 & & 3 & & 1\\\noalign{\smallskip\smallskip}row=4:& 1 & & 4 & & 6 & & 4 & & 1\\\noalign{\smallskip\smallskip}\end{tabular}$

The LaTeX for this triangle is really a pain however this great python script makes the process trivial.

The numbers at the edge of the triangle are all 1, and each number inside the triangle is the sum of the two numbers above it. Write a procedure that computes elements of Pascal's triangle by means of a recursive process.

(define (pascal row col)
  (cond ((and (= row 0) (= col 0)) 1)
        ((= col 0) 1)
        ((and (= row col) 1))
        (else
          (+ (pascal (- row 1) (- col 1))
             (pascal (- row 1) col)))))

Both 1 and 2 hold true for the top node.

The first condition recognizes that the first element on of every row is always 1, hence P(n,0) = 1 for each n
The second condition recognizes that the last element of every row is always 1, hence P(n,n) = 1 for each n
This final condition states that in order to get the next number in the sequence: go up a row and back one column + up a row and stay in the same column P(n,k)=P(n-1,k-1)+P(n-1,k)

> (pascal 0 0)
1
> (pascal 7 2)
21
> (pascal 16 8)
12870

Filed under: lisp, math, SICP No Comments

27May/090

Exercise 1.11 of SICP

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.

The idea is to use a triple of integers a,b,c initialized to f(2) = 2 f(1) = 1 and f(0) = 0, and to repeatedly apply the simultaneous transformations
$a \leftarrow a + 2b + 3c$
a is now f(n)
$b \leftarrow a$
b is now f(n-1)
$c \leftarrow b$
c is now f(n-2)

This is nothing but a slightly more complex Fibonacci iteration example from earlier.
F(n) = F(n-1) + F(n-2) where F(n) = n if n<2
$F_n = F_{n-1} + F_{n-2}$
$a \leftarrow a + b$
a is now F(n)
$b \leftarrow a$
b is now F(n-1)

Recursive implementation of f(n) is straightforward:

(define (f-rec n)
  (if (< n 3)
    n
    (+ (f-rec (- n 1)) (* 2 (f-rec (- n 2))) (* 3 (f-rec (- n 3))))))

> (f-rec 1)
1
> (f-rec 3)
4
> (f-rec 6)
59
> (f-rec 8)
335

Iterative implementation of f(n)

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

> (f 1)
1
> (f 3)
4
> (f 6)
59
> (f 8)
335

Filed under: lisp, math, SICP No Comments

26May/090

Exercise 1.10 of SICP

Exercise 1.10: The following procedure computes a mathematical function called Ackermann's function.

(define (A x y)
  (cond ((= y 0) 0)
        ((= x 0) (* 2 y))
        ((= y 1) 2)
        (else (A (- x 1)
                 (A x (- y 1))))))

What are the values of the following expressions?

> (A 1 10)
1024
> (A 2 4)
65536
> (A 3 3)
65536

Consider the following procedures, where A is the procedure defined above:
Give concise mathematical definitions for the functions computed by the procedures f, g, and h for positive integer values of n. For example, (k n) computes 5n²

(define (f n) (A 0 n))
$f(n) = 2n$
(define (g n) (A 1 n))
(g 3)
(A 1 3)
(A 0 (A 1 2))
(* 2 (A 0 (A 1 1)))
(* 2 (* 2 2)
8
$g(n) = 2^n$
(define (h n) (A 2 n))
(h 3)
(A 2 3)
(A 1 (A 2 2))
From the previous result I know that (h 3) will look like $2^{A(2,2)}$
(A 2 2)
(A 1 (A 2 1))
From the previous result I know that (h 2) will look like $2^{A(2,1)}$
$h(3) = 2^{2^{A(2,1)}} = 2^{2^2}$
$h(n) = 2^{2^{\cdot^{\cdot^{2^2}}}}$
Where the exponent repeats n times.
This operation is called tetration.
(define (k n) (* 5 n n))
$f(n) = 5n^2$

Filed under: lisp, math, SICP No Comments

20May/090

Exercise 1.8 of SICP

(define (square x) (* x x))
(define (great-enough? guess newguess)
  (< (abs (- 1 (/ guess newguess))) 0.0001))
(define (improve guess x)
  (/ (+ (/ x (square guess)) (* 2 guess))
     3))
(define (cube-root guess x)
  (if (great-enough? guess (improve guess x))
    guess
    (cube-root (improve guess x) x)))

> (cube-root 1.0 8)
2.000004911675504
> (cube-root 1.0 1e9)
1000.0162369748963

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

19May/090

Exercise 1.7 of SICP

Exercise 1.7: The good-enough? test used in computing square roots will not be very effective for finding the square roots of very small numbers. Also, in real computers, arithmetic operations are almost always performed with limited precision. This makes our test inadequate for very large numbers. Explain these statements, with examples showing how the test fails for small and large numbers. An alternative strategy for implementing good-enough? is to watch how guess changes from one iteration to the next and to stop when the change is a very small fraction of the guess. Design a square-root procedure that uses this kind of end test. Does this work better for small and large numbers?

Algorithm:
$x_{n+1} = \frac{1}{2}(x_n+\frac{S}{x_n})$

(define (average a b)
  (/ (+ a b) 2))
(define (improve guess x)
  (average guess (/ x guess)))
(define (good-enough? guess x)
  (< (abs (- x (* guess guess))) 0.0001))
(define (sqrt-iter guess x)
  (if (good-enough? guess x)
      guess
      (sqrt-iter (improve guess x)
                 x)))

> (sqrt-iter 1.0 16e-8)
.007819325057615187
> (sqrt-iter 1.0 16e64)
**Infinite Loop***

(define (great-enough? oldguess guess)
  (< (abs (- 1 (/ oldguess guess))) 0.0001))
(define (sqrt-iter2 guess x)
  (if (great-enough? guess (improve guess x))
      guess
      (sqrt-iter2 (improve guess x)
                 x)))

> (sqrt-iter2 1.0 16e-8))
4.000016244484425e-4
< (sqrt-iter2 1.0 16e8)
40000.16244495736

The reason sqrt-iter does so much worse with small numbers is because 1e-4, which is my threshold in good-enough?, happens to be much larger than 16e-8. Once "guess" start converging to the square root, let's say it becomes something close to 4e-3, (* guess guess) now becomes (* 4e-3 4e-3) which in turns becomes 16e-6. 16e-6 and 16e-8 are both much smaller than the threshold. This makes good-enough? becomes a little meaningless. In this particular case guess eventually gets down to about 0.008: (16e-8)-(.008)^2 = (16e-8) - (6.4e-5) = about (-6e-5). Since this is smaller than the threshold the test passes and "guess" get's returned.

The reason sqrt-iter does so poorly with large numbers is because of the way floating point operations work. Eventually "guess" gets somewhere close to 4e32. In the good-enough? procedure the following step gets evaluated (abs (- 16e64 (square 4.0e32))) this turns out to be about 2.33e49. Clearly that is the wrong answer. The true result should be somewhere very close to 0. Since floating point numbers have a finite amount of bits to represent a base and an exponent this means large numbers are susceptible to overflow. Which is exactly what happens here. good-enough? now returns false even though the numbers are actually within the threshold, and improve tries to make the "guess" even better but it makes no difference because the original 4.0e32 was good enough. A few places changed in the decimal's place won't make a difference to numbers this large.

great-enough? avoids both pitfalls by using ratios of numbers that are of similar magnitude. In both small and large numbers for guesses the bit patterns of new guess and old guess are similar and by taking ratios of these similar numbers the overflow problem is avoided. For numbers much smaller than the threshold their ratios are still meaningful because with respect to each other they are not so small , there is no risk of operating outside of threshold sensitivity. The disparity between each iteration of old and new guess would have to be impossibly large in order to operate outside of the threshold. See the growth of each iterative error term bellow to see why this can't theoretically ever be the case.

Error Growth for this algorithm:
$x_{n+1} = \frac{1}{2}(x_n+\frac{S}{x_n})$
$\epsilon_n = \frac{(x_n)^2}{S} - 1$
$(x_n)^2 = S(\epsilon_n+1)$
$(x_{n+1})^2 = \frac{1}{4}\left(x_n^2+2S+\frac{S^2}{x_n^2}\right)$
$\frac{(x_{n+1})^2}{S} = \frac{1}{4}\left((\epsilon_n+1)+2+\frac{1}{(\epsilon_n+1)}\right)$
$\frac{(x_{n+1})^2}{S} = \frac{1}{4}\left((\epsilon_n+1)+2+\frac{1}{(\epsilon_n+1)}\right)$
$\frac{(x_{n+1})^2}{S} = \frac{(\epsilon_n+1)^2+2(\epsilon_n+1)+1}{4(\epsilon_n+1)}$
$\frac{(x_{n+1})^2}{S} = \frac{\left((\epsilon_n+1)+1\right)^2}{4(\epsilon_n+1)}$
$\epsilon_{n+1}=\frac{(x_{n+1})^2}{S}-1 = \frac{\left((\epsilon_n+1)+1\right)^2}{4(\epsilon_n+1)}-1$
$\epsilon_{n+1}=\frac{(\epsilon_n+2)^2}{4(\epsilon_n+1)}-1$

(define (error n eps)
  (if (< (abs eps) 1e-9)
    n
    (error (+ n 1) (- (/ (square (+ eps 2))
                         (* 4 (+ eps 1)))
                      1))))

> (error 1 (- (/ 1.0 16e-8) 1))
> 16
The actual epsilon value is only: 1.6492585075411625e-11
After only 16 iterations the algorithm can return a solution with the accuracy better than 1 part in a billion. The actual value is closer to 1 part in a trillion.

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

22Mar/090

Derivation for sum of squares

My problem with mathematical induction is that a lot of the time it's simple to prove a statement is true (just do enough of the algebra calisthenics) yet just because a statement is true doesn't mean that you learned all that much about a problem. I think everyone that ever took a discrete math class had to churn out the proof for the following statement.

$1 + 2 + \cdots + n = \frac{n(n + 1)}{2}$

Certainly the wikipedia page on induction doesn't fail to include this as an example.
This is the canonical inductive proof. However to me it's not interesting. I think deriving the solution is easy enough. More interesting is deriving the solution to the sum of square or sum of cubes or an arbitrary higher power.

$1^2 + 2^2 + \cdots + n^2 = \frac{n(n+1)(2n+1)}{6}$
$1^3 + 2^3 + \cdots + n^3 = \frac{n^2(n+1)^2}{4}$

I was curious as to how these formulas came to be. It's all well and good to say, "Sure, this formula does indeed work", and it's a different thing to derive it.
For convenience I will write the identities in the following way:

$\displaystyle\sum_{1}^{n}{i} = \frac{n(n + 1)}{2}$

The crux of the derivation for higher powers lies in the fact that:
$\displaystyle\sum_{1}^{n}{(i+1)^2-i^2} = (n+1)^2-1$
to see why this is so, plug in the first few terms of the expansion of (i+1)^2 as well as the last few terms:
$(i+1)^2 = 2^2 + 3^2 + \cdots + (n-3)^2 + (n-2)^2 + (n-1)^2+ n^2 + (n+1)^2$
$i^2 = 1^2 + 2^2 + 3^2 + \cdots + (n-3)^2 + (n-2)^2 + (n-1)^2+ n^2$
Subtracting all the like terms, only two terms remain: -1 and (n+1)^2
Here's the derivation:
$\displaystyle\sum_{1}^{n}{(i+1)^2-i^2} = \displaystyle\sum_{1}^{n}{(i^2+2i+1)-i^2} = \displaystyle\sum_{1}^{n}{2i+1} = (n+1)^2-1$
$\displaystyle\sum_{1}^{n}{2i+1} = \displaystyle2\sum_{1}^{n}{i} + \displaystyle\sum_{1}^{n}{1} = \displaystyle2\sum_{1}^{n}{i} + n = (n+1)^2-1$
$\displaystyle\sum_{1}^{n}{i} = \frac{(n+1)^2-1 - n}{2}$
$\displaystyle\sum_{1}^{n}{i} = \frac{n(n+1)}{2}$

By a similar argument the sum of cubes is derived as well:
$\displaystyle\sum_{1}^{n}{(i+1)^3-i^3} = (n+1)^3-1$
$\displaystyle\sum_{1}^{n}{(i+1)^3-i^3} = \displaystyle\sum_{1}^{n}{3i^2+3i+1} = (n+1)^3-1$
$\displaystyle\sum_{1}^{n}{i^2} = \frac{(n+1)^3-1 - n - \frac{3n(n+1)}{2}}{3}$
After plenty of simplification the simplified form does indeed come out
$\displaystyle\sum_{1}^{n}{i^2} = \frac{n(n+1)(2n+1)}{6}$

The same process can be used to derive sums of cubes and any other power, as long as the sums of previous powers are known.

Filed under: math No Comments

16Jul/080

Things that surprised me

These are concepts from Math,Comp Sci, Physics that have surprised me over the years:

Recursion
My first formal introduction to recursion was through Fibonacci numbers. While the uses of these numbers are quite surprising( i.e. Euclidean algorithm's running time ). I was never really impressed with them. They seem simple enough. What really floored me about recursion, was the solution of the Tower's of Hanoi problem and subsequently the 8 Queens problem. On the surface they seemed so difficult to write a good simple solution for but where made almost trivial with recursion. Recursion still amazes me today.

Mathematical Induction
Induction ties into recursion above, but it needs it's own category in my list. Induction did not even seem fair when I first learned it. It was as if I was cheating. I could prove statements that I knew almost nothing about with minimal effort just by slogging through the algebra steps. At the end I still felt like I haven't really grasped the problem. A good example at the time was the sum of cubes.

Lambda Calculus
If there is one good way to study recursion it's lambda calculus. This field has augmented my understanding of every single item bellow and still continues to. From combinators to formal renaming schemes to deriving natural numbers(more on this later).

Lisp
I have learned as much about Computer Science from this language as I did with everything else I've encountered in this field combined: Functional Programming, OO, macros, closures, higher order functions, bottom-up programming, parsers, syntax. Each of these areas has either been completely revolutionized in my head or planted there, where before was nothing.

Dedekind Cuts
I have never imagined that it's possible to derive Real and Complex numbers out of Integers. What an amazing feat (see: Foundation of Analysis by Landau).

Derivatives and Integrals
These two things are reason I decided to also get a Math degree. I really wanted to understand why and how these work.
I remember what I thought when I saw derivatives in high school like it was yesterday. "So I take a function f(x) change an argument I am passing to it to x+h where h is so small that I can never write down a number y greater than zero that is smaller than h. Subtract f(x) and then divide this sum by another effectively zero "number" and I get rate of change with respect to x. Great!" Once that became clear, I ran into another issue.

The situation got worse with Leibniz notation. I could do things like "cancel out" infinitesimals (dx/dt)*(dt/dy)=dx/dy but I could not do this with derivatives with a higher order than 1. Makes perfect sense!

For integrals it was worse. "I take f(x) multiply it by something that makes no sense by itself and add up all the columns in the interval and BAM!... I got the function I was looking for. Great!"

Taylor's Theorem
I was delighted when I found this was possible. I always wanted to prove it. When I finally got to in Real Analysis and was able to prove it, it was one of the most satisfying experiences of my college career.

Fourier Series/Transform
This started in my Differential Equations course and made worse by my Signals and Systems EE class. In DEq I was told about Laplace transforms. In EE I was told about FTs. I had no idea how these could possibly work. Fourier series was also quite disturbing. If I add a bunch of infinitely smooth functions I can get non-smooth ones. This was a pretty big deal. These transforms are still on my list to prove.

Lagrangian and Hamiltonian Formulation of Newtonian Mechanics
I am still working on these concepts. I was never satisfied by their explanations in class. You'll have to watch the Sussman lecture to understand why.

Noether's theorem
Out of all the physics discovered this in my mind is the most amazing theorem. Even in simplified form I am in awe.

Of course the list has more entries but these are my top 10.

Filed under: lisp, math, physics No Comments

9Jul/080

Countable Rational Numbers

In my Real Analysis class in college I remember seeing the proof that rational numbers are countable via a diagonalisation argument. Proving something is true and implementing a counting algorithm that will assign any given Rational a unique positive integer are two very different things.

This link is class room explanation behind all the math without dumbing anything down. Here's the original amazing paper by virtue of simplicity and elegance. The only things needed to understand it is the most rudimentary understanding of binary trees and Euclid’s gcd algorithm.

I wrote up some scheme code to print out the "parents" of the node, which allows one to assign a unique positive integer to every/any Rational.

;Tree
(define (display-all . args) (display args))
(define (p-parents m n)
  (cond ((or (zero? m) (zero? n)))
        ((= m n) (display-all 1 "\n") 1)
        ((< m n) (display-all m " " n "\n")(p-parents m (- n m)))         ((&gt; m n) (display-all m " " n "\n")(p-parents (- m n) n))))
 
;Print the parent nodes of 49/19
(p-parents 49 19)

Filed under: lisp, math No Comments

8Jul/080

World Magnetic Model

In about a dozen or so math courses I remember having to prove (always by induction) that the sum of N numbers (1+2+3+...+N) is equal to n(n+1)/2 and always thinking, this is pointless. I stand corrected.

At work I was tasked to calculate magnetic declination using java instead of the C implementation NOAA used. The program is pretty much a function: MagDec:(double,double,double)->double (elevation,latitude,longitude)->magnetic_dec. A key component used to calculate the declination is a file with spherical harmonics gathered by satellites. It looks something like this:

          G(m,n)   H(m,n)       G_dot(m,n)  H_dot(m,n)
1  0  -29556.8       0.0        8.0        0.0
1  1   -1671.7    5079.8       10.6      -20.9
2  0   -2340.6       0.0      -15.1        0.0
2  1    3046.9   -2594.7       -7.8      -23.2
2  2    1657.0    -516.7       -0.8      -14.6
3  0    1335.4       0.0        0.4        0.0
3  1   -2305.1    -199.9       -2.6        5.0
3  2    1246.7     269.3       -1.2       -7.0
3  3     674.0    -524.2       -6.5       -0.6
4  0     919.8       0.0       -2.5        0.0
4  1     798.1     281.5        2.8        2.2
4  2     211.3    -226.0       -7.0        1.6
....
....
12 10      -0.1      -0.9        0.0        0.0
12 11      -0.3      -0.4        0.0        0.0
12 12      -0.1       0.8        0.0        0.0

My goal was to create S = G(m,n) - G_dot(m,n) and C = H(m,n) - H_dot(m,n) for each value of m and n.
Since Java IO isn't the most convenient thing in the world, instead of sscanf() I was forced to use java.util.Scanner, which IMO is not nearly as convenient for this type of data organization.
I realized that if I have a flat list it looks something like this:

Col. A   Col. B
G(1,0)  -29556.8
        8.0
G(1,1)  -1671.7
        10.6
G(2,0)  -2340.6
       -15.1
G(2,1)   3046.9
       -7.8
G(2,2)   1657.0
       -0.8

Written in a different way:
G(1) G(2) G(3) ... G(12)
2 3 4 ... 13

which means there are 13(13+1)/2 - 1 coefficients for G. If I am presented with a single column of numbers like Col B. and I want to pick G(2,1) all I need to do is is realize that there are 3(3+1)/2 coefficients up to G(2,2) which means G(2,1) is 3(3+1)/2-1. Of course it would've been a lot easier if I could do the straightforward thing easily like:

WMM = open("WMM.COF")
coeffs = []
try:
   for line in WMM:
       nline = []
       for each in line.split():
           nline.append(float(each))
       coeffs.append(nline)
finally:
   WMM.close()

Filed under: math, physics No Comments

7Jul/080

School Algebra Story Problems

This was taken from Marvin Minsky's Why programming is a good medium for expressing poorly understood and sloppily-formulated ideas.

Mary is twice as old as Ann was when Mary was as old as Ann is now. If Mary is 24 years old, how old is Ann?

My school problems weren't worded in quite this confusing a way.

Filed under: ai, math No Comments

« Newer Entries

Dan's Thoughts Thinking somewhat carefully

Exercise 1.12 of SICP

Exercise 1.11 of SICP

Exercise 1.10 of SICP

Exercise 1.8 of SICP

Exercise 1.7 of SICP

Derivation for sum of squares

Things that surprised me

Countable Rational Numbers

World Magnetic Model

School Algebra Story Problems

Pages

Categories

Blogroll

Archive

Meta