Yet another π formula

[Albert Chan]

(11-24-2024 08:17 PM)Gerson W. Barbosa Wrote:  Mine is slower: ~ 1'33"

Not really. It is comparing apples and oranges.

Your calculations sum 5 terms (ignored the 1/2), and correct with CF factor 6 levels deep.
Summation of infinite alternative series do 0 .. (n=10) = 11 terms, and 0 .. (d=7) = 8 (Δ^k/2^k+1) terms.

To match same number of 'terms', we set n = 5-1 = 4 , d = 6-1 = 5
Note: S=S1-S2 because n is even, otherwise it would be S=S1+S2

lua> S1 = 1/(1*3) - 1/(3*5) + 1/(5*7) - 1/(7*9) + 1/(9*11)
lua> S2 = 0.5/(11*13) * (1 + 2/15*(1 + 3/17*(1 + 4/19*(1 + 5/21*(1 + 6/23)))))
lua> (S1-S2)*4 + 2
3.1415943802536885

This is only half as many good digits as your 5 terms + correction pi estimate.
Even with n=10, d=7, getting pi estimate ≈ 3.141592654 is still not as good.

lua> n = 10
lua> S1 = 0
lua> for i=n,0,-1 do S1=1/((2*i+2)^2-1) - S1 end
lua> S1
0.28633982191661556

For S2, we can do calculations inside out, but I rather do this left to right.
We can Euler-transform it as continued fraction, and do convergents of continued fraction.

Code:

