Calculating Pi

From Antiquity to Modern Times

by George Beck and Michael Trott

Why is π the same in C = 2  π r and A = π (r^2)?

We slice up a circle like a pie and rearrange the sectors to form a figure that is roughly a rectangle. As the sectors become thinner, the rectangle's left and right sides become more vertical and the top and bottom become flatter. We can imagine that in the limit we get an actual rectangle.

Table[PiRectangle[2  n], {n, 2, 10}];

The height of the rectangle equals the radius of the circle, r. Because the top and bottom of the rectangle are made up of all the small arcs that add up to the circumference of the circle, the rectangle has length C/2. Since the rectangle has the same area as the circle, the length times the height is C r/2 = π (r^2). Dividing by r and multiplying by 2 gives C = 2  π r.

PiRectangle[100];

Approximating π Using Areas of Polygons—the Method of Archimedes

Sandwiching π

A circle of radius r has area π (r^2), so the area of a circle of radius 1 is π. Therefore π is sandwiched between the areas of the regular n-sided polygons inscribed and circumscribed around the circle. PiPolygons[n] draws the three figures and prints the corresponding inequality for the areas.

PiPolygons[5];

Finding the Area of the Inscribed Polygon

Break the inscribed polygon into n triangles of height Cos[a/2] and base 2 Sin[a/2], where a is the angle 2 π/n.

