For pi day: A continued fraction expansion of pi

6

Many geeky mathematical people celebrate "pi day" on March 14, because the date is written 3/14 in the US, which is evocative of the decimal representation of π = 3.14....

Most people are familiar with the decimal representation of π. The media occasionally reports on a new computational tour-de-force that approximates π to unimaginable precision. The decimal value of π has been computed to more than 10 trillion digits!

Using the decimal approximation of π is ubiquitous in mathematics, but in elementary school I learned another approximation: π ≈ 22/7. As I got older, I discovered that there were other rational approximations of π, such as 355/113.

Continued fractions

To celebrate π day, I am going to use SAS to estimate π by using its simple continued fraction representation. Every real number x can be represented as a continued fraction:

In this continued fraction, the ai are positive integers. Because all of the fractions have 1 in the numerator, the continued fraction can be compactly represented by specifying only the integers: x = [a0; a1, a2, a3, ...]. Every rational number is represented by a finite sequence; every irrational number requires an infinite sequence. Notice that the decimal representation of a number does NOT have this nice property: 1/3 is rational but its decimal representation is infinite.

The continued fraction expansion of pi

The continued fraction expansion of π is known to more than 15 billion values. The first fifteen digits are
π ≈ [3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1]
What does this mean? Well, the sequence [3; 7] corresponds to the fraction 3 + 1/7 = 22/7. Similarly, the sequence [3; 7, 15, 1] corresponds to the fraction 3 + 1/(7 + 1/(15 + 1/1)) = 355/113.

Trying to compute these fractions by hand is tedious. Why not write a SAS program to evaluate the continued fraction sequence? As the examples show, the computation starts at the end of the sequence. The computation begins by taking the reciprocal of the last number in the sequence. That values is added to the second-to-last number, and the sum is inverted. This process continues until you reach the first value. You can do this computation in the SAS DATA step (post your code in the comments!), but I have written a SAS/IML function to evaluate an arbitrary continued fraction representation. Notice that in SAS you can call the CONSTANT function to obtain a numerical approximation of π and other mathematical constants, which is useful for checking the results:

proc iml;
start ContinuedFraction(f);
   s = f[nrow(f)];              /* last value */
   do i = nrow(f)-1 to 1 by -1; /* traverse sequence in reverse order */
      s = f[i] + 1/s;           
   end;
   return( s );
finish;
 
/* test the function by using the continued fraction expansion for pi */
f = {3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1};
s = ContinuedFraction(f);
pi = constant("pi");    /* call the CONSTANT function to approximate pi */
print s[format=16.14], pi[format=16.14];

As you can see, the first fifteen terms of the continued fraction expansion of π are sufficient to reproduce the first 15 digits of the decimal representation.

Other continued fractions

Here are a few continued fraction expansions for some other irrational numbers:

  • The continued fraction expansion for e ≈ 2.718281828459 shows a fascinating regular pattern: two 1s followed by an even number.
    e = [2; 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1, 1, 10,....]
  • The continued fraction expansion for the square root of 2 ≈ 1.414213562373095 contains a repeating sequence of 2s:
    sqrt2 = [1; 2, 2, 2, 2, ...]
  • The continued fraction expansion for the golden ratio φ ≈ 1.61803398874989 is a repeating sequence of 1s:
    phi = [1; 1, 1, 1, 1, 1, ...]
You can use these sequences to estimate the corresponding decimal representation:
e     = ContinuedFraction({2, 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1, 1, 10});
sqrt2 = ContinuedFraction({1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2});
phi   = ContinuedFraction({1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1});
print e[format=16.14], sqrt2[format=16.14], phi[format=16.14];

The results contain less than 15-digit decimal precision. In order to obtain more precision, you would need to use more terms in the continued fraction representation.

By the way, next year's Pi Day should be special: celebrate on 3/14/15 at 9:26:53 for maximal enjoyment!

Share

About Author

Rick Wicklin

Distinguished Researcher in Computational Statistics

Rick Wicklin, PhD, is a distinguished researcher in computational statistics at SAS and is a principal developer of PROC IML and SAS/IML Studio. His areas of expertise include computational statistics, simulation, statistical graphics, and modern methods in statistical data analysis. Rick is author of the books Statistical Programming with SAS/IML Software and Simulating Data with SAS.

6 Comments

  1. I once tried to get 355/113 as a vanity license plate, but it got rejected. :(
    I thought it would make a deep statement about perfection vs. good enough.

    • Rick Wicklin

      I think it is fascinating that 22/7 and 355/113 both have prime numbers as denominators. So does a later approximation, 833719/265381. Question: Are there infinitely many prime numbers in the denominator of the continued fraction series for pi?

    • gordon Keener on

      I know someone at SAS, in JMP, with the plate "E**IPI&1". I believe '+' was not an allowed character when she got it.

      Tangent: A friend of mine had "/DEV/CAR".

  2. Pingback: Analyzing the first 10 million digits of pi: Randomness within structure - The DO Loop

Leave A Reply

Back to Top