Threaded Mode | Linear Mode

Gerson W. Barbosa · (This post was last modified: 01-04-2021 11:53 PM by Gerson W. Barbosa.)

The alternate sum of the factors of the Wallis product tends to \(\pi\)/4 - 1/2 as \(n\) tends to infinity:

\(\lim_{n\rightarrow \infty } \left [ \frac{1}{1\cdot 3}-\frac{1}{3\cdot 5}+\frac{1}{5\cdot 7}-\frac{1}{7\cdot 9}+\frac{1}{9\cdot 11}-\frac{1}{11\cdot 13}+\frac{1}{13\cdot 15}-\frac{1}{15 \cdot 17}\pm \cdots +\frac{(-1)^{n-1}}{4n^{2}-1}\right ]=\frac{\pi }{4}-\frac{1}{2}\)

This converges to the constant significantly faster than the Wallis Product. For 12 correct digits of \(\pi\) it requires \(\sqrt{10^{12}}\) terms instead of the \(10^{12}\) terms required by
Wallis. That's one million times as fast in this case, but obviously not enough for increasing numbers of digits. Fortunately, again, a correction continued fraction is possible which provides linear convergence, about 25/12 digits per term as the other formula.

\(\pi = \lim_{n\rightarrow \infty } \left [4\left ( \frac{1}{2} + \sum_{k=1}^{n}\frac{(-1)^{k-1}}{4k^{2}-1} \right )+\frac{(-1)^{n}}{2n(n+1)+\frac{3}{2}-\frac{3}{2n(n+1)+\frac{19}{2}-\frac{60}{2n(n+1)+\frac{51}{2}-\frac{315}{\frac{\ddots }{\cdots 2n(n+1)+\frac{8n(n-1)+3}{2}-\frac{n^{2}(4n^{2}-1)}{2n(n+1)+\frac{8n(n+1)+3}{2}}} }}}}\right ]\)

Tentatively, I think I can name this one as "Wasicki's formula" :-) Well, at least while an older equivalent version doesn´t appear around.

The algorithm is very simple, as tested on the HP-75C:

10 INPUT N
15 C=0
20 D=6*N*(N+1)+3/2
25 E=8*N
30 M=4*N*N-1
35 S=1
40 W=0
45 FOR K=1 TO N
50 W=W-S/M
55 C=-M*(M+1)/(4*(C+D))
60 M=M-E+4
65 D=D-E
70 E=E-8
75 S=-S
80 NEXT K
85 DISP S*(4*W+1/(C+D))+2
>run
?5
3.14159265359

I can get one thousand digits in little more than one fourth of a second on my old desktop computer, using that algorithm:

n =  480

3.1415926535 8979323846 2643383279 5028841971 6939937510  (0050)
  5820974944 5923078164 0628620899 8628034825 3421170679  (0100)
  8214808651 3282306647 0938446095 5058223172 5359408128  (0150)
  4811174502 8410270193 8521105559 6446229489 5493038196  (0200)
  4428810975 6659334461 2847564823 3786783165 2712019091  (0250)
  4564856692 3460348610 4543266482 1339360726 0249141273  (0300)
  7245870066 0631558817 4881520920 9628292540 9171536436  (0350)
  7892590360 0113305305 4882046652 1384146951 9415116094  (0400)
  3305727036 5759591953 0921861173 8193261179 3105118548  (0450)
  0744623799 6274956735 1885752724 8912279381 8301194912  (0500)
  9833673362 4406566430 8602139494 6395224737 1907021798  (0550)
  6094370277 0539217176 2931767523 8467481846 7669405132  (0600)
  0005681271 4526356082 7785771342 7577896091 7363717872  (0650)
  1468440901 2249534301 4654958537 1050792279 6892589235  (0700)
  4201995611 2129021960 8640344181 5981362977 4771309960  (0750)
  5187072113 4999999837 2978049951 0597317328 1609631859  (0800)
  5024459455 3469083026 4252230825 3344685035 2619311881  (0850)
  7101000313 7838752886 5875332083 8142061717 7669147303  (0900)
  5982534904 2875546873 1159562863 8823537875 9375195778  (0950)
  1857780532 1712268066 1300192787 6611195909 2164201989  (1000)

