A Brief History of Pi

The date March 14th, or 3/14, is Pi Day in the U.S., so I thought I'd commemorate by sharing a brief history of one of the most important and frequently occurring numbers in mathematics.

Pi (usually denoted by the Greek letter π) is the ratio of the circumference to the diameter of a circle.

π = C / d

This ratio is constant, regardless of the size of the circle. π can also be defined as the ratio of the area to the square of the radius of a circle.

π = A / r²

Time line

c. 2000 BC - Babylonians put the ratio of the circumference of a circle to its diameter at around 3. Later estimates (around 1900 B.C.) put the value at 25/8 (3.125). At around the same time in Egypt, the approximation 256/81 (about 3.16) was used, and in India, 339/108 (about 3.139).

c. 225 BC - Archimedes of Syracuse showed that the value of π was between 3¹⁰/₇₁ and 3¹/₇. He used a method of inscribing a circle inside regular polygons, then estimating π based on their perimeters.

5th century BC - Chinese mathematician and astronomer Zu Chongzhi calculated the value of the ratio of the circumference of a circle to its diameter to be 355/113, which is accurate to seven digits. This is the best approximation to the true value of π of any fraction with a denominator less than 30,000.

c. 1400 - Madhava of Sangamagrama used an infinite series method to estimate π to be 3.14159265359, which is accurate to eleven decimal places.

1424 - Jamshīd al-Kāshī gave an estimate π that is correct to 16 digits.

c. 1600 - Ludolph van Ceulen computed the first 35 digits of π, and had the digits engraved on his tombstone in 1610.

1647 - William Oughtred used π.δ as an abbreviation for "periphery-diameter" in his Clavis mathematicae (The Key to Mathematics).

1706 - William Jones was the first to use the symbol π on its own in his Synopsis Palmariorum Matheseos.

1737 - Leonhard Euler adopted the use of the symbol π, cementing its popularity.

1761 - Johann Lambert proved that π is irrational (π cannot be expressed as the ratio of two integers).

1789 - Jurij Vega calculated 140 digits of π, but only the first 126 were correct.

1881 - Ferdinand Lindemann showed that π is a transcendental number (there is no polynomial with rational coefficients of which π is a root).

1873 - William Shanks computed 707 digits of π, but only the first 527 were correct.

1949 - John von Neumann used ENIAC to compute 2037 digits of π, using about 70 hours of computing time.

1973 - The 1,000,000th digit of π was computed using a CDC 7600 supercomputer.

1999 - Yasumasa Kanada lead a group at Tokyo University that computed π to more than 206 billion places.

2002 - The Tokyo University group beat their own record by computing π to more than 1.24 trillion digits using 600 hours of computing time on a Hitachi SR8000 supercomputer.

2009 - Daisuke Takahashi at the University of Tsukuba in Japan, calculated nearly 2.6 trillion digits of π in 29 hours on a T2K Open Supercomputer. In the same year, Fabrice Bellard computed nearly 2.7 trillion digits of π in a total of 131 days on a 2.93 GHz Core i7 (desktop class) CPU.


Piphilology

The first 50 digits of π are:

3.14159265358979323846264338327950288419716939937510

Piphilology is the art of composing mnemonic phrases for remembering the digits of π. For example, the first nine digits are encoded in the number of letters in each word of the phrase

How I wish I could recollect pi easily today!

Sir James Jeans embedded the first 15 digits in the following phrase (sometimes attributed to Isaac Asimov, there are several variations):

How I want a drink, alcoholic of course, after the heavy lectures involving quantum mechanics!

The first 740 digits of π are encoded in Michael Keith's 1995 variation of Edgar Allen Poe's The Raven. Here's the first section of 42 words:

Poe, E. Near A Raven

Midnights so dreary, tired and weary,
Silently pondering volumes extolling all by-now obsolete lore.
During my rather long nap - the weirdest tap!
An ominous vibrating sound disturbing my chamber's antedoor.
"This," I whispered quietly, "I ignore."
...

