[VA] SRC#006- Pi Day 2020 Special: A New Fast Way to Compute Pi

[Valentin Albillo]

 
Hi all, welcome to my SRC#006- Pi Day 2020 Special:

To commemorate Pi Day 2020 I'll give and comment here a new, extremely fast way to compute Pi to arbitrary precision, millions or billions of digits if desired, but first a little background of some of my recent posts on the subject matter.
          
In the past I've posted here a couple of novel ways to compute Pi, quite unexpected but very slow. The first one, statistical in nature, was (HP-71B BASIC code):

      1   DESTROY ALL @ RANDOMIZE 1 @ FOR K=1 TO 5 @ N=10^K @ S=0
      2   FOR I=1 TO N @ IF NOT MOD(IROUND(RND/RND),2) THEN S=S+1
      3   NEXT I @ P=S/N @ STD @ DISP N, @ FIX 3 @ DISP 5-P*4 @ NEXT K

>RUN

      10          3.000
      100         3.360
      1000        3.076
      10000       3.125
      100000      3.142

where the even/odd nature of the random-valued expression IROUND(RND/RND) is the key, no circles in sight at all despite what many believe to be always the case.

The second one (posted originally here) was also incredibly simple and this was my implementation for the HP-71B as a user-defined function:

      1   DEF FNP(N)
      2       T = N
      3       FOR I = N-1 TO 2 STEP -1
      4             T = CEIL(T/I)*I
      5       NEXT I
      6       FNP = N*N/T
      7   END DEF


These are the function's values for N = 10, 100, ..., 1000000:

      DESTROY ALL @ FIX 5

      FNP(10)          -> 2.94118
      FNP(100)         -> 3.09215
      FNP(1000)        -> 3.13903
      FNP(10000)       -> 3.14133
      FNP(100000)      -> 3.14153
      FNP(1000000)     -> 3.14159

The limit for  N -> Inf   is exactly Pi, of course.

As stated above, these methods are quite unexpected, very rarely seen (if at all), and interesting from a purely academic point of view but utterly impractical for computing Pi as they are unbearably slow, producing just 5-6 digits of Pi after millions of iterations. Even producing as few as 10 digits is probably out of the question.

On the other hand, the method I'll introduce here is equally simple, perhaps even simpler in a way, but capable of producing millions, billions and even trillions of digits of Pi as fast or faster than other well-known methods, assuming of course the necessary hardware resources (fast processors, enough RAM) and multiprecision software are available.

For starters, here is the HP42S version, just 7 steps (12 bytes) of RPN code, particularized to carry out exactly 3 iterations of the algorithm and requiring no user input at all (but make sure to set RAD mode).

      01    3
      02  ENTER
      03  LBL 00
      04    SIN
      05    RCL+ ST L
      06    DSE ST Y
      07  GTO 00

When run in Free42, which allows for 34-digit precision, it will instantly produce this 34-digit value of Pi:

      [R/S] -> 3.14159265359
      [SHOW]-> 3.141592653589793238462643383279502

which is the correct value of Pi when truncated to 34 digits. Should Free42 support more than 34 digits, say 100, 1000 or a million, this same code (trivially modified to perform more than 3 iterations) would produce the value of Pi to the maximum precision supported.

The algorithm used is extremely simple, and it consists of the following steps:

      Step 1: Assign to x some approximation to Pi, say x = 3 will do.

      Step 2: Assign x += sin(x) . This new value will have at least 3 times as many correct digits of Pi as the previous value did.

      Step 3: Repeat Step 2 above until the desired precision is achieved (the number of iterations needed is easily computed in advance, see below).


Now you may be thinking something like: "Hey, you're using a trigonometric function, namely sin(x), to help compute the value of Pi and that's circular [pun intended]. Come to that, you might simply compute Pi = 4*arctan(1), say."

My reply is: first of all, arctan is an inverse function, namely the inverse of tan(x) which, as the inverses of sin(x) and cos(x) do, can produce Pi or rational fractions of Pi (say Pi/2, Pi/4, Pi/3, etc) for rational or algebraic values of their arguments. But sin(x) is not an inverse function like those and won't produce Pi or any fractions of Pi for any real rational or algebraic values of its argument. In other words, there's no rational or algebraic value x for which sin(x) = Pi/4 or Pi/6 or 2/3*Pi, etc.

Secondly, what matters here is not which functions we do use but how fast can we implement them and how cleverly we can use them to compute Pi as fast as possible, so to be able to use this algorithm to compute Pi to millions or billions of digits, apart from suitable hardware we only need a fast sin(x) custom implementation, which only needs to be able to compute sin(x) for x in the very short range [3 .. Pi], say. Pi itself isn't needed in the implementation as no argument reduction is ever performed, which a general implementation of sin(x) would require. Further, our custom implementation can use a number of arithmetic tricks to reduce the argument to a suitably small value which might greatly speed the computation, say using the double-angle (half-angle) formulas three times in succession to reduce the [3 .. Pi] argument range to the still smaller [3/8 .. Pi/8] range.

