Interesting Miscellany

Here are some of the links I've found interesting enough to tweet or retweet in that past several weeks.

Programming
A beautiful obfuscated Pi calculator
Programming language trends of 2009
Tweetable game of life in MATLAB
Factorization of RSA768
Calculating Sines


Science
Leonardo da Vinci's resume
Richard Dawkins: Growing up in the universe
Richard Hamming: You and your research


Math
A quick little puzzle (in probability)
Euler Book Prize (A list to watch in the future. Since the site isn't updated for 2010 yet, the winner was Euler's Gem by David S. Richeson.)
A cryptic math advertisement


Miscellaneous
Succed Blog
RockYou.com's 20 most common passwords
50 Books Every Geek Should Read

SICP Exercise 1.21: Smallest Divisor

From SICP section 1.2.6 Example: Testing for Primality

Exercise 1.21 asks us to use the smallest-divisor procedure to find the smallest divisor of each of the following numbers: 199, 1999, 19999.

The smallest-divisor procedure was given earlier in the section.

(define (smallest-divisor n)
   (find-divisor n 2))

(define (find-divisor n test-divisor)
   (cond ((> (square test-divisor) n) n)
         ((divides? test-divisor n) test-divisor)
         (else (find-divisor n (+ test-divisor 1)))))

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

This procedure makes use of the square procedure that was defined earlier. You'll need to include it if you want to run the code in an interpreter.

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

The smallest-divisor procedure finds the smallest integral divisor (greater than 1) of a given number n in a straightforward way, by testing to see if n is divisible by each integer from 2 to √n. Since the number of steps is between 1 and √n, the order of growth of the algorithm is O(√n).

We can find the solutions to the exercise by simply running the code in an interpreter.

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

As you can see, the first two numbers are prime because they are their own smallest divisor, but the third is divisible by 7.


Related:

For links to all of the SICP lecture notes and exercises that I've done so far, see The SICP Challenge.

SICP Exercise 1.20: GCD

From SICP section 1.2.5 Greatest Common Divisors

Exercise 1.20 asks us to consider the following iterative procedure that uses Euclid's Algorithm to comput the GCD:

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

We're asked to illustrate the process generated by this procedure when using normal-order and applicative-order evaluation rules to evaluate (gcd 206 40). How many remainder operations are actually performed using each evaluation model?

We learned the rules for both applicative-order and normal-order evaluation back in SICP section 1.1.5 The Substitution Model for Procedure Application. We learned that under normal-order evaluation rules, the interpreter fully expands a procedure before reducing it. Operand expressions are substituted for parameters until an expression involving only primitive operators is reached. Under applicative-order rules, arguments are evaluated before operators are applied. This can often avoid multiple unnecessary evaluations of the same expression, which is one of the reasons why Lisp uses applicative-order evaluation.

In section 1.1.6 Conditional Expressions and Predicates (see exercise 1.5) we're told to assume that the evaluation rule for the special form if is the same whether we use normal or applicative order evaluation. "The predicate expression is evaluated first, and the result determines whether to evaluate the consequent or the alternative expression." We'll continue using this assumption until we're told otherwise.

Now let's put all of this together to illustrate both processes for (gcd 206 40).


Normal-order evaluation - fully expand and then reduce.

(gcd 206 40)

(if (= 40 0)
    206
    (gcd 40 (remainder 206 40)))

(gcd 40 (remainder 206 40))

(if (= (remainder 206 40) 0)
    40
    (gcd (remainder 206 40) (remainder 40 (remainder 206 40))))

(if (= 6 0) ;; 1 remainder evaluated
    40
    (gcd (remainder 206 40) (remainder 40 (remainder 206 40))))

(gcd (remainder 206 40) (remainder 40 (remainder 206 40)))

(if (= (remainder 40 (remainder 206 40)) 0)
    (remainder 206 40)
    (gcd (remainder 40 (remainder 206 40))
         (remainder (remainder 206 40)
                    (remainder 40 (remainder 206 40)))))

At this point both of the nested remainder operations that were substituted in for b in the (= b 0) procedure need to be evaluated. I'll combine these and show them as one step.

(if (= 4 0) ;; 2 remainders evaluated
    (remainder 206 40)
    (gcd (remainder 40 (remainder 206 40))
         (remainder (remainder 206 40)
                    (remainder 40 (remainder 206 40)))))

