1 Introduction

Programming is one of the most difficult branches of applied mathematics; the poorer mathematicians had better remain pure mathematicians.

—Edsger Dijkstra

As we know, computers are simple machines that carry out a sequence of elementary steps, albeit very quickly. Unless you have a special-purpose processor, a computer can only compute addition, subtraction, multiplication, and division. If you think about it, you will see that the functions that might interest you when dealing with real or complex numbers can be built up from those four operations. We use many of these functions in nearly every program that we write, so we ought to understand how they are created.

We cannot expect the computer to solve integrals such as

ll− e−x2 dx |2 = π

in order to calculate π, and computing numerical integrals is more complicated than we want to attempt in this class. But, suppose that you wanted to calculate that integral, how would you do it? You might use a formula like

this:

[f (x0) + 2 ⌊n1⌋f (x2j) + 4 ⌊⌋ f (x2j −1) + f (xn)|

which is called the composite Simpson’s rule. If you were to compute 1 the error would be just 9.69935 × 10 − 12

Fortunately, you may recall from your Calculus class, with some conditions a function f (x) can be represented by its Taylor series expansion near some point f (a):

f (x) = f (a) + (x − a)k .

We will make extensive use of infinite series, suitably truncated, to solve many problems of interest.  Since we cannot compute an infinite series, we must be content to calculate a finite number of terms. In general, the more terms that we calculate, the more accurate our approximation. Note: when you see Σ or Π, you should generally think of a for loop.

If you have forgotten (or have never taken) Calculus, do not despair. Attend a laboratory section for review: the concepts required for this assignment are do not extend beyond derivatives.

2 Fundamental Constants

Transcendental numbers, they transcend the power of algebraic methods.

—Leonhard Euler

We live in a world (especially if you are Platonist) that is described by and perhaps governed by mathematics and mathematical objects. Our world is populated by numbers—fundamental constants—that we know to exist, but which we cannot write exactly using our decimal (or any positional) number system. We are thus presented with a pleasant conundrum: How do we calculate these numbers that need in order to pursue science?

There are many little things that remind us of the wonders of the physical world in which we live. One of the most beautiful things in mathematics, which Richard Feynman called“our jewel,”is Euler’s identity:

ei π + 1 = 0

which unites the most fundamental numbers in a single formula. A simple informal proof will show you why this is so.

The tool we will use for our proof is the Taylor series (Brook Taylor, 1685– 1731). Taylor showed that you can write most functions as an infinite sum (there are some restrictions). This allows us to evaluate a function f (x) near a point a by writing it like this:

f (x) = (x − a)k = f (a) + (x − a) + (x − a)2 + (x − a)3 + · · ·

where f \ is the first derivative, f \\ is the second, f (3) is the third, and so forth. For example,

x3 x5 x7 x9

3!     5!     7!     9!

and similarly,

ei x = 1 + i x − − + + − − + + − + · · ·

and look closely at the even numbered terms. You will see that they are the same as cos x, and if you look at the odd numbered terms, you will see that they are the same as i sin x (you simply have to factor out the i). Consequently,

we see that

ei x = cos x + i sin x.

If you let x = π, then

ei π = cos π + i sin π = − 1 + i × 0 = − 1.

Thus, if ei π = − 1 then ei π + 1 = 0. If you are like most scientists, you are left with a feeling of awe.

Can we make use of this jewel to calculate a value of π? We will start with ei π + 1 = 0 and subtract 1 from both sides, and we get ei π = − 1. We take the square root of both sides:

^ei π = ^− 1

and we get e = i . We take the ith root of both sides:

^ie = ^i i

which simplifies to e = ^i i . We are almost finished once we take the natural logarithm of both sides: log(e ) = log(^i i) yielding = log(^i i) which simplifies to π = 2 log(^i i) = 2 log(i−i). So ultimately we find that π is the loga- rithm of the ith root of an imaginary number. We may want to contemplate the words of Benjamin Peirce (1809– 1880), a professor of mathematics at Harvard, who said“Gentlemen, that is surely true, it is absolutely paradoxical; we cannot understand it, and we don’t know what it means. But we have proved it, and therefore we know it is the truth.”

3 Calculating e

Let us change our traditional attitude to the construction of programs. Instead of imagining that our main task is to instruct a computer what to do, let us concentrate rather on explaining to human beings what we want a computer to do.

—Donald Knuth

The number e, also known as Euler’s number (Leonhard Euler, 1707– 1783), is an irrational mathematical constant approximately equal to 2.71828, that appears pervasively in the natural and mathematical worlds. It is the base of the natural logarithm, it is the limit of limn一∞ (1 + )n which was discovered by Jacob Bernoulli in his work on the calculation of compound interest. And, of course, it can be expressed as the Taylor series:

∞ 1            1     1     1      1        1         1          1            1              1                1

e = k! = 1 + 1 + 2 + 6 + 24 + 120 + 720 + 5040 + 40320 + 362880 + 3628800 + · · ·

How many terms must you compute? Fewer than you might expect, since k! grows very fast. You will be deter- mining that experimentally as part of this assignment.

If we are naïve about computing the terms of the series we can quickly get into trouble — the values of k! get large very quickly. We can do better if we observe that:

xk xk−1 x

= ×

k!     (k − 1)! k .

At first, that looks like a recursive definition (and in fact, you could write it that way, but it would be wasteful). As we progress through the computation, assume that we know the previous result. We then just have to compute the next term and multiply it by the previous term. At each step we just need to compute , starting with k = 0! (remember 0! = 1) and multiply it by the previous value and add it into the total. It turns into a simple for or while loop.

25

20

15

10

5

2                                          4                                          6                                          8                                          10

Figure 1: Comparing      for x = 2, 3, 4, 5.

4 Calculating π

The best programs are written so that computing machines can perform them quickly and so that human beings can understand them clearly. A programmer is ideally an essayist who works with traditional aesthetic and literaryforms as well as mathematical concepts, to communicate the way that an algorithm works and to convince a reader that the results will be correct.

—Donald Knuth

Why, you might ask, do we need to calculate π?  In practice, we do not: it is available as part of the of the πr2 and the volume of a sphere is πr3 . The cosmological constant is

8πG

Λ = ρ .

Heisenberg’s uncertainty principle is given by

4π .

Einstein’s general relativity field equation is

Rµν − gµν+ Λgµν = Tµν ,

and so forth. No matter where you look, π is pervasive in the physical world.

So why do we calculate it?  Well, suppose the person who typed up the π fall into the area of experimental mathematics.

Experimental mathematics is an approach to mathematics in which computation is used to investigate mathe- matical objects and identify properties and patterns. It has been defined in a discussion by J. Borwein, P. Borwein, R. Girgensohn and S. Parnes as“that branch of mathematics that concerns itself ultimately with the codification

and transmission of insights within the mathematical community through the use of experimental (in either the Galilean, Baconian, Aristotelian or Kantian sense) exploration of conjectures and more informal beliefs and a care- ful analysis of the data acquired in this pursuit.”

In the subsequent sections, we will present a number of functions p of n that are series that can be used to approximate π . In each of these formulæ there is a p(n) where the limit

lim p(n) = π .

We know, both through experiment (which usually comes first) and mathematical proof, that these do result in π . The question then is: Which of these compute it most efficiently? Efficiency in this case is measured in the number of terms (or factors, in the case of Wallis and Viète).

4. 1 The Leibniz Formula

The Leibniz formula for π (Gottfried Wilhelm von Leibniz, 1646– 1716), given by

p(n) = 4 k(−)1(k) = 4 ( 1 − + − + − ··· ) = ((−1)n H + − H − ) + π

converges extremely slowly, since the error term involves harmonic numbers, and consequently is not reasonable for calculation. It is no more than the Taylor series expansion of tan−1 x with x = 1.

4.2 The Madhava Series

The Madhava series (Mdhava of Sangamagrma, c. 1340 – c. 1425) is also related to tan−1 x,

1(k) = ^3tan− 1 =

but it gives us a more rapidly converging series, given by:

p(n) = ^12 1(k) = ^12 l 3−n−1 ( (−1)n Φ ( − , 1, n + ) + π3n+ )| .

Which leaves us with the problem of calculating ^12. You may be thinking:“Can’t I just use the library?”  The answer, of course, is no, you will need to compute ^12 on your own. Do you need to worry about Φ? No, you just need to understand that in the limit, the Φ term (called the Lerch transcendent) goes to zero and that the remaining

term π 3 −n−13n+ = π = π

4.3 The Wallis Series

John Wallis (1616– 1703) was an English clergyman and mathematician who is given partial credit for the develop- ment of infinitesimal calculus. He gave us a series that is purely multiplicative (instead of terms, we have factors):

n 4k2 πΓ(n + 1)2

It is easy to calculate, but how rapidly does it converge? What is this Γ function? Do we need to compute it? How do we know this converges to π? Let’s factor out π, then take the limit and note that

Γ(1 + n)2

n一∞ Γ ( + n) Γ ( + n)

This tells us that as you compute more terms, you get closer to the true value of π .

π - Leibniz n

0.010

0.005

Figure 2: Leibniz sequence 4 zk(n)=0 1(k)

π - Madhava n

25

-5. × 10-10

-1. × 10-9

-1.5 × 10-9

Figure 3: Madhava sequence ^12 zk(n)=0 (k(3)k

4.4 Euler’s Solution

The Basel problem is a problem in mathematical analysis with relevance to number theory, first posed by Pietro Mengoli in 1650 and solved by Leonhard Euler in 1734. The Basel problem asks for the precise summation of the reciprocals of the squares of the natural numbers, i.e. the sum of the infinite series

= + + + · · · = H ,

which again involves harmonic numbers. Euler’s solution showed that the solution is π2 /6, but his method gave

us this series:

p(n) = 46

which also requires us to calculate the square root, of this case an unknown (until we calculate it) number.

4.5 The Bailey-Borwein-Plouffe Formula

The Bailey-Borwein-Plouffe formula (BBP formula) is a formula for π . It was discovered in 1995 by Simon Plouffe and is named after the authors of the article in which it was published, David H. Bailey, Peter Borwein, and Plouffe. The formula that they discovered is remarkably simple:

p(n) = 16 −k( − − − ).

And if you desire to reduce it to the least number of multiplications, you can rewrite it in Horner normalform:

p(n) =z(n) 16 −k × (k(120k + 151) + 47)

4.6 Viète’s Formula

Named after François Viète, Viète’s formula is a infinite product of nested radicals that can be used for calculations of π, though it should be noted that methods found before this specific formula are known to produce greater accuracy.

Viète’s formula can be written as follows:

2 ^2 ^2 + ^2 ^2 + ^2 + ^2

=                ×                                   × ···

π 2           2                    2

Or more simply,

2 =u(∞) ak

π k = 1  2

where a1 = ^2 and ak = ^2 + ak−1 for all k > 1.