03-14-2022, 08:03 PM
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 ?
- 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 )
- 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*N2), and the corrected evaluation is now:
\(\pi\) ~ PN(N) - \(\pi\)/2 * ( 1/N - 1/(4*N2) )
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*N2) )
which, if desired, could be easily converted to the form \(\pi\) ~ PN(N) * ( 1 - 1/(2*N) + 3/(8*N2) + ...) 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:
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.