Runtime: 0.27 seconds

Edited to correct a mistake, per Valentín Albillo’s note below.

Valentin Albillo · 01-04-2021, 11:21 PM

(01-04-2021 08:41 PM)Gerson W. Barbosa Wrote: This converges to the constant significantly faster than the Wallis Product. For 12 correct digits of \(\pi\) it requires \(\sqrt{10^{12}}\) terms instead of the \(10^{12}\) terms required by
Wallis. That's 1000 times as fast, [...]

Nope. The square root evaluates to 10⁶ so it's 1,000,000 times as fast, not 1000 as you say.

Regards.
V.

Gerson W. Barbosa · 01-04-2021, 11:36 PM

(01-04-2021 11:21 PM)Valentin Albillo Wrote:
(01-04-2021 08:41 PM)Gerson W. Barbosa Wrote: This converges to the constant significantly faster than the Wallis Product. For 12 correct digits of \(\pi\) it requires \(\sqrt{10^{12}}\) terms instead of the \(10^{12}\) terms required by
Wallis. That's 1000 times as fast, [...]

Nope. The square root evaluates to 10⁶ so it's 1,000,000 times as fast, not 1000 as you say.

Regards.
V.

You are right, thanks! I was thinking of the n=1000 example I had just checked on my smartphone when I wrote that. I’d like to have checked the sum of the first 1000000 terms on Free42, but my program would fail short of n = 100000.
Hopefully the rest is correct.

Best regards,

Gerson.

Albert Chan · 01-05-2021, 10:50 PM

(01-04-2021 08:41 PM)Gerson W. Barbosa Wrote: The alternate sum of the factors of the Wallis product tends to \(\pi\)/4 - 1/2 as \(n\) tends to infinity:

\(\lim_{n\rightarrow \infty } \left [ \frac{1}{1\cdot 3}-\frac{1}{3\cdot 5}+\frac{1}{5\cdot 7}-\frac{1}{7\cdot 9}+\frac{1}{9\cdot 11}-\frac{1}{11\cdot 13}+\frac{1}{13\cdot 15}-\frac{1}{15 \cdot 17}\pm \cdots +\frac{(-1)^{n-1}}{4n^{2}-1}\right ]=\frac{\pi }{4}-\frac{1}{2}\)

This converges to the constant significantly faster than the Wallis Product.
For 12 correct digits of \(\pi\) it requires \(\sqrt{10^{12}}\) terms instead of the \(10^{12}\) terms required by Wallis.

Where does this limit comes from ?

CAS can get pi/4 - 1/2, but the "proof" felt cheating.

Cas> expand(sum(1/((4*k-3)*(4*k-1))-1/((4*k-1)*(4*k+1)),k = (1 .. ∞))) → 1/4*π-1/2

---

I find this remarkable formula for approximating the sum of alternating series.

Example: above pi formula, for 1000 terms

To reduce errors, do the sum in pairs, and backwards.

lua> s = 0
lua> for k=500,1,-1 do c=4*k-1; s=s+4/(c*(c*c-4)) end
lua> (s + 0.5) * 4
3.141592154089668

Add correction, from next 3 terms:

lua> function term(k) return (-1)^(k-1)/(4*k*k-1) end
lua> f1, f2, f3 = term(1001), term(1002), term(1003)
lua> s = s + (7*f1 - f2 - 2*f3)/12 -- alternate series correction

lua> (s + 0.5) * 4
3.141592653588056
lua> pi
3.141592653589793

Gerson W. Barbosa · 01-06-2021, 12:53 AM

(01-05-2021 10:50 PM)Albert Chan Wrote:
(01-04-2021 08:41 PM)Gerson W. Barbosa Wrote: The alternate sum of the factors of the Wallis product tends to \(\pi\)/4 - 1/2 as \(n\) tends to infinity:

\(\lim_{n\rightarrow \infty } \left [ \frac{1}{1\cdot 3}-\frac{1}{3\cdot 5}+\frac{1}{5\cdot 7}-\frac{1}{7\cdot 9}+\frac{1}{9\cdot 11}-\frac{1}{11\cdot 13}+\frac{1}{13\cdot 15}-\frac{1}{15 \cdot 17}\pm \cdots +\frac{(-1)^{n-1}}{4n^{2}-1}\right ]=\frac{\pi }{4}-\frac{1}{2}\)

