A Fast Self-correcting π Algorithm

by   Tsz-Wo Sze, et al.
The Apache Software Foundation

We have rediscovered a simple algorithm to compute the mathematical constant π=3.14159265⋯. The algorithm had been known for a long time but it might not be recognized as a fast, practical algorithm. The time complexity of it can be proved to be O(M(n)log^2 n) bit operations for computing π with error O(2^-n), where M(n) is the time complexity to multiply two n-bit integers. We conjecture that the algorithm actually runs in O(M(n)log n). The algorithm is self-correcting in the sense that, given an approximated value of π as an input, it can compute a more accurate approximation of π with cubic convergence.


page 1

page 2

page 3

page 4


A faster algorithm for Cops and Robbers

We present an algorithm of time complexity O(kn^k+2) deciding whether a ...

Coloring Fast Without Learning Your Neighbors' Colors

We give an improved randomized CONGEST algorithm for distance-2 coloring...

On the number of error correcting codes

We show that for a fixed q, the number of q-ary t-error correcting codes...

On Euclidean Methods for Cubic and Quartic Jacobi Symbols

We study the bit complexity of two methods, related to the Euclidean alg...

Foundation for a series of efficient simulation algorithms

Compute the coarsest simulation preorder included in an initial preorder...

A Fast Tree Algorithm for Electric Field Calculation in Electrical Discharge Simulations

The simulation of electrical discharges has been attracting a great deal...

RD-Optimized Trit-Plane Coding of Deep Compressed Image Latent Tensors

DPICT is the first learning-based image codec supporting fine granular s...

1 Introduction

The computation of the mathematical constant has drawn a great attention from mathematicians and computer scientists over the centuries [3, 14]. The known asymptotically fastest algorithms for computing run in

bit operations with error , where is the time complexity to multiply two -bit integers. The AGM algorithms [6, 15, 4] are the only examples before this paper. If the recent result in [10] is correct,

The known asymptotically fastest algorithms run in

The Chudnovsky algorithm [8], which runs in

is a popular implementation choice. The computer program, y-cruncher, implemented the Chudnovsky algorithm has been used to compute to 31.4 trillion digits [18, 11].

In this paper, we present an algorithm. From the best of our knowledge, the algorithm is new although the proofs are simple. It is self-correcting in the sense that, given an approximated value of as an input, it can compute a more accurate value of with cubic convergence.

Our algorithm and the AGM algorithms share the same time complexity. Similar to the Chudnovsky algorithm, our algorithm uses binary splitting. According to [7, 19], binary splitting has advantages over AGM including:

  1. the implicit constants for binary splitting are smaller than the ones for AGM;

  2. binary splitting can be speeded up by simultaneously summing up many terms at once but it is difficult to speed up AGM;

  3. AGM has very poor memory locality.

Moreover, the AGM iteration is not self-correcting so that full precision is required throughout. In contrast, the intermediate results can be truncated in our algorithm. For example, suppose the current step has computed in decimal places and the next target precision is decimal places for some . Then the current result can be truncated to . Thus, our algorithm potentially runs faster than the AGM algorithms in practice. We have the following hypothesis.

Hypothesis 1.

Algorithm 3 is faster than the AGM algorithms by a constant factor.

Our algorithm is presented in the next section. We discuss the verification problem in Section 3. Finally, we show a family of sequences converging to in Section 4.

2 The Computational Problem

Let be an approximated value of and


be the error with


for some fixed . By the Taylor series


it is easy to see that


Note that


Finally, we obtain a better approximated value of


such that the error

becomes cubic by (2.1), (2.2), (2.4), (2.5) and (2.6). We have proved the following theorem.

Theorem 2 (Cubic Convergence).

Let be an approximated value of such that , Then, , where .

We present our algorithm and then prove its time complexity.

Algorithm 3 (Self-correcting Computation).

The input is a positive integer . This algorithm returns such that .

  1. Let and .

  2. For , use -bit precision to compute

  3. Return .

Theorem 4.

Algorithm 3 runs in bit operations.


The main ingredient of Algorithm 3 is to compute . Consider as a rational number since it is a truncated value of . For a -bit rational number with , both and can be computed using the binary splitting sine algorithm [12, 13]. Therefore, and can be computed in . Then use the doubling formulae

to compute ), ) and, finally, . The time complexity to compute is . The time complexity of Algorithm 3 is . ∎

Note that the AGM sine algorithm [5], which runs in , cannot be used here since it requires as an input. Note also that the binary splitting algorithm can be used to compute directly [12]. However, the time complexity is .

2.1 A Numerical Example

The following example has been computed by PARI/GP [17] and GMP [9]. We simply have used the sine function provided by PARI/GP. In the table below, is the approximated value of in iteration , is an upper bound of the error and is the precision in , where

0 3 3 0.141120008059867222100744802808110
1 3.141 10 0.000592653555099468066916718249636
2 3.1415926535 30 0.000000000089793238462643383279382

We have the following sequence converging to ,

3 The Decision Problem

Let with decimal places be a computed value of . How to verify if the digits are correct? In other words, verify if


It is interesting to ask if the decision problem, i.e. checking (3.1) for a given in decimal places, is easier than the computational problem, i.e. computing in decimal places. An algorithm deciding (3.1) asymptotically faster than computing has not been discovered.

