What is the fastest way to compute digits of $\pi$ using Mathematica?

Question

There are a lot of ways to calculate digits of $\pi$ using Mathematica. The most naïve way I can think of is

N[π, 100000000]

Of course, there are a lot of fast classic formulas (Chudnovsky, Ramanujan) to achieve this goal. I'm wondering what is the fastest way to calculate digits of $\pi$ using Mathematica. The reason that this question may be interesting is that Mathematica has a lot of unique features that can make this calculation faster (or can improve known classic ways of calculating $\pi$).

What are your ideas for calculating digits of $\pi$ using Mathematica in the fastest way possible?

A good answer will involve:

• Explanation of the reason of choosing a particular formula / algorithm.
• Why this particular Method is optimal (in Mathematica at least).
• Optional: Why this Method suits best the use of Mathematica compared with other languages.

---

Note: As @J.M. points, Mathematica implements the Chudnovsky formula for the default calculation.

Answer 1

Here is an adaptation of MATLAB code from Trefethen, Ten Digit Algorithms (2005), based on Borwein & Borwein, The Arithmetic-Geometric Mean and Fast Computation of Elementary Functions (1984) that calculates $\pi$ via the AGM method.

ClearAll[npi];
npi[digits_] := Block[{two, iter},
   iter[{x_, y_, p_}] :=
    With[{s = Sqrt@x},
     {(s + 1/s)/2,
      (y*s + 1/s)/(1 + y),
      p*(1 + x)/(1 + y)}];
   two = SetPrecision[2, 1 + digits];
   With[{y = Sqrt@(Sqrt@two), eps = 10.^(-digits/2)},
    NestWhile[
      iter,
      {(y + 1/y)/2, y, two + Sqrt@2},
      Abs[Last@#1 - Last@#2] > eps &,
      2]
    ]
   ];

Examples:

ClearSystemCache[]    (* clears cached values of Pi *)
digits = 10^6;
N[Pi, digits]; // AbsoluteTiming
pi = Last@npi[digits]; // AbsoluteTiming
pi - Pi // Abs
(*
  {0.393234, Null}
  {6.24918, Null}
  0.*10^-1000000
*)

ClearSystemCache["Numeric"]    (* clears cached values of Pi *)
digits = 7000;
N[Pi, digits]; // AbsoluteTiming
pi = Last@npi[digits]; // AbsoluteTiming
pi - Pi // Abs
(*
  {0.001243, Null}
  {0.008355, Null}
  0.*10^-7000
*)

The AGM algorithm is asymptotically quadratically convergent. Below are the number of digits of accuracy and the ratio with the previous step. It takes 12 iterations to reach 7000 digits and 19 to reach a million digits (in fact, about 1.43 million); in general it will take around Log2[digits] - 1 iterations.

Answer 2

direct implementation of Chudnovsky formula for reference:

a[0] = 1;
a[k_] := a[k] = 
         a[k - 1] (-(((-1 + 2 k) (-5 + 6 k) (-1 + 6 k) 
             (13591409 +545140134 k))/
               (10939058860032000 k^3 (-531548725 + 545140134 k))))
 (pi1 = N[((426880 Sqrt[10005])/(13591409 Sum[ a[k], {k, 0, 500}])),10000]) 
       // AbsoluteTiming // First
 (pi0 = N[Pi, 10000]) // AbsoluteTiming // First

0.15705

0.00156379

 Last@RealDigits[pi0 - pi1]

-7105

Note running it again (using the saved a[i] ) doesn't save that much time:

0.129899

interestingly N[Pi, 10000] gets considerably faster on repeated evals, must be caching something.

Edit: the above is sped up a good bit if we immediately numerically eval each a[k] ( just start with a[0]=N[1,10000] ). With that I get 0.0193068 sec for 7000 digits. Only a factor of 10 off the bulitin..not too bad.

---

here is a version that lets you dial in a specified number of digits:

numdigits = 7000;
(pi1 = N[((426880 Sqrt[
           10005])/(13591409 (NestWhileList[
              Function[{k}, {k, #[[2]] (-(((-1 + 2 k) (-5 + 
                    6 k) (-1 + 6 k) (13591409 + 
                    545140134 k))/
       (10939058860032000 k^3 (-531548725 + 545140134 k))))}][#[[1]] + 1] &,
        {0,N[1, numdigits + 1]}, 
          (-RealDigits[#[[2]]][[2]] < numdigits) &][[All, 2]] // Total))), 
             numdigits + 1];) // AbsoluteTiming // First

0.259337

Last@RealDigits[pi1 - Pi]

7000

It is a good bit slower than the first form though.

Answer 3

Here's another method, based on the fact that $\pi$ is a fixed point of $x + \sin x$ and converges quickly if $x$ starts moderately close to $\pi$. This is clear from the Taylor series or even from this picture:

I don't know that many digits of $\pi$ off the top of my head, but I know enough that one iteration will get me MachinePrecision. From that we can get a (more than) million digits in 11 more iterations. If you know a few more digits of $\pi$, e.g. 3.14159265358979324, then you need only 10 iterations and a 1/4 sec. less time.

ClearSystemCache[]
digits = 10^6;
myPi = Nest[
    Function[x,
     # + Sin@# &@SetPrecision[x, 1 + Min[digits, 3 Precision[x]]]
     ],
    # + Sin@# &@3.14159654, 
    Ceiling@Log[3, digits/16.]
    ]; // AbsoluteTiming
(*  {0.766425, Null}  *)

myPi - Pi
(*  0.*10^-1000001  *)

Some may know this as a calculator trick: Take $x$ to be the calculator's $\pi$ and then find $\sin x$. Add the two numbers by hand to get $\pi$ to more digits than the calculator has, provided the calculator has a good sine routine. (Note: my iPhone calculators's $\pi$ carries extra digits or is a special symbol. The sine of the calculator $\pi$ is zero, and the calculator seems to have quad precision, at least for the significand. The exponent is bounded by two digits, like many calculators.)