This converges to the constant significantly faster than the Wallis Product.
For 12 correct digits of \(\pi\) it requires \(\sqrt{10^{12}}\) terms instead of the \(10^{12}\) terms required by Wallis.

Where does this limit comes from ?

Sum(n=1, inf, (-1)^(n-1)/(4n^2-1)) (WolframAlpha)

Albert Chan · 01-06-2021, 01:32 AM

Quote:\(\lim_{n\rightarrow \infty } \left [ \frac{1}{1\cdot 3}-\frac{1}{3\cdot 5}+\frac{1}{5\cdot 7}-\frac{1}{7\cdot 9}+\frac{1}{9\cdot 11}-\frac{1}{11\cdot 13}+\frac{1}{13\cdot 15}-\frac{1}{15 \cdot 17}\pm \cdots +\frac{(-1)^{n-1}}{4n^{2}-1}\right ]=\frac{\pi }{4}-\frac{1}{2}\)

Proof:

atan(x) = x - x^3/3 + x^5/5 - x^7/7 + ...
atan(1) = pi/4 = 1 - 1/3 + 1/5 - 1/7 + ...

2*I = 2/(1*3) - 2/(3*5) + 2/(5*7) - 2/(7*9) + ...
= (1-1/3) - (1/3-1/5) + (1/5-1/7) - (1/7-1/9) + ...
= (1-1/3+1/5-1/7+...) + (-1/3+1/5-1/7+...)
= pi/4 + (pi/4 - 1)

→ I = pi/4 - 1/2

toml_12953 · 01-06-2021, 02:10 AM

(01-04-2021 08:41 PM)Gerson W. Barbosa Wrote: The algorithm is very simple, as tested on the HP-75C:

10 INPUT N
15 C=0
20 D=6*N*(N+1)+3/2
25 E=8*N
30 M=4*N*N-1
35 S=1
40 W=0
45 FOR K=1 TO N
50 W=W-S/M
55 C=-M*(M+1)/(4*(C+D))
60 M=M-E+4
65 D=D-E
70 E=E-8
75 S=-S
80 NEXT K
85 DISP S*(4*W+1/(C+D))+2
>run
?5
3.14159265359

I can get one thousand digits in little more than one fourth of a second on my old desktop computer, using that algorithm:

[size=x-small][font=Courier]n = 480

How do you get all those digits? There's only one DISP statement. What is the precision of the interpreter on your PC?

Gerson W. Barbosa · 01-06-2021, 02:41 AM

(01-06-2021 02:10 AM)toml_12953 Wrote:
(01-04-2021 08:41 PM)Gerson W. Barbosa Wrote: The algorithm is very simple, as tested on the HP-75C:

10 INPUT N
15 C=0
20 D=6*N*(N+1)+3/2
25 E=8*N
30 M=4*N*N-1
35 S=1
40 W=0
45 FOR K=1 TO N
50 W=W-S/M
55 C=-M*(M+1)/(4*(C+D))
60 M=M-E+4
65 D=D-E
70 E=E-8
75 S=-S
80 NEXT K
85 DISP S*(4*W+1/(C+D))+2
>run
?5
3.14159265359

I can get one thousand digits in little more than one fourth of a second on my old desktop computer, using that algorithm:

[size=x-small][font=Courier]n = 480

How do you get all those digits? There's only one DISP statement. What is the precision of the interpreter on your PC?

Sorry, I forgot to mention I was using Decimal BASIC at maximum precision (1000 digits).

ttw · 01-06-2021, 03:44 AM

One problem for all methods is that the last few iterations may dominate in that these must be done to really large precision. Some years ago (about the time of the first trillion digits computation) it was noted that the last division of sizes around 1,000,000, leading to the 1,000,000,000,000 digits takes most of the time in the computation. It completely dominated (at that time) the computation of the previous terms.

Fast (fourth order) methods were not much faster than first order or second order computations.