(gcd (remainder 40 (remainder 206 40))
     (remainder (remainder 206 40)
                (remainder 40 (remainder 206 40))))

(if (= (remainder (remainder 206 40)
                  (remainder 40 (remainder 206 40)))
       0)
       (remainder 40 (remainder 206 40))
       (gcd (remainder (remainder 206 40)
                       (remainder 40 (remainder 206 40)))
            (remainder (remainder 40 (remainder 206 40))
                       (remainder (remainder 206 40)
                                  (remainder 40 (remainder 206 40))))))

Here again, multiple nested remainder operations that were substituted in for b in the (= b 0) procedure need to be evaluated. I'll combine all four and show them as one step.

(if (= 2 0) ;; 4 remainders evaluated
    (remainder 40 (remainder 206 40))
    (gcd (remainder (remainder 206 40)
                    (remainder 40 (remainder 206 40)))
         (remainder (remainder 40 (remainder 206 40))
                    (remainder (remainder 206 40)
                               (remainder 40 (remainder 206 40))))))

(gcd (remainder (remainder 206 40)
                (remainder 40 (remainder 206 40)))
     (remainder (remainder 40 (remainder 206 40))
                (remainder (remainder 206 40)
                           (remainder 40 (remainder 206 40)))))

(if (= (remainder (remainder 40 (remainder 206 40))
                  (remainder (remainder 206 40)
                             (remainder 40 (remainder 206 40)))) 0)
    (remainder (remainder 206 40)
               (remainder 40 (remainder 206 40)))
    (gcd (remainder (remainder 40 (remainder 206 40))
                    (remainder (remainder 206 40)
                               (remainder 40 (remainder 206 40))))
         (remainder (remainder (remainder 206 40)
                               (remainder 40 (remainder 206 40)))
                    (remainder (remainder 40 (remainder 206 40))
                               (remainder (remainder 206 40)
                                          (remainder 40
                                                     (remainder 206 40)))))))

There are seven nested remainder operations to evaluate at this point.

(if (= 0 0) ;; 7 remainders evaluated
    (remainder (remainder 206 40)
               (remainder 40 (remainder 206 40)))
    (gcd (remainder (remainder 40 (remainder 206 40))
                    (remainder (remainder 206 40)
                               (remainder 40 (remainder 206 40))))
         (remainder (remainder (remainder 206 40)
                               (remainder 40 (remainder 206 40)))
                    (remainder (remainder 40 (remainder 206 40))
                               (remainder (remainder 206 40)
                                          (remainder 40
                                                     (remainder 206
                                                                40)))))))

(remainder (remainder 206 40)
           (remainder 40 (remainder 206 40)))

Finally, we can evaluate these four remainder operations to come up with the answer (GCD) and the total.

2 ;; 18 total remainders evaluated

Applicative-order evaluation - evaluate arguments and then apply operands. This one is a little bit easier to follow. With each substitution, we just evaluate the operands that we need before applying the operators at each step.

(gcd 206 40)

(if (= 40 0)
    206
    (gcd 40 (remainder 206 40))))

(gcd 40 (remainder 206 40))

(gcd 40 6) ;; 1st remainder evaluated

(if (= 6 0)
    40
    (gcd 6 (remainder 40 6)))

(gcd 6 (remainder 40 6))

(gcd 6 4) ;; 2nd remainder evaluated

(if (= 4 0)
    6
    (gcd 4 (remainder 6 4)))

(gcd 4 (remainder 6 4))

(gcd 4 2) ;; 3rd remainder evaluated

(if (= 2 0)
    4
    (gcd 2 (remainder 4 2)))

(gcd 2 (remainder 4 2))

(gcd 2 0) ;; 4th remainder evaluated

(if (= 0 0)
    2
    (gcd 0 (remainder 2 0)))

2

As you can see, using applicative-order evaluation, remainder was only evaluated four times, compared to 18 times using normal-order evaluation.


Related:

For links to all of the SICP lecture notes and exercises that I've done so far, see The SICP Challenge.

Unsolved: Odd Perfect Numbers

A perfect number is a number that is the sum of its proper divisors (i.e., including 1 but excluding itself).

