[VA] SRC #010 - Pi Day 2022 Special

[Valentin Albillo]

                
Today it's March, 14 aka \(\pi\) Day, so Happy \(\pi\) Day to all of you !, and

Welcome to my SRC #010 - \(\pi\) Day 2022 Special, intended to once again commemorate this most famous of constants, \(\pi\).

After posting many threads over the years about \(\pi\), it would seem difficult to find new, interesting appearances of it but actually that's not the case at all, \(\pi\) is inexhaustible and to prove the point let me introduce a new appearance for your enjoyment. Here you are !:

Note: Unlike some previous SRCs, this one is not intended as a challenge to the readers to produce results on their own. Instead, this will be article-like: Here I give outright all the details and my commented results (which were obtained by my own original research, so they've never been available anywhere on the Internet so far), and you can read and enjoy them at once without delay. Come to that, you might try to reproduce and even expand my results using your own calcs if you feel like it. But you saw them here first !

An unexpected infinite product for \(\pi\) and related questions.

Consider this infinite product P(x), for 0 ≤ x < ∞ :

    

Now let me answer the following 4 Questions 4 by first writing a little bit of RPN code for a programmable HP calc, and then using it to do some sleuthing. For speed and accuracy considerations I'll use an ancient version 2.2 (2019) of Free42 without using any of its extended instruction set, so that the code will run unmodified in a physical HP-42S (albeit at reduced speed and accuracy).

First of all, I need to write code to evaluate P(x), and as I can't use an infinite number of terms, I'll store the number of terms to use, N, in variable "N" so that I'll be able to see how it does affect the accuracy of the results obtained, which will prove useful for the sleuthing afterwards.

Thus, this 46-byte program will evaluate P(x) for a given x, assuming that the number N of terms to use has been previously stored in variable "N"

      01  LBL "PX"    10  1          19  STOx 03
      02  STO 02      11  RCL 00     20  DSE 00
      03  RCL "N"     12  X^2        21  DSE 04
      04  STO 00      13  STO 01     22  GTO 00
      05  STO 04      14  1/X        23  RCL 02
      06  1           15  -          24  SQRT
      07  STO- 04     16  RCL 01     25  RCLx 02
      08  STO 03      17  Y^X        26  RCLx 03
      09  LBL 00      18  RCLx 02    27  END

Let's try computing a few values with PX using just N = 10 terms, for speed. We get, for instance:

      FIX 6, 10, STO "N", 1, XEQ "PX" -> 0.000091 ... etc., ...
      _________________________________________________
         N        x = 1         x = 2         x = 3
      ¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­​¯­¯­¯­¯­¯­¯­¯­¯­¯
        10      0.000091      0.131395      9.279782

The Four Questions

a.  Is there a value of x for which P(x) equals \(\pi\) ?

To answer this question I'll first create a wrapper program and then use the [SOLVER] to solve the equation:
    
The wrapper program is this trivial 23-byte piece of code:

      01  LBL "PXEQ"      05  XEQ "PX"
      02  MVAR "N"        06  PI
      03  MVAR "X"        07  -
      04  RCL "X"         08  END

and solving for increasing values of N = 10, 100, 1000, 10000 , we get:

      FIX 8, SOLVER -> Select Solve Program, [PXEQ],
      10, [N][X] -> 2.70596645, 100, [N][X] -> ... etc., ...
      _______________________
         N         x
      ¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­­­¯
        10       2.70596645
        100      2.71814726
        1000     2.71828047
        10000    2.71828181

and it's fairly obvious to any math-inclined person that the root is clearly converging to e = 2.718281828..., the base of the natural logarithms and its inverse, the exponential function, which is thus the answer to the first Question, and so we have the promised unexpected infinite product for \(\pi\) announced in the title, the awesome expression:

    

which beautifully relates \(\pi\) and e. Now, if you remember my last year's SRC #009, I gave there a "trick" expression for \(\pi\) as a function of e, namely:

      \(\pi\) = 4 * ( Arctan e - Arctan \(\frac{e   -   1}{e   +   1}\) )

the catch being that e isn't necessary here at all, an infinity of other values will do, e.g. your age, or your phone number, or your friend's. Can't it be the same case here, that the above infinite product P(x) will evaluate to \(\pi\) for arguments x other than e ? This leads me to the second Question ...

    
b.  We know now that P(e) = \(\pi\). Are there any other such arguments or is e unique ?

In order to answer this question, I'll do some sleuthing after reformulating it as two other related questions, namely:
    
• What's the value of P(x < e) ?
• What's the value of P(x > e) ?
  

Well, using the PX program above for various increasing values of N = 10, 100, 1000, 10000, we get the following, in SCI 4:
      ________________________________________________________________________________​____
         N        x = 1         x = 2         x = 2.7     x = e    x = 2.8      x = 3
      ¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­​¯­¯­¯­¯­¯­¯­¯­¯­¯¯¯¯¯¯­¯­¯­¯­¯­¯¯¯¯¯¯­¯­¯­¯­¯­¯¯¯¯¯­¯­¯­¯­¯­¯¯­¯¯­¯¯­¯
        10      9.0733E-5     1.3140E-1     3.0696E0    3.2950E0  4.4970E0    9.2798E0
        100     7.1239E-44    1.2771E-13    1.6024E0    3.1573E0  6.1958E1    6.3592E4
        1000    9.6769E-435   1.4664E-133   3.6744E-3   3.1432E0  2.3300E13   2.2159E43
        10000   2.1637E-4343  6.1049E-1333  1.5436E-29  3.1417E0  1.3775E129  6.1138E428

and we can clearly see that for arguments x < e the value of P(x) goes to 0, while for arguments x > e it goes to ∞, so the answer to the second question is:

      e  is indeed the only argument which makes this infinite product evaluate to \(\pi\).

Just for fun, if you own some HP calc which has graphics capabilities, try plotting P(x) for x = 0 to 2*e in steps of e/10, for various values of N (say, 10, 100, 1000, ...). Post a screen capture of the plot, if possible. That said, time for third Question ...

                           
c.  Now, fixing x as e, the question is: How many correct digits of \(\pi\) (give or take a few ulps) do we get when using N = 10, 100, 1000, ..., terms ?

To help answer this question, and to speed the computation, now that x is not an argument anymore because it's fixed as e, I've written a new version of PX, a 59-byte program now called PN because it only depends on the number of terms, N, and also optimized for speed by displaying a Wait... message while the program runs (avoids the scrolling symbol) and using stack registers, not as easy to understand as PX but faster:

    01  LBL "PN"    10  E^X        19  X^2         28  X<> ST T
    02  "Wait..."   11  ENTER      20  ENTER       29  ISG ST Y
    03  AVIEW       12  SQRT       21  1/X         30  LBL 00
    04  STO 03      13  RCLx ST Y  22  RCL- 00     31  DSE 02
    05  STO 02      14  STO 01     23  +/-         32  GTO 00
    06  2           15  Rv         24  X<>Y        33  CLST
    07  1           16  LBL 00     25  Y^X         34  CLD
    08  STO 00      17  X<>Y       26  RCLx ST T   35  RCL 01
    09  STO- 02     18  ENTER      27  STOx 01     36  END

Let's use it to obtain these data, in FIX 5 :

    FIX 5, 10, XEQ "PN" -> 3.29501 ... etc., ...
    _______________________
     N             PN
    ¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­­­¯
    10           3.29501
    100          3.15726
    1000         3.14316
    10000        3.14175
    100000       3.14161
    1000000      3.14159

My program benefits from using Free42 Decimal's 34-digit precision, which helps cater for any cumulative rounding errors, and as can be seen PN does indeed converge very slowly to \(\pi\) as N goes to ∞, and its value for N = 1,000,000 comes out as:

    [SHOW] (and hold) -> 3.14159 42243 85727 33446 22511 05879 403

computed accurately to at least 27 digits or better, but giving an estimated value of \(\pi\) accurate to only 7 correct digits, save 2 units in the last place (ulp).

Thus, we have that this infinite product computed to N terms does indeed converge to \(\pi\) but at an excruciatingly slow speed, needing a million terms to get just about 7 digits. And this brings us to the fourth and final Question: Can we do something about it ?

                
d.  Finally, the sleuthing part: Can we do something about the extremely slow convergence ?

Note: I did my own original research so don't search the Internet for any of this 'cause it's not there, this is the first time this info appears on the Internet. It's not rocket science but no one published any of this before.

To try and speed the convergence, first of all I went on to estimate the error, using PN to compute the following values and then subtracting \(\pi\) from each to obtain the errors, which I recognized as being very close to \(\pi\)/2 = 1.57079632679... divided by a million (i.e., N):
    ________________________________________________________________________________​____________________
       N        PN(N)                                       PN(N) - \(\pi\)           Observations
    ¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­­­¯¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯​­¯­¯­¯­¯­¯­­­¯¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­­­¯¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯​­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­¯­­¯­¯­¯­­­¯¯­¯­¯­¯­
    999998    3.14159 42243 88868 93182 82237 27239 246   1.57079907569 E-6   = \(\pi\)/2 * 10.00001750003 E-7
    999999    3.14159 42243 87298 13157 44454 09443 986   1.57079750489 E-6   = \(\pi\)/2 * 10.00000750001 E-7
    1000000   3.14159 42243 85727 33446 22511 05879 403   1.57079593410 E-6   = \(\pi\)/2 *  9.99999750000 E-7
    1000001   3.14159 42243 84156 54049 16628 07700 354   1.57079436330 E-6   = \(\pi\)/2 *  9.99998750002 E-7
    1000002   3.14159 42243 82585 74966 25768 41450 776   1.57079279251 E-6   = \(\pi\)/2 *  9.99997750005 E-7

thus, all the errors seem to be very close to \(\pi\)/2 * 1E-6, with the smallest one occurring for N = 1000000, where the error is \(\pi\)/2 * 9.99999750000 E-7. Noticing this, I then used \(\pi\)/2 * 1E-6 as a correction term to be applied to the computed PN(1000000), obtaining the following:

    PN(1000000) - \(\pi\)/2 * 1E-6 = 3.14159 26535 89400 53956 56318 74557 711

and subtracting \(\pi\), the absolute error now is ~ 3.92699 E-13, which means we've got about 14 correct digits (save 4 ulp), where previously we had just 7 correct digits (save 2 ulp). In other words, applying this extremely simple correction term, \(\pi\)/2 * 1/N, essentially duplicates the number of correct digits.

Can we do better ?  Yes, we can. Observing the errors using just PN(N) above for N = 999998 to N = 1000002, we notice that not only are they of the form

    \(\pi\)/2 * 1E-6 ~ \(\pi\)/2 * 1/N ,

but the actual differences with respect to that value also have a very regular form:

    ..750003, ..750001, ..750000, ..750002, ..750005,

which suggests a *second* correction term to cater for the ..75 difference. To cut to the chase, after a few trivial arithmetic operations the second correction term is immediately found to be in absolute value equal to \(\pi\)/2 * 1/(4*N²), and the corrected evaluation is now:

    \(\pi\) ~ PN(N) - \(\pi\)/2 * ( 1/N - 1/(4*N²) )

and as \(\pi\) appears on both the LHS and the RHS, we proceed to isolate \(\pi\) at the LHS, which gives:

    \(\pi\) ~ PN(N) / ( 1 + 1/(2*N) - 1/(8*N²) )

which, if desired, could be easily converted to the form  \(\pi\) ~ PN(N) * ( 1 - 1/(2*N) + 3/(8*N²) + ...) by polynomial division, but the above expression will do for now, as we do not have enough additional terms to do an accurate polynomial division anyway.

This short additional code applies both correction terms to the output of PN(N). First, change 36 END to 36 STOP and then include after it the following lines:

      37  RCL 03       41  X^2      45  -     49  END
      38  RCL+ 03      42  8        46  1
      39  1/X          43  x        47  +
      40  RCL 03       44  1/X      48  /

This adds just 17 bytes to PN and executes instantly but as we'll see in a moment, it greatly increases the number of correct digits. To use it, simply:

      N (number of terms), XEQ "PN" -> (shows computed PN(N) and pauses), R/S -> (shows corrected value)

When particularized for N = 1,000,000, the corrected evaluation gives:

         PN(N) = 3.14159 42243 85727 33446 22511 05879 403 ( 7 correct digits save 2 ulp )
     Corrected = 3.14159 26535 89793 23864 73305...
             \(\pi\) = 3.14159 26535 89793 23846 26433...
         Error ~ 1.84687 E-19 ( i.e. 20 correct digits save ~ 2 ulp )

This means we've got essentially 20 correct digits using just two simple, inexpensive correction terms, while the original uncorrected PN(1000000) gave us only about 7 correct digits. Let's check the results for other values of N, for instance:
         
• For N = 100,000:
  
           PN(N) = 3.14160 83615 13791 56287 28512 11516 805 ( 6 correct digits save 2 ulp )
       Corrected = 3.14159 26535 89793 52207...
               \(\pi\) = 3.14159 26535 89793 23846...
           Error ~ 2.83612 E-16 ( 17 correct digits save 2 ulp )
   
• For N = 10,000:
  
           PN(N) = 3.14174 97292 95765 50614 37729 49086 661 ( 5 correct digits save 2 ulp )
       Corrected = 3.14159 26535 90076 819...
               \(\pi\) = 3.14159 26535 89793 238...
           Error ~ 2.83581 E-13 ( 14 correct digits save 3 ulp )

And it seems that using the two correction terms we've obtained, we empirically have:

    New #correct digits = 3 * Old #correct digits - 1

Thus, while using just one correction term duplicates the number of correct digits, using two correction terms essentially triples the precision obtained, i.e. :

    N = 10,000     (5 correct digits)  -> 3 * 5 - 1 = 14 correct digits
    N = 100,000    (6 correct digits)  -> 3 * 6 - 1 = 17 correct digits
    N = 1,000,000  (7 correct digits)  -> 3 * 7 - 1 = 20 correct digits

and of course the number of correct digits could be increased even further by simply obtaining additional correction terms, either empirically as I've done here or better yet, analytically.

Matter of fact, I've managed to obtain two additional terms empirically, but giving details here would make this already humongous exposition even much longer, so that's left as an exercise for the interested reader.

That's all. Any and all constructive and on-topic comments are most welcome and appreciated.
  
V.