The online stuff nice (but it's base 16) and requires a trillion digit base conversion at the end.

Gerson W. Barbosa · (This post was last modified: 01-07-2021 12:58 PM by Gerson W. Barbosa.)

(01-06-2021 01:32 AM)Albert Chan Wrote:
Quote:\(\lim_{n\rightarrow \infty } \left [ \frac{1}{1\cdot 3}-\frac{1}{3\cdot 5}+\frac{1}{5\cdot 7}-\frac{1}{7\cdot 9}+\frac{1}{9\cdot 11}-\frac{1}{11\cdot 13}+\frac{1}{13\cdot 15}-\frac{1}{15 \cdot 17}\pm \cdots +\frac{(-1)^{n-1}}{4n^{2}-1}\right ]=\frac{\pi }{4}-\frac{1}{2}\)

Proof:

atan(x) = x - x^3/3 + x^5/5 - x^7/7 + ...
atan(1) = pi/4 = 1 - 1/3 + 1/5 - 1/7 + ...

2*I = 2/(1*3) - 2/(3*5) + 2/(5*7) - 2/(7*9) + ...
= (1-1/3) - (1/3-1/5) + (1/5-1/7) - (1/7-1/9) + ...
= (1-1/3+1/5-1/7+...) + (-1/3+1/5-1/7+...)
= pi/4 + (pi/4 - 1)

→ I = pi/4 - 1/2

Great! The continued fraction part is left unproved, however. But that would not be easy, I presume.

Albert Chan · (This post was last modified: 01-09-2021 08:35 PM by Albert Chan.)

(01-07-2021 12:57 PM)Gerson W. Barbosa Wrote: Great! The continued fraction part is left unproved, however. But that would not be easy, I presume.

"Wasicki's formula" had first term converge to pi, correction converge to 0. QED

As to the gain of 25/12 digits per term, it tested OK even with 10,000 digits precision. Smile

Code:

def test_pi(n, cf=0, c=0, s=0):

    cf += n         # extra cf "terms"

    b = 2*n*(n+1) + 4*cf*(cf+1) + mpfr(1.5)

    for k in range(cf,n,-1):

        t = mpfr(4*k*k-1)

        c = k*k*t/(b-c)

        b = b - 8*k

    for k in range(n,0,-1):

        t = mpfr(4*k*k-1)

        s = 1/t - s

        c = k*k*t/(b-c)

        b = b - 8*k

    return (s+0.5)*4 + (-1)**n/(b-c)

>>> from gmpy2 import *
>>> get_context().precision = 34000 # > 10000 dec. digits
>>> pi = const_pi()
>>>
>>> for n in range(1000,4001,1000): print n, format(pi - test_pi(n), 'g')
...
1000 2.19966e-2091
2000 2.92237e-4181
3000 3.88382e-6271
4000 5.16203e-8361

Gerson W. Barbosa · 01-08-2021, 01:05 AM

(01-07-2021 09:56 PM)Albert Chan Wrote:
(01-07-2021 12:57 PM)Gerson W. Barbosa Wrote: Great! The continued fraction part is left unproved, however. But that would not be easy, I presume.

"Wasicki's formula" had first term converge to pi, correction converge to 0. QED

Were the correction term simply 1/(2n*(n + 1)) then this proof would still hold. Perhaps I should say “prove that the continuous fraction is equivalent to the optimal correction polynomial expression to the series”. That would be really hard, at least for me. If you find a valid proof then the formula should be renamed as “Wasicki-Chan Formula” or more appropriately “Chan-Wasicki Formula”. Anyone can easily find such formulae. Anyone or anything – even Mathematica and W/A can, I think. But only a few can prove them. Anyway nowadays these formulae are just useless mathematical curiosities.

(01-07-2021 09:56 PM)Albert Chan Wrote: As to the gain of 25/12 digits per term, it tested OK even with 10,000 digits precision.

Thank you for performing these tests. Much appreciated!

Gerson W. Barbosa · 01-08-2021, 01:24 AM

(01-07-2021 09:56 PM)Albert Chan Wrote:
(01-07-2021 12:57 PM)Gerson W. Barbosa Wrote: Great! The continued fraction part is left unproved, however. But that would not be easy, I presume.

"Wasicki's formula" had first term converge to pi, correction converge to 0. QED

Were the correction term simply 1/(2n*(n + 1)) then this proof would still hold. Perhaps I should say “prove that the continuous fraction is equivalent to the optimal correction polynomial expression to the series”. That would be really hard, at least for me. If you find a valid proof then the formula should be renamed as “Wasicki-Chan Formula” or more appropriately “Chan-Wasicki Formula”. Anyone can easily find such formulae. Anyone or anything – even Mathematica and W/A can, I think. But only a few can prove them. Anyway nowadays these formulae are just useless mathematical curiosities.

(01-07-2021 09:56 PM)Albert Chan Wrote: As to the gain of 25/12 digits per term, it tested OK even with 10,000 digits precision.

Thank you for performing these tests. Much appreciated!

Albert Chan · (This post was last modified: 11-06-2021 04:21 PM by Albert Chan.)

Hi, Gerson W. Barbosa

I don't know how to quantify "optimal" correction. May be you want to proof this instead ?

\(\pi = \displaystyle\left( 2 + 4\sum_{k=1}^{n}\frac{(-1)^{k-1}}{4k^{2}-1} \right )
+ \lim_{k\rightarrow \infty }\left(\frac{(-1)^{n}} {b\;-\;\frac{3}{\large b+8\;-\;
\frac{60}{\frac{\ddots }{b+4k(k-1)\;-\;\frac{k^2 (4k^2-1)}{ b+4k(k+1) }}}}}\right)
,\quad b = 2n(n+1) + {3\over2}\)

This is not a proof, but I am convinced that above is true. Nice work Smile

First, try n=0, i.e. summation term = 0.
Relied only on limit, convergence rate is bad, but it is indeed approach pi.

>>> from gmpy2 import *
>>> get_context().precision = 1000 # ~ 300 dec. digits
>>> pi = const_pi()
>>> for i in range(7): print i, format(pi - test_pi(n = 0, cf = 10**i), 'g')
...
0 0.297148
1 0.0675375
2 0.00772836
3 0.000784124
4 7.85271e-05
5 7.85385e-06
6 7.85397e-07

Now, try combination of summation and limit.

>>> for i in range(7): print i, format(pi - test_pi(n = 9, cf = 10**i), 'g')
...
0 -5.30332e-21
1 -1.07719e-25
2 -1.47115e-39
3 -8.06555e-58
4 -9.64229e-77
5 -9.81688e-96
6 -9.83452e-115

>>> for i in range(7): print i, format(pi - test_pi(n = 10, cf = 10**i), 'g')
...
0 4.18442e-23
1 5.76716e-28
2 6.21078e-43
3 4.82312e-63
4 5.98619e-84
5 6.11764e-105
6 6.13095e-126

Note: test_pi() from previous post, modified to allow extra cf terms, beyond n.
Note: unlike summation term, CF correction convergence is one-sided.

Gerson W. Barbosa · (This post was last modified: 01-10-2021 03:12 PM by Gerson W. Barbosa.)

Hello, Albert Chan,

(01-09-2021 09:22 PM)Albert Chan Wrote: First, try n=0, i.e. summation term = 0.
Relied only on limit, convergence rate is bad, but it is indeed approach pi.

>>> from gmpy2 import *
>>> get_context().precision = 1000 # ~ 300 dec. digits
>>> pi = const_pi()
>>> for i in range(7): print i, format(pi - test_pi(n = 0, cf = 10**i), 'g')
...
0 0.297148
1 0.0675375
2 0.00772836
3 0.000784124
4 7.85271e-05
5 7.85385e-06
6 7.85397e-07

I did that about seven years ago ( https://www.hpmuseum.org/forum/post-9194.html#pid9194 ). Then I noticed that optimal result (best approximation using the least number of terms) is obtained when the number of terms of the series is equal to the number of terms of the continued fraction.

ln(2) = 1 - 1/2 + 1/3 - 1/4 + ... + (-1)^(n + 1)/n + 1/((2*n+1) + CF)

CF: numeratorₖ = k*k; denominatorₖ = 2*n + 1; k = 1 to ∞

n = 2
k = 1

ln(2) ~ 1 - 1/2 + (-1)^(1 + 1)/((2*2+1)+1*1/(2*2+1))
~ 9/13
~ 0.69(23076923076923076923076923076923)

n = 2
k = 10^6

ln(2) ~ 0.693147180559945309417232121458(3016)

n = 2
k = 10^7

ln(2) ~ 0.6931471805599453094172321214581766

Here we need 10^7 terms of the continued fraction to correct ln(2) ~ 1/2 to 34 digits. On the other hand we get the same result with n = 22 and k = 22.

Perhaps the correction CF for the Wallis Product is easier to prove, since the equivalent polynomial expressions of the first few terms of the CF match those described in this paper. More details in the beginning of this post.

Albert Chan · (This post was last modified: 11-07-2021 11:09 AM by Albert Chan.)

This pi formula is very similar to zeta(2) alternating sum formula

(11-05-2021 03:55 PM)Albert Chan Wrote:
(10-31-2021 03:40 PM)Albert Chan Wrote: Let s = (-1)^n, x=n*n+n+1, T's = triangular number

ζ(2) ≈ (2 - 2/2^2 + 2/3^2 - ... + s*2/n^2)

- s/(x - 1/(x+4*T1 - 2^4/(x+4*T2 - 3^4/(x+4*T3 - 4^4/(x+4*T4 - ...
...
We can spend more time summing the alternating series, less time doing CF corrections.
In the end, we may come out ahead.

As an experiment, we cut down CF terms in half. (p pairs alternating series, do p CF terms)
Bonus: summing alternating series pairwise also made final sum more accurate.

Apply the same idea here, we have this for estimating PI

Code:

OPTION ARITHMETIC DECIMAL_HIGH

10 INPUT  PROMPT "p (pairs) = ":p

   LET x = (8*p+4)*p + 1.5

   LET a = 0

   LET b = 0

   FOR i = p TO 2 STEP -1

      LET t = i*i

      LET a = a + 4/((i-0.75)*(16*t-1))

      LET b = t*(4*t-1) / (x+4*(t+i)-b)

   NEXT i

   LET r = 1/(x-3/(x+8-b)) + a + 46/15   

   PRINT "Accurate digits ="; 1-LOG10(MAX(ABS(PI-r),1E-1000))  

   GOTO 10

END

p (pairs) = 100
Accurate digits = 308.79106444764398
p (pairs) = 101
Accurate digits = 311.85336518362854
p (pairs) = 300
Accurate digits = 921.2386113713281
p (pairs) = 301
Accurate digits = 924.30082763127136

For p=326, it reached 1000 digits full precision.

Albert Chan · (This post was last modified: 11-07-2021 12:45 PM by Albert Chan.)

Wallis product part is *very slowly* converging; probably not a good idea to cut down CF terms.

Another issue is to limit rounding erors of Wallis product.
To improve accuracy, we keep track of partial product, less 1.

(1+a)*(1+1/t) - 1 = a + (1+a)/t

To recover pi, we just add back 1 (after loops finished)

pi ≈ (1+a)*(1+b)*8/3 = (a*b + b + (a-1/8)) * 8/3 + 3 (*)

Code:

OPTION ARITHMETIC DECIMAL_HIGH

10 INPUT  PROMPT "n = ":n

   LET x = 8*n+4

   LET a = 0

   LET b = 0

   FOR i = n TO 2 STEP -1

      LET t = 4*i*i-1

      LET a = a + (1+a)/t

      LET b = t / (x+b)

   NEXT i

   LET b = 2/(x-1+3/(x+b)) ! CF correction

   LET r = (a*b+b+(a-0.125)) * 8/3 + 3

   PRINT "Accurate digits ="; 1-LOG10(MAX(ABS(PI-r),1E-1000))

   GOTO 10

END

n = 100
Accurate digits = 211.18296908289588
n = 101
Accurate digits = 213.27288787760364
n = 300
Accurate digits = 629.16112325948984
n = 301
Accurate digits = 631.25100446831984

For n=478, it reached 1000 digits full precision.

Previous version take more total terms (3p=3*326=978) vs here (2n=2*478=956)
But, pairwise sum can be greatly simplified → Previous version is faster, by about 30%

(*) By Sterbenz lemma, term (a-1/8) is exact

n<2: it never enter the loop, thus a = 0
n=2: a = 16/15 - 1 = 1/15, within [1/16,1/4]
With more products, a is increasing. limit(a, n=inf) = 3/8*pi - 1 ≈ 1/5.6 ≤ 1/4