I just "discovered" an interesting interpolation for the

Wallis product formula for pi. I don't think this deserves a new thread, and this seems like a good place for it.

The Wallis product formula is

[math]\prod_{n=1}^{\infty} \frac{(2n)(2n)}{(2n-1)(2n+1)} = \frac{2}{1} \cdot \frac{2}{3} \cdot \frac{4}{3} \cdot \frac{4}{5} \cdot \frac{6}{5} \cdot \frac{6}{7} \cdot \frac{8}{7} \cdot \frac{8}{9} \cdots = \frac{\pi}{2}[/math] Note that successive partial products alternate above and below the limit. So it's likely that interpolating between them will give a good pi approximation.

Also note that terms can be grouped into pairs in two ways:

[math]\frac{n^2}{n^2-1} \mbox{ for even n}[/math] or

[math]\frac{n^2-1}{n^2} \mbox{ for odd n}[/math]

The second form is useful because we can divide n² out to give (1-1/n²).

I have found that by calculating upto a term of the form

[imath]\frac{2n}{2n-1}[/imath] & interpolating between the last two products, that a much better approximation to pi can be found. Let p be a partial product with the next term being

[imath]r = \frac{m+2}{m+1}[/imath]Next we find a & g, the arithmetic & geometric means respectively of p & pr.

Then a - g is very close to pi - a. So 2a - g is a good pi approximation.

Here's a short bc program that illustrates this technique. Note that m

must be even for this to work properly.

Code: Select all

`scale=30;m=5000;p=4*m/(m+1);for(i=3;i<m;i+=2)p-=p/(i*i);r=(m+2)/(m+1);p*(1+r-sqrt(r))`

The above code yields 3.141592653589793630690807651016

Proof of why this method works will be left as an exercise for the reader.

Of course, this algorithm is rather slow compared to the Arithmetic-Geometric Mean algorithm (and other formulae discovered by Ramanujan), or even some good arctan-based formulae, but apart from one square root most of the arithmetic is very simple.