Computing pi

Collection of algorithm implementations for computing digits of 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:

(1)

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:

(2)

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 :

(3)

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: