10
$\begingroup$

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.

$\endgroup$
  • 2
    $\begingroup$ Note that the implementation notes say that the Chudnovsky formula is used for evaluation. $\endgroup$ – J. M.'s technical difficulties Feb 24 '16 at 18:17
  • $\begingroup$ Thanks! That's important information indeed. But I'm curious if there are other sophisticated ways of perform the calculation. $\endgroup$ – Dargor Feb 24 '16 at 18:18
  • 2
    $\begingroup$ Do you want to pick out single digits or the whole expansion to a certain number of digits? Does it have to be base-10? $\endgroup$ – Kellen Myers Feb 24 '16 at 19:42
  • $\begingroup$ Base 10 would be better than other base. And I want the expansion to certain number of digits. $\endgroup$ – Dargor Feb 24 '16 at 20:17
8
$\begingroup$

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.

Mathematica graphics

| improve this answer | |
$\endgroup$
  • $\begingroup$ I upvoted partly because of your effort to translate one of the ten-digit algorithms. :) $\endgroup$ – J. M.'s technical difficulties Feb 25 '16 at 18:01
  • $\begingroup$ @J.M. Thanks! The Borweins' paper is nice, too. As I recall, there are a couple of judiciously placed typos that help keep you on your toes. :) $\endgroup$ – Michael E2 Feb 25 '16 at 18:14
  • $\begingroup$ I will accept this answer because it's the fastest so far ( 0.008355 vs 0.157986 with 7000 digits ) but both are great answers!. $\endgroup$ – Dargor Mar 1 '16 at 17:39
11
$\begingroup$

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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Maybe it can be compiled for better performance. $\endgroup$ – Dargor Feb 24 '16 at 20:18
  • 3
    $\begingroup$ compile i think only works at machine precision, so you can't simply compile this implementation. $\endgroup$ – george2079 Feb 24 '16 at 22:58
  • $\begingroup$ It would be faster to check the a[k] if they are getting sufficiently small instead of hardcoding an upper limit like 500 so that there are no needless evaluations. $\endgroup$ – J. M.'s technical difficulties Feb 25 '16 at 5:02
  • $\begingroup$ A check would be needed if you want to specify the desired precision. I don't think you'll make it faster if the desired precision is 7105 digits though :) $\endgroup$ – george2079 Feb 25 '16 at 14:32
3
$\begingroup$

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:

enter image description here

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.)

| improve this answer | |
$\endgroup$
  • $\begingroup$ I don't believe I've seen this before. Thanks! $\endgroup$ – Mr.Wizard Jan 10 '19 at 11:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.