The self-correcting step in Algorithm 3 can be used for verification. Split

into higher digits and lower digits for some such that

is expected to have a few more correct digits than . Check if all the digits in match . The time complexity is . If Hypothesis 1 is correct, the decision problem is slightly easier than the computational problem by a constant factor since the decision problem takes only the last iteration but the computational problem needs iterations.

In practice, after is computed in decimal places by an algorithm, a different algorithm or the same algorithm with a different set of parameters is used to verify the result.

The result mentioned in the introduction has decimal digits111Note that . and 26,090,362,246,629 hexadecimal digits [18, 11]. The computation used the Chudnovsky algorithm. For verification, the Bailey–Borwein–Plouffe (BBP) formula [1] and also the Bellard’s improved BBP formula [2] were used to compute 48 hexadecimal digits starting at the 26,090,362,246,601st position. There were 29 hexadecimal digits,

agreed in all three results from Chudnovsky, BBP and Bellard.

In 2010, we computed the two quadrillionth bit of [16] using Bellard’s formula. Two computations at two different bit positions,

were executed. There were 256 bits agreed in both computations.

4 Convergent Sequences

We extend Theorem 2 to show a family of sequences converging to below.

Lemma 5.

For , define

where . Then,


We will show


Then Theorem 2 implies

Now we show (4.1). If , we are done. Assume . There exists the least integer such that . If not, let be the least upper bound of . Let

Since is bounded above by , we have and for all . Since is increasing with the least upper bound ,


However, (4.2) is contradiction since if , is increasing; otherwise, for large enough if . Therefore,

Since is the least integer, . If , we have ; otherwise, and . In both cases,

For any with , define

if and only if

for some integer . We show a more general theorem below.

Theorem 6 (Convergent Sequences).

For any such that

For , define



If , it is trivial. Assume .

Let and

We have . For , define

It is obvious that, for ,

If , Lemma 5 implies . Then

Suppose . Let so that . For , define

Lemma 5 implies . Since for all ,


  • [1] David Bailey, Peter Borwein, and Simon Plouffe. On the rapid computation of various polylogarithmic constants. Mathematics of Computation, 66(216):903–913, apr 1997.
  • [2] Fabrice Bellard. A new formula to compute the n’th binary digit of pi. http://bellard.org/pi/pi_bin.pdf, 1997.
  • [3] Lennart Berggren, Jonathan Borwein, and Peter Borwein. Pi: A Source Book. Springer, New York, NY, USA, 3rd edition, 2004.
  • [4] Jonathan Borwein and Peter Borwein. Pi and the AGM: A Study in Analytic Number Theory and Computational Complexity. John Wiley & Sons, Toronto, Canada, 3rd edition, 1987.
  • [5] Richard Brent. Fast multiple-precision evaluation of elementary functions. J. ACM, 23(2):242–251, April 1976.
  • [6] Richard Brent. Multiple-precision zero-finding methods and the complexity of elementary function evaluation. Analytic Computational Complexity, pages 151–176, 1976.
  • [7] Richard Brent and Paul Zimmermann. Modern Computer Arithmetic. Cambridge University Press, New York, NY, USA, 2010.
  • [8] David Chudnovsky and Gregory Chudnovsky. The computation of classical constants. Proceedings of the National Academy of Science USA, 86:8178–8182, 1989.
  • [9] Torbjörn Granlund and the GMP development team. GNU MP: The GNU Multiple Precision Arithmetic Library, 6.1.2 edition, 2016. http://gmplib.org/.
  • [10] David Harvey and Joris Van Der Hoeven. Integer multiplication in time O(n log n). Preprint, March 2019.
  • [11] Emma Haruka Iwao. Pi in the sky: Calculating a record-breaking 31.4 trillion digits of archimedes’ constant on google cloud. https://cloud.google.com/blog/products/compute/calculating-31-4-trillion-digits-of-archimedes-constant-on-google-cloud, 2019.
  • [12] Ekatherina Karatsuba. Fast evaluation of transcendental functions. Probl. Peredachi Inform, 27:76–99, 1991.
  • [13] Ekatherina Karatsuba. Fast Computation of Some Special Integrals of Mathematical Physics, pages 29–40. Springer US, Boston, MA, 2001.
  • [14] Donald Knuth. The Art of Computer Programming, Volume 2: Seminumerical Algorithms. Addison-Wesley, Reading, MA, USA, 3rd edition, 1997.
  • [15] Eugene Salamin. Computation of

    using arithmetic-geometric mean.

    Mathematics of Computation, 30(135):565–570, jul 1976.
  • [16] Tsz-Wo Sze. The two quadrillionth bit of pi is 0! distributed computation of pi with apache hadoop. Proceedings - 2nd IEEE International Conference on Cloud Computing Technology and Science, CloudCom 2010, pages 727 – 732, 01 2011.
  • [17] The PARI Group, Univ. Bordeaux. PARI/GP version 2.11.2, 2018. available from http://pari.math.u-bordeaux.fr/.
  • [18] Alexander Yee. Google cloud topples the pi record. http://www.numberworld.org/blogs/2019_3_14_pi_record/, 2019.
  • [19] Alexander Yee. y-cruncher - frequently asked questions. http://numberworld.org/y-cruncher/faq.html, 2019.