function dfc2f_euler(b0)    -- b0 + a0*(1+a1*(1+a2*(1+a3*(...
    local f, a = dfc2f(b0)
    return function(a2)
        a, a2 = a2, a
        return a2 and f(-a,1+a) or f(a,1)
    end
end

lua> f = dfc2f_euler(0) -- euler-transform S2 as CF
lua> S2 = f(0.5/((2*n+3)*(2*n+5))) -- 0th term
lua> for k=1,7 do S2 = f( (k+1) / (2*(n+k)+5) ) end
lua> S2
0.000941658270090828
lua> (S1-S2)*4 + 2
3.141592654586099

We would need to push S2 d from 7 to 20 to reach full precision

lua> for i=8,20 do S2 = f( (i+1) / (2*(n+i)+5) ) end
lua> S2
0.0009416585191672306
lua> (S1-S2)*4 + 2
3.141592653589793

[Gerson W. Barbosa]

(11-25-2024 01:47 PM)Albert Chan Wrote:  

(11-24-2024 08:17 PM)Gerson W. Barbosa Wrote:  Mine is slower: ~ 1'33"

Not really. It is comparing apples and oranges.

Indeed apples and apples (or oranges and oranges). 1'33" is the time it took for Thomas Klemm’s program to run on my HP-15C. I only meant my HP-15C is slower than his.
I thought of writing an HP-15C program to compare both methods, but I didn’t.
Sorry for the confusion!

[Thomas Klemm]

(11-25-2024 04:45 PM)Gerson W. Barbosa Wrote:  I only meant my HP-15C is slower than his.

To be honest, I was using the following emulator:

Quote:HP-15C
Version 1.0.2 (1001)

Hewlett-Packard Voyager Series
Retina Graphics Interface & Enhancements
Copyright © 2022 Mark H. Shin

based on nonpareil for MacOS
Copyright 2005-2012 Maciej Bartosiak

based on Nonpareil 0.77 for Linux
Copyright 2004-2006 Eric Smith

I had assumed that it was running at the same time as the original.
It is apparently running a bit faster, but still in the same ballpark.

[Albert Chan]

Σ((-1)^k/(2k+1), k=0..inf) = atan(1) = pi/4
Σ((-1)^k/(2k+1)/(2k+3), k=0..inf) = (pi/4 - 0.5/(1)) / 1 = pi/4 - 1/2
Σ((-1)^k/(2k+1)/(2k+3)/(2k+5), k=0..inf) = ((pi/4 - 1/2) - 0.5/(1*3)) / 2 = pi/8 - 1/3
...

With this pattern, we can build formula for left hand side, S = pi / c1 - c2

lua> N = 40 -- number of consecutive odd numbers for LHS denominator
lua> c1, c2, t = 4, 0, 1 -- start from atan(1) = pi/4
lua> for i=2,N do c1=c1*(i-1); c2=(c2+0.5/t)/(i-1); t=t*(2*i-1) end

pi = c1*c2 + c1*S

lua> c1*c2 -- close enough to pi, we only need tiny correction.
3.141592653588783
lua> c1*c2 + c1 * (1-1/(2*N+1))/t -- 2 terms for S
3.141592653589793

[Albert Chan]

From previous post, if we ignore S correction, with goal of c1 * c2 = pi, we get this.

\(\displaystyle \frac{\pi}{2} = \sum _{k=0}^∞ \frac {k!}{(2k+1)!!}\)

Formula approximately doubled its accuracy after each sum.
This explained why previous post, with N=40, c1*c2 get about 40 bits accuracy of pi.

---

To speed up convergence, we can "double" the last term.

pi/2 = 1 + 1/3 + 2/15 + 6/105 + 24/945 + ...
       = (2-1) + (2/3-1/3) + (4/15-2/15) + (12/105-6/105) + (48/945-24/945) + ...
       = 2 - (1-2/3) - (1/3-4/15) - (2/15-12/105) - (6/105-48/945) - ...

k!/(2k+1)!! - 2*(k+1)!/(2k+3)!! = k!/(2k+3)!!

Again, we have a pattern here

Code:

pi/2 = 0 + sum(k!/(2k+1)!!, k=0..inf)
pi/2 = 2 - sum(k!/(2k+3)!!, k=0..inf)             -- 0 + 2/1 = 2
pi/2 = 4/3 + 3!! * sum(k!/(2k+5)!!, k=0..inf)     -- 2 - 2/3 = 4/3
pi/2 = 26/15 - 5!! * sum(k!/(2k+7)!!, k=0..inf)   -- 4/3 + 2/5 = 26/15
pi/2 = 152/105 + 7!! * sum(k!/(2k+9)!!, k=0..inf) -- 26/15 - 2/7 = 152/105
...

Ignore the sum part, we goes full circle

pi/2 = 2/1 - 2/3 + 2/5 - 2/7 + ... = 2 * atan(1)

[Albert Chan]

(11-25-2024 01:47 PM)Albert Chan Wrote:  To match same number of 'terms', we set n = 5-1 = 4 , d = 6-1 = 5
Note: S=S1-S2 because n is even, otherwise it would be S=S1+S2

lua> S1 = 1/(1*3) - 1/(3*5) + 1/(5*7) - 1/(7*9) + 1/(9*11)
lua> S2 = 0.5/(11*13) * (1 + 2/15*(1 + 3/17*(1 + 4/19*(1 + 5/21*(1 + 6/23)))))
lua> (S1-S2)*4 + 2
3.1415943802536885

FYI, above convergence acceleration scheme is EXACTLY the same as my previous post.

pi/2 = (2 - 2/3 + 2/5 - 2/7 + 2/9) - 9!! * sum(k!/(2k+11)!!, k=0..inf)

lua> s1 = (2 - 2/3 + 2/5 - 2/7 + 2/9)
lua> s2 = 1/11*(1 + 1/13*(1 + 2/15*(1 + 3/17*(1 + 4/19*(1 + 5/21*(1 + 6/23))))))
lua> (s1-s2)*2
3.141594380253689

[EdS2]

I like this double factorial approach. Here it is in BBC Basic:

Code:

PRINT "Calculate Pi with double factorial"
L=21:N=1:D=1:S=1
FOR K=1 TO L
N=N*K
D=D*(2*K+1)
S=S+N/D
PRINT 2*(S+N/D)
NEXT

[AnnoyedOne]

Or just use 3 (close enough for most real-life applications) or see

https://www.piday.org/million/

A1

[C.Ret]

(11-26-2024 02:19 PM)AnnoyedOne Wrote:  https://www.piday.org/million/

Thank you for the link!

I just found my date of birth in the decimal of PI.
Does that mean anything?
How do I find my date of death?

[AnnoyedOne]

(11-26-2024 03:42 PM)C.Ret Wrote:  How do I find my date of death?

No idea. But you might be eating pi(e) at the time

A1

[Albert Chan]

Hi, EdS2

Here are improvements over your BBC Basic code.

1. k! / (2k+1)!!, numerator and denominator get huge very quickly
    We can let T = N/D, eliminate possible overflow issue.
2. 2k+1 replaced by 2k+p with p=1, for usage flexibility (see example below)
3. Instead of add last term for better pi estimate, I use Aitken extrapolated result.
4. Added plain 2*S pi estimate for comparison

Code:

10 print "Calculate Pi with double factorial"
20 L=21:P=1:T=1/P:S=T
30 for K=1 to L
40 D=2+P/K:T=T/D:S=S+T
50 print 2*(S+T/(D-1)), 2*S
60 next

BTW, my original Aitken code is this.

Code:

40 T0=T:T=T*K/(2*K+P):S=S+T
50 print 2*(S+T*T/(T0-T))

But I like D = 2+P/K ≈ 2 version better.

It clearly showed T=T/D shrink to about half, thus sum accuaracy doubled per step.
Also, it showed why add last term T to S improve estimate. Correction T/(D-1) ≈ T

(11-26-2024 01:28 PM)Albert Chan Wrote:  pi/2 = (2 - 2/3 + 2/5 - 2/7 + 2/9) - 9!! * sum(k!/(2k+11)!!, k=0..inf)

lua> s1 = (2 - 2/3 + 2/5 - 2/7 + 2/9)
lua> s2 = 1/11*(1 + 1/13*(1 + 2/15*(1 + 3/17*(1 + 4/19*(1 + 5/21*(1 + 6/23))))))
lua> (s1-s2)*2
3.141594380253689

For this example, we change 1 line (run on DOS ubasic)
Note: we start from k=1, thus k=0 correction goes to S.

20 L=6:P=11:T=-1/P:S=2-2/3+2/5-2/7+2/9+T
run
Calculate Pi with double factorial
 3.1427128427128427129    3.1438783438783438784
 3.1417266494189571114    3.1420135420135420136
 3.1416139416139416141    3.1416844593315181551
 3.1415967039496451264    3.1416151787668815534
 3.1415935285904326155    3.1415986833943490292
 3.1415928614981610796    3.1415943802536883708
OK

[EdS2]

Nice, thanks!