For example, the smallest perfect number is 6.

1 + 2 + 3 = 6

The next two perfect numbers are

1 + 2 + 4 + 7 + 14 = 28
1 + 2 + 4 + 8 + 16 + 31 + 62 + 124 + 248 = 496

See sequence A000396 in OEIS for a list of the first ten perfect numbers.

On a related note, a number that is smaller than the sum of its proper divisors is an abundant number, while a number greater than the sum of its proper divisors is called deficient. The smallest abundant number is 12, whose proper divisors sum to 16.

All of the perfect numbers that have been discovered so far are even.

The problem is to find an odd perfect number, or prove that no odd perfect numbers exist.

If there is an odd perfect number, we already know a lot about it. It must be at least 300 digits long, with at least 75 prime factors, at least 9 of them distinct. Its largest prime factor must be greater than 100,000,000. Its second largest prime factor is greater than 10,000, and its third largest is greater than 100.

SICP Exercise 1.19: Computing Fibonacci numbers

From SICP section 1.2.4 Exponentiation

Exercise 1.19 asks us to complete a procedure for computing Fibonacci numbers in a logarithmic number of steps. The following code is given:

(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
                    <??>      ; compute p'
                    <??>      ; compute q'
                    (/ count 2)))
         (else (fib-iter (+ (* b q) (* a q) (* a p))
                         (+ (* b p) (* a q))
                         p
                         q
                         (- count 1)))))

We're reminded of the transformation of the state variables a and b in the original fib-iter procedure from section 1.2.2: a ← a + b and b ← a. If these state changes are labeled transformation T, then applying T repeatedly for n iterations starting with a = 1 and b = 0 produces the pair a = Fib(n + 1) and b = Fib(n). So the Fibonacci numbers are produced by the nth power of the transformation T, or Tⁿ, starting with the pair (1, 0).

We are then asked to consider the family of transformations T_pq which transforms the pair (a, b) according to the following rules:

a ← bq + aq + ap
b ← bp + aq

We can verify by quick substitution that the original transformation T is just a special case of T_pq, where p = 0 and q = 1.

a ← b(1) + a(1) + a(0)
a ← b + a

b ← b(0) + a(1)
b ← a

We are asked to show that if we apply T_pq twice, the effect is the same as using a single transformation T_p'q' of the same form, and compute p' and q' in terms of p and q. This will give us an explicit way to square these transformations, which we can use to compute Tⁿ using successive squaring, just like the fast-expt procedure from exercise 1.16.

We can apply T_pq twice by defining new variables and using substitution. Let's define a₁ and b₁ as the results of applying transformation T_pq once

a₁ = bq + aq + ap
b₁ = bp + aq

The next step is to define a₂ and b₂ and apply the tranformation a second time, this time using a₁ and b₁ in place of a and b.

a₂ = b₁q + a₁q + a₁p
b₂ = b₁p + a₁q

Now that we have a system of equations defined, we can use substitution on our way to simplifying.

a₂ = (bp + aq)q + (bq + aq + ap)q + (bq + aq + ap)p
b₂ = (bp + aq)p + (bq + aq + ap)q

The second equation is shorter, so it should be easier to manipulate. Remember, we're trying to find p' and q', so we need to rewrite the equation to fit the form

b₂ = bp' + aq'

where p' and q' can be computed in terms of q and p.

b₂ = (bp + aq)p + (bq + aq + ap)q
= (bpp + apq) + (bqq + aqq + apq)
= bpp + apq + bqq + aqq + apq
= (bpp + bqq) + (2apq + aqq)
= b(pp + qq) + a(2qp + qq)

From this we can see that p' and q' can be computed using the following equations:

p' = p² + q²
q' = 2pq + q²

Manipulating the equation for a₂ in the same way, we can verify those results. This time we're trying to fit the form

a₂ = bq' + aq' + ap'

The groupings required for this manipulation are made even easier by the fact that we now already know the formulas for p' and q'.

a₂ = (bp + aq)q + (bq + aq + ap)q + (bq + aq + ap)p
= (bpq + aqq) + (bqq + aqq + apq) + (bpq + apq + app)
= bpq + aqq + bqq + aqq + apq + bpq + apq + app
= (bpq + bpq + bqq) + (apq + apq + aqq) + (app + aqq)
= b(pq + pq + qq) + a(pq + pq + qq) + a(pp + qq)
= b(2pq + qq) + a(2pq + qq) + a(pp + qq)

