1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/security/nss/lib/freebl/mpi/doc/pi.txt Wed Dec 31 06:09:35 2014 +0100
1.3 @@ -0,0 +1,53 @@
1.4 +This file describes how pi is computed by the program in 'pi.c' (see
1.5 +the utils subdirectory).
1.6 +
1.7 +Basically, we use Machin's formula, which is what everyone in the
1.8 +world uses as a simple method for computing approximations to pi.
1.9 +This works for up to a few thousand digits without too much effort.
1.10 +Beyond that, though, it gets too slow.
1.11 +
1.12 +Machin's formula states:
1.13 +
1.14 + pi := 16 * arctan(1/5) - 4 * arctan(1/239)
1.15 +
1.16 +We compute this in integer arithmetic by first multiplying everything
1.17 +through by 10^d, where 'd' is the number of digits of pi we wanted to
1.18 +compute. It turns out, the last few digits will be wrong, but the
1.19 +number that are wrong is usually very small (ordinarly only 2-3).
1.20 +Having done this, we compute the arctan() function using the formula:
1.21 +
1.22 + 1 1 1 1 1
1.23 + arctan(1/x) := --- - ----- + ----- - ----- + ----- - ...
1.24 + x 3 x^3 5 x^5 7 x^7 9 x^9
1.25 +
1.26 +This is done iteratively by computing the first term manually, and
1.27 +then iteratively dividing x^2 and k, where k = 3, 5, 7, ... out of the
1.28 +current figure. This is then added to (or subtracted from) a running
1.29 +sum, as appropriate. The iteration continues until we overflow our
1.30 +available precision and the current figure goes to zero under integer
1.31 +division. At that point, we're finished.
1.32 +
1.33 +Actually, we get a couple extra bits of precision out of the fact that
1.34 +we know we're computing y * arctan(1/x), by setting up the multiplier
1.35 +as:
1.36 +
1.37 + y * 10^d
1.38 +
1.39 +... instead of just 10^d. There is also a bit of cleverness in how
1.40 +the loop is constructed, to avoid special-casing the first term.
1.41 +Check out the code for arctan() in 'pi.c', if you are interested in
1.42 +seeing how it is set up.
1.43 +
1.44 +Thanks to Jason P. for this algorithm, which I assembled from notes
1.45 +and programs found on his cool "Pile of Pi Programs" page, at:
1.46 +
1.47 + http://www.isr.umd.edu/~jasonp/pipage.html
1.48 +
1.49 +Thanks also to Henrik Johansson <Henrik.Johansson@Nexus.Comm.SE>, from
1.50 +whose pi program I borrowed the clever idea of pre-multiplying by x in
1.51 +order to avoid a special case on the loop iteration.
1.52 +
1.53 +------------------------------------------------------------------
1.54 + This Source Code Form is subject to the terms of the Mozilla Public
1.55 + # License, v. 2.0. If a copy of the MPL was not distributed with this
1.56 + # file, You can obtain one at http://mozilla.org/MPL/2.0/.