The area is \(TraditionalForm
      \`n\ ×\ 1/2\ ×\ base\ ×\ height\  = \ 
        \(n/2\ \((2\ Sin[a/2]\ Cos[a/2])\) = n/2\ Sin[a]\)\).

Finding the Area of the Circumscribed Polygon

For the circumscribed polygon, the n triangles each have height 1 and base 2 Tan[a/2]. The area is  \(TraditionalForm
      \`n×1/2\ ×2×Tan[a/2] = n\ Tan[a/2]\).

PiTriangles[7];

Vieta's Formula for π

We continue the discussion on approximating π by calculating the areas of a circle using inscribed and circumscribed regular polygons. We illustrate Vieta's formula, developed in 1593, the oldest exact result derived for π.

The Formula

Vieta's formula expresses (2/π) as an infinite product of nested square roots.

Quality of Approximation

The graph in the middle is the approximation from Vieta's formula, which calculates the area of a regular polygon of (2^(n + 1)) sides. In contrast, the upper and lower graphs show the estimates to π based on the areas of the n-sided regular polygons circumscribed and inscribed around a circle of radius 1.

PiSqueeze[100];

The second column shows the number of sides of the regular polygons, while the third shows the approach to π. The numbers in the fourth column are the differences between π and the partial products of the infinite product.

TableForm[
      Table[{n, (2^(n + 1)), 
          tmp = (2/NProduct)[Cos[(π/(2^i))], {i, 2, n + 1}], π - tmp},  
        {n, 15}]]

Why is It True?

To derive Vieta's formula we need to combine several results from calculus and trigonometry. We will illustrate these results without proof.

Finding a Limit

When x is large, Sin[x] and x are rather different.

Plot[{x, Sin[x]}, {x, -10, 10}, AspectRatio -> Automatic]

But when x is small, Sin[x] and x are approximately the same.

Plot[{x, Sin[x]}, {x, -1, 1}, AspectRatio -> Automatic]

Mathematically, this is expressed in terms of limits.

If t is any fixed number, (t/(2^n)) tends to 0 as n gets large.

Limit[(t/(2^n)), n -> ∞]

So sin[t/(2^n)]/(t/(2^n)) will tend to 1.

Limit[Sin[(t/(2^n))]/((t/(2^n))), n -> ∞]

Therefore t sin[t/(2^n)]/(t/(2^n))tends to t.

Limit[t Sin[(t/(2^n))]/((t/(2^n))), n -> ∞]

This limit can be rewritten more simply.

Limit[(2^n) Sin[(t/(2^n))], n -> ∞]

Applying a Trigonometric Rule

We apply and keep reapplying the following rule to Sin[t].

NestList[(# /. SinRule) &, Sin[t], 5] // TableForm

Each of these products is equal to Sin[t]. We showed earlier that the product of the first and last terms in such a product, (2^n) Sin[t/(2^n)], tends to t as n tends to infinity.

Taking the limit and dividing by t, we get a formula of Euler:

We substitute (π/2) for t.

% /. t -> (π/2)

Applying Another Trigonometric Rule

All that is left to do is to verify that the cosines are nested square roots. Mathematically this would be proved by induction. We illustrate the result for the fifth term by repeatedly applying another trigonometric rule.

CosRule = Cos[z_] :> (Sqrt[(1/2) (Cos[2 z] + 1)]);

Cos[(π/(2^6))] //. CosRule

This simplifies to the form of nested square roots at the very beginning.

Simplify[%]

Computing π with Series

A Classical Sum

This is Leibniz' series, developed in 1674.

(∑_(j = 0)^∞)((((-1))^j)/(2 j + 1))

Summing up to an finite limit gives an approximation for π. Here are the differences between π/4 and the first 15 partial sums.

TableForm[
      Table[{n, 
          N[(π/4) - (∑_(j = 0)^n)((((-1))^j)/(2 j + 1))]}, 
        {n, 1, 15}]]

Two Elegant Series that Converge Slowly

This series is due to Euler, developed in 1748.

(∑_(j = 1)^∞) (1/(j^2))

The series doesn't converge particularly fast, so it is not useful in practice for calculating π.

TableForm[
      Table[{n, N[π - (Sqrt[6 (∑_(j = 1)^n) (1/(j^2))])]}, 
        {n, 1, 15}]]

Here is a similar series.

(∑_(j = 0)^∞) (1/(((2 j + 1))^2))

This series also converges too slowly to be used to calculate many digits of π.

TableForm[
      Table[{n, 
          N[π - (Sqrt[8 (∑_(j = 0)^n) (1/(((2  j + 1))^2))])]}, 
        {n, 1, 15}]]

A Series Based on (tan^(-1))

Historically the calculations of the first few tens of digits of π were based on various identities involving the (tan^(-1)) function. Here is such an identity, developed by Machin in 1706.

(tan^(-1)) has the following series representation.

This sum converges much faster.

TableForm[
      Table[{n, 
          N[N[π - 
                (16 (∑_(j = 0)^n)((((-1))^j)/(2 j + 1)) 
                          ((((1/5)))^(2  j + 1)) - 
                    4 (∑_(j = 0)^n)((((-1))^j)/(2 j + 1)) 
                          ((((1/239)))^(2  j + 1))), 22]]}, {n, 1, 15}]]

Another (tan^(-1)) Identity

In the spirit of the last series, Lehmer in 1938 used the following identity to calculate π .

This sum converges even better.

TableForm[
      Table[{n, 
          N[N[π - 
                (88 (∑_(j = 0)^n)((((-1))^j)/(2 j + 1)) 
                          ((((1/28)))^(2  j + 1)) + 
                    8 (∑_(j = 0)^n)((((-1))^j)/(2 j + 1)) 
                          ((((1/443)))^(2  j + 1)) - 
                    20 (∑_(j = 0)^n)((((-1))^j)/(2 j + 1
                                )) ((((1/1393)))^(2  j + 1)) - 
                    40 (∑_(j = 0)^n)((((-1))^j)/(2 j + 1
                                )) ((((1/11018)))^(2  j + 1))), 
              22]]}, {n, 1, 15}]]

A Formula Based on Modular Forms

In 1914 Ramanujan gave the following formula for π based on modular identities.

The sum of the series can be expressed in terms of hypergeometric functions.

((2 (Sqrt[2]))/9801) (∑_(j = 0)^∞)(((4 j)! (26390 j + 1103))/(((4^(4 j)) ((j!)^4)) (99^(4 j))))

This series converges much faster than the previous ones, making it necessary to increase the value of $MaxExtraPrecision.

TableForm[
      Table[{n, 
          N[N[π - 
                (1/(((2 (Sqrt[2]))/9801) (∑_(j = 0)^n)(((4 j)! (26390 j + 1103))/(((4^(4 j)) ((j!)^4)) (99^(4 j)))))), 22]]}, {n, 1, 15}]]

Another Formula Based on Modular Forms

In 1989 D. V. Chudnovsky and G. V. Chudnovsky gave another formula for π based on modular identities.

This sum converges even faster than the previous one. Each new term gains roughly 14 more digits of accuracy.

TableForm[
      Table[{n, 
          N[N[π - 
                (1/(12 (∑_(j = 0)^n)(((((-1))^j) (6 j)! (545140134 j + 13591409))/(((j!)^3) (3 j)! ((((640320^3)))^(j + (1/2))))))), 22]]}, {n, 1, 10}]]

Computing π by Iterative Processes

An Iteration Formula

In 1989 Borwein expressed π  as the limit of a sequence defined by an iteration, (lim_(n -> ∞))(Subscript[a, n]) = 1/π.

Here the first term of the sequence.

(1/(Subscript[a, 1])) // Simplify

It is an algebraic number that is a root of the following polynomial.

RootReduce[%][[1]][x]

The series approximations for π give roughly a fixed number of digits per term. In comparison, the number of correct digits in iteration formulas grows quadratically with the iteration number.

TableForm[N[Table[(1/(Subscript[a, n])) - π, {n, 0, 3}], 22]]

The quantity \(TraditionalForm
      \`|1\/π - a\_n|\) is bounded by \(TraditionalForm
      \`16\ 4\^n\ e\^\(\(-2\)\ π\ 4\^n\)\).

TableForm[
      N[Table[16 (4^n) (E^(-2 π (4^n))), {n, 0, 5}]]]

An Improved Iteration Formula

In 1996 Bailey gave an improved iteration formula for the calculation of π. Again the result is that  (lim_(n -> ∞))(Subscript[a, n]) = 1/π.

Here is the first term of this sequence.

(1/(Subscript[a, 1])) // Simplify

It is an algebraic number that is a root of the following polynomial.

RootReduce[%][[1]][x]

Here is how close the sequence is to π for the first few iterations.

TableForm[N[Table[(1/(Subscript[a, n])) - π, {n, 0, 3}], 22]]

The quantity \(TraditionalForm
      \`|1\/a\_n - π|\) is bounded by \(TraditionalForm
      \`16\ 5\^n\ e\^\(\(-π\)\ \ 5\^n\)\).

TableForm[
      N[Table[16 (5^n) (E^(-π  (5^n))), {n, 0, 5}]]]

Implementation

Usage Messages

Function Definitions