Now that we've verified the formulas for p' and q' in terms of p and q, we can use them to complete the procedure we were given.

(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
                    (+ (* p p) (* q q))     ; compute p'
                    (+ (* 2 p q) (* q q))   ; compute q'
                    (/ count 2)))
         (else (fib-iter (+ (* b q) (* a q) (* a p))
                         (+ (* b p) (* a q))
                         p
                         q
                         (- count 1)))))

You can test several known values in an interpreter to verify that this actually works.

> (fib 0)
0
> (fib 1)
1
> (fib 2)
1
> (fib 3)
2
> (fib 5)
5
> (fib 10)
55
> (fib 20)
6765

Related:

For links to all of the SICP lecture notes and exercises that I've done so far, see The SICP Challenge.

SICP Exercise 1.18: Iterative multiplication

From SICP section 1.2.4 Exponentiation

Exercise 1.18 asks us to use the results from the two previous exercises (1.16 and 1.17) to devise a procedure that generates an iterative process for multiplying two integers in terms of adding, doubling and halving, and which uses a logarithmic number of steps.

The key idea from exercise 1.16 was to keep an additional state variable a, and manipulate the state transition so that the product (a * bⁿ) remained unchanged. We can modify our solution to exercise 1.17 to use the same idea. Since we're now dealing with multiplication instead of exponentiation, the invariant quantity will be (a * b + p) where a and b are the two operands, and p will be our state variable. Initially p will be 0, but by the end of the iterative process it will hold the product (a * b).

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

(define (double x)
   (+ x x))

(define (halve x)
   (/ x 2))

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

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

Related:

For links to all of the SICP lecture notes and exercises that I've done so far, see The SICP Challenge.

SICP Exercise 1.17: Multiplication by repeated addition

From SICP section 1.2.4 Exponentiation

Exercise 1.17 asks us to design a multiplication procedure analogous to fast-expt that uses a logarithmic number of steps. We're told that in addition to adding, we can use procedures for doubling or halving an integer argument.

(define (double x)
   (+ x x))

(define (halve x)
   (/ x 2))

We'll also use the same even? procedure we've seen before:

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

The following multiplication procedure that takes a linear number of steps is given:

(define (* a b)
   (if (= b 0)
       0
       (+ a (* a (- b 1)))))

This procedure defines multiplication in terms of the algorithm we all learned in grade school: add a to itself b times.

The exercise also makes reference to the following fast-expt procedure, which we can use as a guide:

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

This procedure makes use of the following observations:

bⁿ = (b^n/2)² if n is even
bⁿ = b * b^n-1 if n is odd

We can define similar rules for integer multiplication.

a * b = 2 * (a * b/2) if be is even
a * b = a + a * (b - 1) if b is odd

The first rule tells us when to use our procedures double and halve. The second rule just lets us take an odd b and make it even in preparation for the next iteration. Translating these rule into a Scheme procedure give us:

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

The procedure takes a logarithmic number of steps. If you double the size of the input parameter b, the procedure takes only one more step to compute the product.


Related:

For links to all of the SICP lecture notes and exercises that I've done so far, see The SICP Challenge.

Guest article on Programming Praxis

I was excited to be asked to be the guest author for an article on one of my favorite programming blogs, Programming Praxis. If you're not familiar with the site, Programming Praxis provides weekly practice exercieses to "sharpen your saw" on. It was inspired in part by Project Euler, but the problems on Programming Praxis focus more on developing programming skills and usually require less advanced math than those on Project Euler.

To follow up earlier exercises on Calculating Pi and Logarithms, my article is on Calculating Sines using two different methods, the Taylor series and the Triple-angle formula. The sample solutions are written in Scheme, but readers are encouraged to submit solutions in any language they choose (see the HOWTO on posting source code to post your own solution). Since I'm still a relative newbie to Scheme (I'm currently still in chapter one of SICP), I'll be interested to see what alternate solutions you can come up with.

The article is also currently featured on Planet Scheme.

SICP Exercise 1.16: Fast Exponentiation

From SICP section 1.2.4 Exponentiation