Keith followed this poem (or more appropriately, piem) in 1996 with Cadaeic Cadenza, a short story which encodes the first 3,835 digits of π, and in early 2010 he published Not a Wake, a novelette that recounts "A dream embodying π's digits fully for 10,000 decimals."


Additional References:

History of Pi by Petr Beckmann

The Universal Book of Mathematics by David Darling

For a much more detailed record of the computation of the digits of π, see the Chronology of computation of π.

SICP Exercise 1.28: The Miller-Rabin test

From SICP section 1.2.6 Example: Testing for Primality

(Note: The explanation of the Miller-Rabin test given in SICP is slightly different from every other source I've read. I don't know if this was done purposefully by the authors for the sake of simplicity, or if the algorithm has simply evolved since SICP's publication. For the sake of consistency, I'm solving the exercise using the explanation given in the text. If anyone is interested, Programming Praxis has a short explanation of the algorithm and an implementation in Scheme. Naturally, you can also read about it on Wikipedia.)

Exercise 1.28 introduces a variation on the Fermat test that isn't fooled by Carmichael numbers called the Miller-Rabin test. We start with an alternate form of Fermat's Little Theorem. In exercise 1.24 we saw that if n is prime and a is any positve integer < n,

aⁿ ≡ a (mod n)

If we divide both sides of the congruency by a, we get

a^n-1 ≡ 1 (mod n)

When we check this congruency using the expmod procedure, at the squaring step we need to also check to see if we've discovered a "non-trivial square root of 1 modulo n," i.e., a number not equal to 1 or (n - 1) whose square is equal to 1 modulo n. It's possible to prove that if a non-trivial square root of 1 exists, then n is not prime. It's also possible to prove that if n is composite, then "for at least half the numbers a < n, computing a^n-1 in this way will reveal a nontrivial square root of 1 modulo n." (I'm not sure if this is just a very conservative estimate or if it was the best information available at the time that SICP was written. Every other source I've read claims that at least 3/4 of the bases you choose will reveal a nontrivial square root of 1 modulo n.)

We're asked to modify the expmod procedure to check for non-trivial square roots of 1 modulo n, and to use it to implement the Miller-Rabin test with a procedure analogous to the fermat-test that we wrote before. We can check the new procedure by testing various known primes and non-primes.

The first thing we need to do is modify the square procedure so that it checks for a non-trivial square root of 1 modulo n. We'll write a new procedure called square-check that checks for two conditions:

1. a number not equal to 1 or (n - 1)
2. whose square is equal to 1 modulo n

If both of these conditions are met, we'll signal that we've found a non-trivial square root by returning 0. If these conditions are not met, we'll perform the remainder and squaring steps that expmod did before. Here's the complete code listing:

(require (lib "27.ss" "srfi"))

(define (square-check x m)
  (if (and (not (or (= x 1) (= x (- m 1))))
           (= (remainder (* x x) m) 1))
      0
      (remainder (* x x) m)))

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

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

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

As with the fermat-test procedure from exercise 1.27, prime numbers should always pass. Composite numbers should have a high likelihood of failing, even the Carmichael numbers (so if you repeatedly test these values enough times, eventually you will see some false positives). You can test several primes and composites in a Scheme interpreter. Here's the output I got when I tested the first six Carmichael numbers.

> (miller-rabin-test 561)
#f
> (miller-rabin-test 1105)
#f
> (miller-rabin-test 1729)
#f
> (miller-rabin-test 2465)
#f
> (miller-rabin-test 2821)
#f
> (miller-rabin-test 6601)
#f

Since the Miller-Rabin test is a refinement of the Fermat test that we saw earlier, it's natural to assume that you can substitute our new procedure for the old one. This is a correct assumption. To upgrade the fast-prime? procedure we saw in exercise 1.24, just change the call to fermat-test into a call to miller-rabin-test as follows:

