Sunday, March 14, 2010

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 / r2

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 310/71 and 31/7. 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.


The first 50 digits of π are:


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 π.

Thursday, March 4, 2010

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,

an ≡ a (mod n)

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

an-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 an-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 "" "srfi"))

(define (square-check x m)
(if (and (not (or (= x 1) (= x (- m 1))))
(= (remainder (* x x) m) 1))
(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))
(remainder (* base (expmod base (- exp 1) 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)
> (miller-rabin-test 1105)
> (miller-rabin-test 1729)
> (miller-rabin-test 2465)
> (miller-rabin-test 2821)
> (miller-rabin-test 6601)

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


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