Exercise 1.16 asks us to design an exponentiation procedure that evolves an iterative process using successive squaring and uses a logarithmic number of steps. A hint tells us to use the observation that (b^n/2)² = (b²)^n/2 and to keep, along with the exponent n and base b, an additional state variable a, and to define the state transformation in such a way that the product a * bn is unchanged from state to state. The value of a should start at 1, and should contain the final answer by the end of the process.

Earlier in section 1.2.4, we were given the procedure:

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

which tests whether an integer is even. This will come in handy again for this exercise, since we'll need to adjust a by a different amount when n is even than when n is odd.

We'll also use the following simple procedure for squaring a number:

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

With those building blocks in place, we can define our iterative procedure.

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

We call the procedure expt-iter with the arguments base b, exponent n, and the state variable a. The first thing we do is check to see if the exponent is 0. If so, we've reached the end of the procedure and just return the value of a. Next we check to see if n is even. If it is, we use the observation that (b^n/2)² = (b²)^n/2 to transform state information by squaring b and dividing n by 2. Finally, if n is odd we transform state information by reducing n by 1 and multiplying a by the base b.

The only thing that's left to do is define a procedure that calls expt-iter with the correct initial values.

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

Here are a few sample runs of the procedure:

> (fast-expt 2 2)
4
> (fast-expt 2 3)
8
> (fast-expt 2 10)
1024
> (fast-expt 3 3)
27
> (fast-expt 5 5)
3125

You can try running the procedure with large values of n to verify that it really is fast.


Related:

For links to all of the SICP lecture notes and exercises that I've done so far, see The SICP Challenge.

Math visualization: Multiplying Two Binomials

In my earlier post Math visualization: (x + 1)², I showed how to draw a picture that proved the equality

(x + 1)² = x² + 2x + 1

An alert reader commented that (x + 1)² is a specific instance of the binomial squared

(x + a)² = x² + 2ax + a²

Here's a visualization of that equality:



It turns out that the illustration can be generalized a little bit further to cover multiplying two binomials that form a rectangle (not just a perfect square).

(x + a)(x + b) = x² + ax + bx + ab

If you want to see how this idea can be extended into the third dimension, check out Montessori World's page on the Binomial Cube.

Unsolved: Twin Primes Conjecture

Prime numbers, as most school children can tell you, are those numbers that are evenly divisible by only themselves and 1. The first few primes are

2, 3, 5, 7, 11, 13, 17...

We've known that there are infinitely many primes since c. 300 BC when it was proven by Euclid of Alexandria. Euclid used a very simple and elegant proof by contradiction.

Euclid's proof of the infinitude of primes

First, assume that there are a finite number of primes, p1, p2, p3, ..., pn.

Now let

Q = (p1 * p2 * ... * pn) + 1

That is, Q is equal to all of the primes multiplied together plus one.

By the Fundamental Theorem of Arithmetic, Q is either prime or it can be written as the product of two or more primes. However, none of the primes in our list evenly divides Q. If any prime in our list did evenly divide Q, then that same prime would also evenly divide 1, since

Q - (p1 * p2 * ... * pn) = 1

This contradicts the assumption that we had listed all the primes. So no matter how many primes we start with in our list, there must always be more primes. (Note that this proof does not claim that Q itself is prime, just that there must be some prime not in the initial list.)


Twin primes

Twin primes are those pairs of numers that have a difference of two, and that are both prime.

3, 5
5, 7
11, 13
17, 19
...

The Twin prime conjecture, which dates back to the 18th century, simply states that there are infinitely many twin primes.

Given the simplicity of Euclid's proof of the infinitude of primes, it's tempting to hope for an equally simple proof to the Twin prime conjecture. Needless to say, such a proof has not been found.


Other facts about twin primes:

Other than (3, 5), all twin primes have the form (6n - 1, 6n + 1).

In 1919 Viggo Brun showed that the sum of the reciprocals of the twin primes converges to a definite number, now known as Brun's constant (approximately 1.902160578).

In 1994, while in the process of estimating Brun's constant by calculating the twin primes up to 10¹⁴, Thomas Nicely discovered the infamous Pentium bug.

The largest known twin primes (as of January 2010) are a pair of 100,355 digit primes with the values 65516468355 * 2³³³³³³ ± 1.