Beating Bellard's Formula
By Robert Smith
Fabrice Bellard came up with a computationally efficient formula for calculating the nth hexadecimal digit of $\pi$ without calculating any of the previous n−1. It’s called Bellard’s formula. It wasn’t the first of its kind, but in terms of computational efficiency, it was a substantial improvement over the original, elegant Bailey-Borwein-Plouffe formula. Due to the trio’s discovery, these formulas are often called BBP-type formulas.
Over the years, numerous BBP-type formulas have been discovered. In fact, Bailey gives us a recipe to search for them using integer-relation algorithms. In simple terms, we can just guess formulas, and run a computation to see if it likely equals $\pi$ with high confidence. If we do find one, then we can use it as a conjecture to prove formally.
Like Bellard and many others, I ran a variant of Bailey’s recipe, effectively doing a brute-force search, highly optimized and in parallel. The search yielded another formula that is computationally more efficient than Bellard’s formula. The identity is as follows:
$$ \pi = \sum_{k=0}^{\infty} \frac{1}{4096^k} \left( \frac{1}{6k+1} - \frac{2^{-5}}{6k+3} + \frac{2^{-8}}{6k+5} + \frac{2}{8k+1} - \frac{2^{-5}}{8k+5} + \frac{2^{-1}}{12k+3} - \frac{2^{-4}}{12k+7} - \frac{2^{-8}}{12k+11} \right). $$
It converges at a rate of 12 bits per term. We will prove convergence, and then prove the identity itself (with a little computer assistance). As it turns out, an equivalent form of this formula was already discovered, which we will discuss as well. Finally, we’ll show a very simple implementation in Common Lisp.
Proof of convergence
Write the series as $S := \sum_{k=0}^{\infty} 4096^{-k}R(k)$. Since $R(k)\in O(1/k)$, convergence is dominated by the geometric term $4096^{-k}$:
$$ \lim_{k \to \infty} \left\vert \frac{R(k+1)}{4096^{k+1}} \middle/ \frac{R(k)}{4096^{k}} \right\vert = \frac{1}{4096}. $$
By the ratio test, the series converges absolutely. Since $4096 = 2^{12}$, each additional term contributes exactly 12 bits of precision.
Bellard’s formula converges at 10 bits per term and requires the evaluation of 7 fractions. The above converges at 12 bits per term, and requires the evaluation of 8 fractions. So while we require 20% fewer terms, each term requires about 14% more arithmetic. So, net-net, this formula is approximately 5–6% more efficient.
Proof of identity via a definite integral
Consider $1/(nk+j) = \int_{0}^{1} x^{nk+j-1} dx$. For positive integers $n$ and $b$, we get
$$ \sum_{k=0}^{\infty} \frac{1}{b^k}\cdot\frac{1}{nk+j} = \sum_{k=0}^{\infty} \int_{0}^{1} \left(\frac{x^n}{b}\right)^k x^{j-1} dx. $$
We can swap the sum and integral via the Lebesgue dominated convergence theorem, since the power series $\sum (x^n/b)^k$ converges uniformly for $x \in [0, 1]$ and $b > 1$. Using this and summing the geometric series gives:
$$ \int_{0}^{1} x^{j-1} \sum_{k=0}^{\infty} \left(\frac{x^n}{b}\right)^k dx = \int_{0}^{1} \frac{x^{j-1}}{1 - x^n/b} dx. $$
We now apply this to $S$ termwise with $b=4096=2^{12}$:
$$ S = \int_0^1 \left( \frac{x^{0}}{1 - \frac{x^6}{2^{12}}} - 2^{-5} \frac{x^{2}}{1 - \frac{x^6}{2^{12}}} + 2^{-8} \frac{x^{4}}{1 - \frac{x^6}{2^{12}}} + 2 \frac{x^{0}}{1 - \frac{x^8}{2^{12}}} - 2^{-5} \frac{x^{4}}{1 - \frac{x^8}{2^{12}}} + 2^{-1} \frac{x^{2}}{1 - \frac{x^{12}}{2^{12}}} - 2^{-4} \frac{x^{6}}{1 - \frac{x^{12}}{2^{12}}} - 2^{-8} \frac{x^{10}}{1 - \frac{x^{12}}{2^{12}}} \right) dx. $$
At this point, you could try to algebra your way through, expanding, using the substitution $x=2u$, etc. ultimately yielding a nice denominator $(u^2\pm 2u+2)(u^6-64)(u^{12}-1)$. Maybe compute some residues. Or, just CAS your way through.
% fricas
FriCAS Computer Algebra System
Version: FriCAS 2025.12.23git built with sbcl 2.5.2.1852-1f3beec71
Timestamp: Wed Mar 4 12:41:38 EST 2026
-----------------------------------------------------------------------------
Issue )copyright to view copyright notices.
Issue )summary for a summary of useful system commands.
Issue )quit to leave FriCAS and return to shell.
-----------------------------------------------------------------------------
(1) -> f := (1/(1 - x^6/4096))
- (1/32)*x^2/(1 - x^6/4096)
+ (1/256)*x^4/(1 - x^6/4096)
+ 2*1/(1 - x^8/4096)
- (1/32)*x^4/(1 - x^8/4096)
+ (1/2)*x^2/(1 - x^12/4096)
- (1/16)*x^6/(1 - x^12/4096)
- (1/256)*x^10/(1 - x^12/4096);
Type: Fraction(Polynomial(Fraction(Integer)))
(2) -> normalize(integrate(f, x = 0..1))
3 1 11 19 1
(2) 2 atan(-) - 2 atan(-) + 2 atan(--) + 2 atan(--) + 2 atan(-)
2 2 24 48 4
Type: Expression(Fraction(Integer))
So now we just need to show the arctans all collapse to $\pi$. Recall the identity
$$ \tan^{-1} a \pm \tan^{-1} b = \tan^{-1}\left(\frac{a\pm b}{1\mp ab}\right). $$
The sum of the first four terms can be calculated easily in Common Lisp:
% sbcl --no-inform
* (defun combine (a b) (/ (+ a b) (- 1 (* a b))))
COMBINE
* (reduce #'combine '(3/2 -1/2 11/24 19/48))
4
So we have $2\big(\tan^{-1}4 + \tan^{-1}(1/4)\big)$, and with our final elementary trig identity $\tan^{-1} (a/b) = \pi/2 - \tan^{-1} (b/a)$, we find $S = \pi$.
A new discovery?
Of course, I was excited to find this formula, but after some internet spelunking, it turns out it had already been discovered by Géry Huvent and Boris Gourévitch, perhaps independently. Gourévitch doesn’t credit Huvent as he does with other formulas, but he does say “[…] furthermore, we can obtain BBP formula […] by using what Gery Huvent calls the denomination tables […].” Daisuke Takahashi cites Huvent’s website in this 2019 paper published in The Ramanujan Journal. In all cases, they write the formula in the following way:
$$ \frac{1}{128} \sum _{k=0}^{\infty} \frac{1}{2^{12k}}\left( \frac{768}{24 k+3}+\frac{512}{24k+4}+\frac{128}{24 k+6}-\frac{16}{24 k+12}-\frac{16}{24 k+14}-\frac{12}{24 k+15}+\frac{2}{24 k+20}-\frac{1}{24 k+22}\right), $$
which is structurally equivalent to $S$.
Despite having been known already, this formula doesn’t appear to be well known. As such, I hope this blog post brings more attention to it.
Simple implementation
Here is a simple implementation of digit extraction using BBP-type formulas in Common Lisp:
(defun %pow2-mod (exponent modulus)
(cond
((= modulus 1) 0)
((zerop exponent) 1)
(t
(let ((result 1)
(base (mod 2 modulus))
(e exponent))
(loop :while (plusp e) :do
(when (oddp e)
(setf result (mod (* result base) modulus)))
(setf base (mod (* base base) modulus)
e (ash e -1)))
result))))
(defun %scaled-frac-of-power-two (exponent denom)
(cond
((>= exponent 0)
(let ((residue (%pow2-mod exponent denom)))
(floor (ash residue *precision-bits*) denom)))
(t
(let ((effective-bits (+ *precision-bits* exponent)))
(if (minusp effective-bits)
0
(floor (ash 1 effective-bits) denom))))))
(defun %series-scaled-frac (bit-index bbp-series k-step global-shift alternating-p)
;; A series is a list of series terms. A series term is a quadruple
;; (SIGN SHIFT DENOM-MULTIPLIER DENOM-OFFSET) representing the summand
;; SIGN * 2^SHIFT / (DENOM_MULTIPLIER * k + DENOM_OFFSET).
(let* ((modulus (ash 1 *precision-bits*))
(max-shift (loop :for term :in bbp-series :maximize (second term)))
(k-max (max 0 (ceiling (+ bit-index ; conservative bound
global-shift
max-shift
*precision-bits*
*guard-bits*)
k-step))))
(loop :with acc := 0
:for k :from 0 :to k-max :do
(let ((k-sign (if (and alternating-p (oddp k)) -1 1))
(k-factor (* k-step k)))
(dolist (term bbp-series)
(destructuring-bind (term-sign shift den-mul den-add) term
(let* ((denom (+ den-add (* den-mul k)))
(exponent (+ bit-index global-shift shift (- k-factor)))
(piece (%scaled-frac-of-power-two exponent denom))
(signed (* k-sign term-sign)))
(when (plusp piece)
(setf acc (mod (+ acc (* signed piece)) modulus)))))))
:finally (return acc))))
(defun %nth-hex-from-series (n terms k-step global-shift alternating-p)
(let* ((bit-index (* 4 n)))
(ldb (byte 4 (- *precision-bits* 4))
(%series-scaled-frac bit-index
terms
k-step
global-shift
alternating-p))))
This implementation uses Lisp’s arbitrary precision integer arithmetic. A “real” implementation would use more efficient arithmetic, but this will suffice for some basic testing. Now we can write functions to use the Bellard formula and the new formula:
(defparameter +bellard-terms+
'((-1 5 4 1)
(-1 0 4 3)
(+1 8 10 1)
(-1 6 10 3)
(-1 2 10 5)
(-1 2 10 7)
(+1 0 10 9)))
(defun bellard-nth-hex (n)
(%nth-hex-from-series (* 4 n) +bellard-terms+ 10 -6 t))
(defparameter +new-terms+
'((+1 0 6 1)
(-1 -5 6 3)
(+1 -8 6 5)
(+1 1 8 1)
(-1 -5 8 5)
(+1 -1 12 3)
(-1 -4 12 7)
(-1 -8 12 11)))
(defun new-nth-hex (n)
(%nth-hex-from-series (* 4 n) +new-terms+ 12 0 nil))
Let’s make sure they agree for the first 1000 hex digits:
CL-USER> (loop :for i :below 1000
:always (= (bellard-nth-hex i) (new-nth-hex i)))
T
And now let’s look at timing comparisons. Here’s a little driver:
(defun compare-timings (n)
(flet ((time-it (f n)
(sb-ext:gc :full t)
(let ((start (get-internal-real-time)))
(funcall f n)
(- (get-internal-real-time) start))))
(loop :repeat n
:for index := 1 :then (* 10 index)
:for bellard := (time-it #'bellard-nth-hex index)
:for new := (time-it #'new-nth-hex index)
:do (format t "~v,' D: new is ~A% faster than bellard~%" n index
(round (* 100 (- bellard new)) bellard)))))
And the results if the timing up to the one millionth hexadecimal digit:
CL-USER> (compare-timings 7)
1 : new is 81% faster than bellard
10 : new is 7% faster than bellard
100 : new is 6% faster than bellard
1000 : new is 5% faster than bellard
10000 : new is 4% faster than bellard
100000 : new is 3% faster than bellard
1000000: new is 4% faster than bellard
As predicted, though imperfect a test, it’s consistently faster across a few orders of magnitude.