Computing Pi again

November 15, 2023

Some time ago, I wrote about computing Pi using dedicated approximating formulae and continued fractions. TIL that the Corman Lisp distribution has an example of computing $\Pi$ to any number of decimal places. The code (MIT license) is shown below:

(defun pi-atan (k n)
	(do* ((a 0)
      	  (w (* n k))
		  (k2 (* k k))
		  (i -1))
		((= w 0) a)
		(setq w (truncate w k2))
		(incf i 2)
		(incf a (truncate w i))
		(setq w (truncate w k2))
		(incf i 2)
		(decf a (truncate w i))))

(defun calc-pi (digits)
	(let* ((n digits)
		   (m (+ n 3))
		   (tenpower (expt 10 m)))
					(+ (pi-atan 18 (* tenpower 48))
					   (pi-atan 57 (* tenpower 32)))
					(pi-atan 239 (* tenpower 20)))

It works, and it’s kinda fast for $n < 10^5$!

CL-USER(4): (calc-pi 1000)


As you may already have guessed, it relies on the use of arctan, which we covered in the aforementioned post.

The following expansion is used: $\sum_{n=0}^{\infty}\frac{(-1)^nx^{2n+1}}{2n+1}$ for $-1<x<1$, and we look for combination of arctan such that $\text{arctan}(1)=\pi/4$. In other words, the formula can be rewritten as:

$$ \pi = 4\sum_{n=0}^{\infty}\frac{(-1)^nx^{2n+1}}{2n+1}\left( \text{something} \right). $$

Smaller “something” means faster convergence.1

There are some magic numbers in the above code: 18, 57, etc. Where do they come from? There are many formulae, but this one is from Gauss:

$$ \pi = \sum_{n=0}^{\infty}\frac{(-1)^n}{k}\left( 48\left(\frac{1}{18}\right)^k + 32\left(\frac{1}{57}\right)^k - 20\left(\frac{1}{239}\right)^k\right), $$

with $k=2n+1$. Other formulae are available on this French blog post.

Note that in the GitHub repository, you will also find some sieves for prime numbers.

♪ Monodrama • Sarabande

  1. Machin’s formula ($\pi/4 = 4\text{arctan}(1/5) - \text{arctan}(1/239)$) is probably the easiest to use as it amounts to $\pi = \sum_{n=0}^{\infty}\frac{(-1)^n}{2n+1}\left( 4\left(\frac{1}{5}\right)^{2n+1} - \left(\frac{1}{239}\right)^{2n+1}\right)$. ↩︎

See Also

» Generating power sets in Lisp » Trying out cl-lsp » Welford algorithm in Newlisp » Fitting an OLS model in Newlisp » Plotting with Newlisp