The question is then: how fast can sin(x) be computed ? In the paper "Fast Multiple-Precision Evaluation of Elementary Functions" by Richard P. Brent (Australian National University, Canberra, Australia) it is shown that sin(x) (and some other elementary functions as well) can be evaluated with relative error 0(2^-n), in O(M(n)log(n)) operations as n tends to infinity, for any floating-point argument x in a suitable finite interval, where M(n) is the number of single-precision operations required to multiply n-bit integers.

Thus, the main advantages of this algorithm can be summarized like this:

1) It's extremely simple, each iteration requires the computation of just one sin(x) function and one in-place addition, with no other operations or variables involved. Futher, the starting value for the very first iteration is simply the constant 3, no need to compute many-million-digit values for the starting parameters as other methods require.

2) It's extremely fast, the sin(x) can be performed in O(M(n)log(n)) operations, which is about as fast as it gets, and it can be optimized to work for just a very short range of x.

3) It converges cubically, i.e., each iteration provides a value of Pi which has at least 3 times as many correct digits as the previous iteration, so the number of digits obtained grows exponentially and thus computing millions, billions or trillions of digits requires very few iterations, namely:

      13 iterations provide ~ 2 million digits
      19 iterations provide ~ 1.4 billion (1e9) digits
      25 iterations provide ~ 1 trillion (1e12) digits

4) It is self-correcting, if you're going to compute a million digits of Pi you don't need to use one-million-digit operations from the very first iteration. You can specify the level of precision increasingly for each iteration, namely 5 digit, 13-digit, 35-digit, 102-digit, 303-digit and 905-digit precision for iterations 1..6, respectively, which obviously increases speed enormously.

How does it compare with the best-known methods ? In the paper "The Life of Pi: From Archimedes to Eniac and Beyond" of Jonathan M. Borwein, Frsc, Faa, we find this quartic method described:




Being quartic, this method requires less iterations to achieve a given number of digits, namely 20 iterations to provide a trillion (1e12) digits vs. 25 iterations for the cubic method described here, 25 iterations to provide a quadrillion (1 e15) digits vs. 31, and 30 iterations to provide a quintillion (1e18 digits) vs. 37. As can be seen, there's no large increase in the number of iterations for practically feasible computations even as large as a quintillion digits so the speed advantage of being a quartic method instead of a cubic one is not crucial.

Also, the method here requires just 1 sin and 1 addition per iteration and just one variable, x. On the other hand, Borwein's method requires about 2 raising-to-the-4th power operations, 3 multiplications, 1 squaring, 3 subtractions, 4 additions, 1 division, 1 fourth-root, and at least storing 3 intermediate values, all of it per iteration. Further all operations must be performed at full precision from the very first one, plus computing the two irrational initial values to full precision as well before even starting.

This means that, in practice and for up to at least a quadrillion/quintillion digits, the method described here could be significantly faster if a suitably optimized, fast custom version of sin(x) is used. Now let's see the method in action by conducting an interactive session using a multiprecision environment. We'll use the environment's native built-in sin(x) as we aren't interested here in speed but in showing the accuracy obtained for an increasing number of iterations:

> X=3

3 (starting value, 1 correct digit)

> X+=sin(X)

3.141120008... (1st iteration, 4 correct digits)

> X+=sin(X)

3.1415926535721955... (2nd iteration, 11 correct digits)

> X+=sin(X)

3.1415926535897932384626433832795019759... (3rd iteration, 33 correct digits)

> X+=sin(X)

3.141592653589793238462643383279502884197169399375105820974944592307816406286208​
99862803482534211706785726... (4th iteration, 100 correct digits)

> X+=sin(X)

3.141592653589793238462643383279502884197169399375105820974944592307816406286208​99862803482534211706798214808651328230664709384460955058223172535940812848111745​02841027019385211055596446229489549303819644288109756659334461284756482337867831​
6527120190914564856692346034861045432664821339360726024914127340000...
(5th iteration, 301 correct digits)

> X+=sin(X)

3.141592653589793238462643383279502884197169399375105820974944592307816406286208​99862803482534211706798214808651328230664709384460955058223172535940812848111745​02841027019385211055596446229489549303819644288109756659334461284756482337867831​65271201909145648566923460348610454326648213393607260249141273724587006606315588​17488152092096282925409171536436789259036001133053054882046652138414695194151160​94330572703657595919530921861173819326117931051185480744623799627495673518857527​24891227938183011949129833673362440656643086021394946395224737190702179860943702​77053921717629317675238467481846766940513200056812714526356082778577134275778960​91736371787214684409012249534301465495853710507922796892589235420199561121290219​60864034418159813629774771309960518707211349999998372978049951059731732816096318​59502445945534690830264252230825334468503526193118817101000313783875288658753320​
8381420617177669147303592553972... (6th iteration, 903 correct digits)


and so on and so forth.


Happy  Pi  Day !

V.