(define (fast-prime? n times)
  (cond ((= times 0) true)
        ((miller-rabin-test n) (fast-prime? n (- times 1)))
        (else false)))

Related:

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

"Goes to" Considered Harmful

A few months ago the following question was asked (half jokingly) on Stack Overflow,

What is the name of this operator: “-->”?

The question includes the following C++ code:

#include <stdio.h>
int main()
{
    int x = 10;
    while( x --> 0 ) // x goes to 0
    {
        printf("%d ", x);
    }
}

The output of this program is to simply count down from 9 to 0, printing each value. As many others pointed out, "-->" isn't a single operator at all, but two operators smashed together. The condition of the while loop above can be rewritten with proper spacing as:

while( x-- > 0 ) // x-- is greater than  0
{
    ...
}

Joshua Bloch commented on the goes-to operator on Twitter yesterday.



He links to the results of a Google Code search for instances of "-->" in the wild. The results are specific to the C programming language, but "-->" turns up in similar searches for C++ and Java as well (and probably any other language that allows it).

What's the harm?

So what's so harmful about this bit of cleverness? The main problem is that it's too clever. Any time you use standard operators in a non-standard way (even if the language specification doesn't strictly forbid it) you're negatively impacting the readability of your code. In this specific case, the "-->" operator isn't documented anywhere. It's just not reasonable to expect other programmers to know what it means.

On top of that, this situation isn't best handled by a while loop to begin with. If you know ahead of time how many iterations your loop needs to make, use a for loop. Kernighan and Richie covered this a very long time ago in The C Programming Language. It's still true today.

The for is preferable when there is a simple initialization and increment, since it keeps the loop control statements close together and visible at the top of the loop.

This, of course, is also true for a simple decrementing loop. The meaning of the following code should be obvious to any junior programmer.

int i;
for( i = 9; i >= 0; i-- )
{
    printf("%d ", i);
}

Writing code that's "too clever" for other programmers to understand immediately just obfuscates your code needlessly. You'll seem much cleverer to your colleagues if you write your code in a clear and readable style to begin with. Besides that, using the common idioms of your language is probably the quickest way to reduce your WTFs/minute score at your next code review.

SICP Exercise 1.27: Carmichael numbers

From SICP section 1.2.6 Example: Testing for Primality

Exercise 1.27 asks us to demonstrate that the first six Carmicheal numbers (561, 1105, 1729, 2465, 2821, and 6601) really do fool the Fermat test. We're to write a procedure that takes an integer n and tests whether aⁿ is congruent to a modulo n for every value less than n, then test the procedure on each of the Carmichael numbers listed.

The fermat-test procedure that we first saw in exercise 1.24 is a good starting point.

(define (fermat-test n)
   (define (try-it a)
     (= (expmod a n n) a))
   (try-it (+ 1 (random-integer (- n 1)))))

Let' start by modifying this procedure so that we pass it the value to test, instead of generating a random value.

(define (fermat-test n a)
   (= (expmod a n n) a))

Now we just need to define a new procedure that will call this one for every value less than the one we want to test. If fermat-test passes for all values tested we should return true, but if it fails for any value we should return false. Here is the complete code listing for the exercise.

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

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

(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 (fermat-test n a)
   (= (expmod a n n) a))

(define (fermat-full n)
   (define (iter a)
     (cond ((= a 1) #t)
           ((not (fermat-test n a)) #f)
           (else (iter (- a 1)))))
   (iter (- n 1)))

Prime numbers will legitimately pass this test, and Carmichael numbers will fool it. All other composites should fail. You can test out a few primes and composites in a Scheme interpreter to make sure it works as described. Here's the output for the first six Carmichael numbers.

> (fermat-full 561)
#t
> (fermat-full 1105)
#t
> (fermat-full 1729)
#t
> (fermat-full 2465)
#t
> (fermat-full 2821)
#t
> (fermat-full 6601)
#t

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.26: Explicit multiplication vs. squaring

From SICP section 1.2.6 Example: Testing for Primality

Exercise 1.26 asks us to consider yet another variation on the fast-prime? procedure from exercise 1.24. This time the expmod procedure has been modified to use an explicit multiplication instead of calling square:

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

This change causes fast-prime? to run more slowly than the original prime? test by transforming the O(log n) process into a O(n) process. We're asked to explain this transformation.

A cursory inspection of the code makes it seem like the explicit multiplication will cause twice as many calls to expmod because each input argument to * will be evaluated separately, instead of only once when expmod is passed as the single argument to square. This isn't enough to account for the reported slow down.

Let's take a closer look at the process generated by each version of expmod. (Since the only difference between the two versions of expmod is in the branch where exp is even, I'm using powers of 2 for that argument to better illustrate the difference in the resulting processes. I should point out that this is a worst case scenario.)

Here is what the process generated by the original expmod procedure using square might look like:

(expmod 2 8 2)
   (expmod 2 4 2)
       (expmod 2 2 2)
           (expmod 2 1 2)
               (expmod 2 0 2)
               1
           0
       0
   0
0

This is a pretty straightforward linear recursive process. Now let's look at the process for expmod when explicit multiplication is used. (This is only part of the process diagram, but it's enough to see what's going on.)


Right away, it's easy to see what went wrong here. By replacing the call to square with explicit multiplication, we've transformed expmod from a linear recursive process (like the factorial example) into a tree recursion (like the original Fibonacci example). This causes the number of recursive calls to expmod to grow exponentially instead of simply doubling. What was once a O(log n) process is now O(log 2ⁿ), which simplifies to O(n).


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.25: A closer look at expmod

From SICP section 1.2.6 Example: Testing for Primality

Exercise 1.25 asks us to take a closer look at the expmod procedure that we used in exercise 1.24 to compute the exponential of a number modulo another number:

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

We're asked to consider a more straightforward implementation of the same function. This version makes use of the fast-expt procedure from section 1.2.4, which I've also included here.

(define (expmod base exp m)
   (remainder (fast-expt base exp) m))

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

Will this version perform just as well in the fast prime tester from exercise 1.24? The easiest way to answer that question is to drop it in and test it out. Starting out with relatively small values, you can see that the new procedure doesn't perform very well at all.

> (timed-prime-test 1999)

1999 *** 138.0
> (timed-prime-test 10007)

10007 *** 681.0
> (timed-prime-test 100003)

100003 *** 29059.0

These values are several orders of magnitude smaller than those used in the previous exercise, but the new procedure is taking a much longer time to complete. Why is that?

The answer is buried in a footnote to the expmod code (#46).

The reduction steps in the cases where the exponent e is greater than 1 are based on the fact that, for any integers x, y, and m, we can find the remainder of x times y modulo m by computing separately the remainders of x modulo m and y modulo m, multiplying these, and then taking the remainder of the result modulo m. For instance, in the case where e is even, we compute the remainder of be/2 modulo m, square this, and take the remainder modulo m. This technique is useful because it means we can perform our computation without ever having to deal with numbers much larger than m.

The important point is that the original expmod procedure uses successive squaring to perform its computations without ever having to deal with numbers larger than m. Simple arithmetic with very large numbers is much slower than arithmetic with 32-bit integers. That's because the time complexity for arithmetic operations is based on the number of bits in the operands. The intermediate results during computation in the fast-expt procedure get very big, very quickly (the final result is growing exponentially, after all). Since we're only interested in the remainder, modular arithmetic provides a much sleeker solution, as explained in the footnote.


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.24: The Fermat Test

From SICP section 1.2.6 Example: Testing for Primality

Exercise 1.24 asks us to once again modify the timed-prime-test procedure from exercise 1.22, this time using fast-prime? (which uses the Fermat test), then test the improvement using the timing statistics we gathered before.

The Fermat test is a probabilistic method for testing primality, meaning that if a given number passes the Fermat test, the best we can say is that the number is very probably prime. It's based on Fermat's Little Theorem which states:

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.

So to test n for primality, we select a number, a, which is less than n, and raise it to the nth power. We then divide aⁿ by n to get the remainder. If n is prime, then the remainder will be equal to a, the base we selected.

Fermat's Little Theorem says that this property will always hold true if n is prime, but it doesn't say anything about when n isn't prime. If aⁿ mod n is not equal to a, then we know for certain that n is not prime. However, there are values of n and a that will pass the Fermat test even when n is not prime. If we test enough values of a for a given n, we can increase our confidence that n is prime. Unfortunately, there are extremely rare values of n called Carmichael numbers that pass the Fermat test for any value of a. There are variations on the Fermat test that cannot be fooled. We'll look at one of those variations, the Miller-Rabin test, in a later exercise.

The complete code listing below shows what modifications were necessary to use the Fermat test in the timed-prime-test procedure. The fast-prime? procedure takes two arguments, n, and times. The first argument is the number to test, the second argument tells the procedure how many random values to test with. The more random values we use, the higher our confidence that n is prime, so we should test a lot of values. I've (rather arbitrarily) chosen to test 100 random values.

One other change that you may notice is the inclusion of a library module in the first line of code. When testing the original code listing from the book, I found that Scheme's primitive random procedure has a limit of 4,294,967,087. This won't do for the values we're testing, so I replaced it with the random-integer procedure, which doesn't have this limitation. I found out about this procedure and library through the Scheme Cookbook.

(require (lib "27.ss" "srfi"))

(define (timed-prime-test n)
   (newline)
   (display n)
   (start-prime-test n (current-inexact-milliseconds)))

(define (start-prime-test n start-time)
   (cond ((fast-prime? n 100)
          (report-prime (- (current-inexact-milliseconds) start-time)))))

(define (report-prime elapsed-time)
   (display " *** ")
   (display elapsed-time))

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

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

(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 (fermat-test n)
   (define (try-it a)
     (= (expmod a n n) a))
   (try-it (+ 1 (random-integer (- n 1)))))

(define (fast-prime? n times)
   (cond ((= times 0) true)
         ((fermat-test n) (fast-prime? n (- times 1)))
         (else false)))

Test Results

We need to run timed-prime-test again with the new modifications using the same set of inputs that we found in the previous exercises, then compare the results to see how much we improved the run time. The following table compares the original data we gathered in exercise 1.22, the improved algorithm we used in exercise 1.23 (improved time column), and the new results we got using the Fermat test (all three time columns are in milliseconds). The ratio column is the ratio of the improved values from exercise 1.23 to the fermat time column.

prime	original time	improved time	fermat time	ratio
10000000019	141	103	178	0.58
10000000033	172	101	187	0.54
10000000061	154	112	173	0.65
100000000003	516	310	181	1.71
100000000019	493	295	186	1.59
100000000057	527	305	189	1.61
1000000000039	1627	861	189	4.56
1000000000061	1559	884	195	4.53
1000000000063	1549	873	197	4.43
10000000000037	5014	2671	208	12.84
10000000000051	4932	2668	211	12.64
10000000000099	4855	2697	208	12.97
100000000000031	15745	8531	233	36.61
100000000000067	16022	8264	231	35.77
100000000000097	15861	8530	225	37.91
1000000000000037	48950	26346	244	108.0
1000000000000091	48836	26210	255	102.8
1000000000000159	49008	26256	257	102.2

Analysis

The timing numbers from the Fermat test start out looking pretty poor compared to what we've seen previously. This has mostly to do with the completely arbitrary number of random values I chose to test each prime with. As we increase the size of the numbers we're testing, you can see that the time using the Fermat test barely increases at all. We can verify that the time increase is logarithmic by observing that as the numbers under test increase by a factor of 10, the ratio column increases by a factor of roughly three. This logarithmic characteristic is due to the fact that the expmod procedure used in fermat-test has a logarithmic time complexity.


Related:

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