+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html)
+--- Forum: General Forum (/forum-4.html)
+--- Thread: Yet another π formula (/thread-16122.html)
Pages: 1 2
RE: Yet another π formula - Albert Chan - 11-25-2024 01:47 PM
(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/2k+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
endlua> 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
RE: Yet another π formula - Gerson W. Barbosa - 11-25-2024 04:45 PM
(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!
RE: Yet another π formula - Thomas Klemm - 11-25-2024 06:47 PM
(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.
RE: Yet another π formula - Albert Chan - 11-25-2024 10:01 PM
Σ((-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
RE: Yet another π formula - Albert Chan - 11-26-2024 01:32 AM
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)
RE: Yet another π formula - Albert Chan - 11-26-2024 01:28 PM
(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
RE: Yet another π formula - EdS2 - 11-26-2024 02:09 PM
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)
NEXTRE: Yet another π formula - AnnoyedOne - 11-26-2024 02:19 PM
Or just use 3 (close enough for most real-life applications) or see
https://www.piday.org/million/
A1
RE: Yet another π formula - C.Ret - 11-26-2024 03:42 PM
(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?
RE: Yet another π formula - AnnoyedOne - 11-26-2024 05:37 PM
(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
RE: Yet another π formula - Albert Chan - 11-26-2024 05:57 PM
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 nextBTW, 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
RE: Yet another π formula - EdS2 - 11-27-2024 09:14 AM
Nice, thanks!