Computing pi
This article goes over the algorithms and some optimizations techniques to calculate pi that I implemented out of curiosity and purely for learning purposes.
For the source code see code repository.
For highly optimized software written by others see y-cruncher, which is actually used for world record computations.
Bailey–Borwein–Plouffe formula
The Bailey–Borwein–Plouffe formula is given by:
This formula is particularly useful when calculating in hexadecimal
or binary, as the
factor guarantees the amount of correct digits
up to the cutting point of the sum. Adding an additional term may be needed to
deal with carry overs.
The implementation provided in the repository is quite straightforward. Thus, it is not further explained here.
Binary Splitting
Even if the computer used for calculations supports multiplications in constant time, for numbers beyond the size of a computer register bignums must be used, for which the computation time depends in the size of the number. To reduce the average size of the integers multiplied in each step, on can use the binary splitting technique instead of a direct term-by-term summation.
We consider sums of the form:
that can be expressed calculated as the ratio of some and
, where the terms
,
,
are
integer functions of
.
If we expand the sum and products we have
and find the following expressions for the values of ,
and
Let now be an integer such that
. Then
and we have the recursion formula for ,
and
given
by
with the base cases
A simple implementation of binary splitting in software would be
(define (binary-splitting n base-case) (define (iteration a b) (cond [(= a (- b 1)) (base-case b)] [(> b a) (let ([m (fxdiv (fx+ a b) 2)]) (let-values ([(p1 q1 r1) (iteration a m)] [(p2 q2 r2) (iteration m b)]) (values (+ (* p1 q2) (* p2 r1)) (* q1 q2) (* r1 r2))))] [else (error 'iteration "a >= b")])) (iteration 0 n))
where base-case is a procedure that computes (values p q r).
Additional information:
Chudnovsky algorithm
The Chudnovsky algorithm uses the following series for calculating :
Our goal is to rewrite this formula into the form of (2).
We start by separating the term and noting that:
Thus, we have
For generating the product inside each term, we recall some properties of the product, here applied to the cases relevant to our series:
Hence, we may rewrite part of each term as
and the values of ,
and
can be identified as
A final issue is the square root in front of the sum, as we do not wish to introduce floating point arithmetic and limit the numbers of digits of pi we can compute. To solve this, we observe first that
and then recall the useful Taylor series
for . We can also evaluate this series by using the binary
splitting technique. For this, we rewrite the series as
and directly identify the corresponding ,
and
:
Additional information:
