michael@0: This file describes how pi is computed by the program in 'pi.c' (see
michael@0: the utils subdirectory).
michael@0: 
michael@0: Basically, we use Machin's formula, which is what everyone in the
michael@0: world uses as a simple method for computing approximations to pi.
michael@0: This works for up to a few thousand digits without too much effort.
michael@0: Beyond that, though, it gets too slow.
michael@0: 
michael@0: Machin's formula states:
michael@0: 
michael@0: 	 pi := 16 * arctan(1/5) - 4 * arctan(1/239)
michael@0: 
michael@0: We compute this in integer arithmetic by first multiplying everything
michael@0: through by 10^d, where 'd' is the number of digits of pi we wanted to
michael@0: compute.  It turns out, the last few digits will be wrong, but the
michael@0: number that are wrong is usually very small (ordinarly only 2-3).
michael@0: Having done this, we compute the arctan() function using the formula:
michael@0: 
michael@0:                        1      1       1       1       1     
michael@0:        arctan(1/x) := --- - ----- + ----- - ----- + ----- - ...
michael@0:                        x    3 x^3   5 x^5   7 x^7   9 x^9
michael@0: 
michael@0: This is done iteratively by computing the first term manually, and
michael@0: then iteratively dividing x^2 and k, where k = 3, 5, 7, ... out of the
michael@0: current figure.  This is then added to (or subtracted from) a running
michael@0: sum, as appropriate.  The iteration continues until we overflow our
michael@0: available precision and the current figure goes to zero under integer
michael@0: division.  At that point, we're finished.
michael@0: 
michael@0: Actually, we get a couple extra bits of precision out of the fact that
michael@0: we know we're computing y * arctan(1/x), by setting up the multiplier
michael@0: as:
michael@0: 
michael@0:       y * 10^d
michael@0: 
michael@0: ... instead of just 10^d.  There is also a bit of cleverness in how
michael@0: the loop is constructed, to avoid special-casing the first term.
michael@0: Check out the code for arctan() in 'pi.c', if you are interested in
michael@0: seeing how it is set up.
michael@0: 
michael@0: Thanks to Jason P. for this algorithm, which I assembled from notes
michael@0: and programs found on his cool "Pile of Pi Programs" page, at:
michael@0: 
michael@0:       http://www.isr.umd.edu/~jasonp/pipage.html
michael@0: 
michael@0: Thanks also to Henrik Johansson <Henrik.Johansson@Nexus.Comm.SE>, from
michael@0: whose pi program I borrowed the clever idea of pre-multiplying by x in
michael@0: order to avoid a special case on the loop iteration.
michael@0: 
michael@0: ------------------------------------------------------------------
michael@0:  This Source Code Form is subject to the terms of the Mozilla Public
michael@0:  # License, v. 2.0. If a copy of the MPL was not distributed with this
michael@0:  # file, You can obtain one at http://mozilla.org/MPL/2.0/.