Re: wp34s program Message #21 Posted by Gerson W. Barbosa on 16 July 2011, 11:27 p.m., in response to message #20 by Gerson W. Barbosa
Quote:
I've come up with a variation of Archimedes' method that will produce about 1.8 digits per iteration.
It is essentially the Archimedes' method, only an attempt has been made to improve the approximations. (He did it this way).
Using trigonometry, the Archimedes' lower and upper bounds for pi are, respectively:
a = x*sin(pi/x)
b = x*tan(pi/x)
where x=3*2^k, k=1, 2,...
Let's consider the following limit:
x*tan(pi/x)  pi
lim  = 2
x>inf pi  x*sin(pi/x)
then the following weighted arithmetic mean of the bounds gives successively better approximations to pi at each iteration:
2*a + b
wam =  ~ pi
3
Likewise, the following weighted geometric mean appears to hold:
wgm = (a^2*b)^(1/3) ~ pi
also
wam(x)  pi 9
lim  = 
x>inf wgm(x)  pi 4
this leads to
9*wgm(x)  4*wam(x)
pi = , when x equals infinity
5
or
27*(a^2*b)^(1/3)  4*(2*a + b)
pi ~ 
15
The demonstration may not be rigorously correct, nevertheless we can use the result to write the following wp34s program:
001 LBL 77
002 STO 00
003 3
004 SQRT
005 STO+ X
006 FILL
007 4
008 RCL/ L
009 /
010 *
011 SQRT
012 DSE 00
013 SKIP 01
014 SKIP 06
015 STO 01
016 
017 STO+ X
018 FILL
019 RCL 01
020 BACK 10
021 RCL* X
022 RCL L
023 Rv
024 *
025 CUBERT
026 2
027 7
028 *
029 RCL Z
030 RCL+ X
031 RCL+ Z
032 4
033 *
034 
035 1
036 5
037 /
038 RTN
1 XEQ 77 > 3.141460911773500
2 XEQ 77 > 3.141590947963119
3 XEQ 77 > 3.141592628123497
4 XEQ 77 > 3.141592653196343
5 XEQ 77 > 3.141592653583662
6 XEQ 77 > 3.141592653589791
An equivalent C++ program is:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
int main(int argc, char *argv[])
{
char i, m;
int n = 6;
double gm, hm, wam, wgm; // start with inscribed and circumscribed
scanf("%d", &m); // hexagons (diameter of the circumference = 1)
system("cls"); // hm: harmonic mean
hm = 2*sqrt(3); // gm: geometric mean
gm = 3*hm/4; // 2*sqrt(3): perimeter of the circumscribed
printf("\ ngon\t wgm\t wam\t\t(9*wgm4*wam)/5 \n\n"); // hexagon
for (i = 1; i <= m; i++) // gm(3/4*2*sqrt(3),2*sqrt(3)) = 3, this is the
{ // perimeter of the inscribed hexagon
gm = sqrt(gm*hm); // wam: weighted arithmetic mean of hm and gm
wam = (2*gm+hm)/3; // wgm: weighted geometric mean of hm and gm
wgm = pow(gm*gm*hm,1/3.); // weights of gm and hm: 2 and 1, respectively
printf("%12d %.16f %.16f %.16f\n",n,wgm,wam,(9*wgm4*wam)/5); // the weighted arithmetic mean of wgm and wam
hm *= 2*gm/(gm+hm); // (weights 9 and 4, respectively) will approximate
n += n; // pi at a rate of 1.8 correct digits per iteration.
}
printf("\n");
system("pause");
return 0;
}
Output for nine iterations:
ngon wgm wam (9*wgm4*wam)/5
6 3.1473451902649443 3.1547005383792515 3.1414609117734984
12 3.1419279179993587 3.1423491305446571 3.1415909479631199
24 3.1416132628330500 3.1416390562199923 3.1415926281234960
48 3.1415939364016969 3.1415955404083897 3.1415926531963425
96 3.1415927336837215 3.1415928338087955 3.1415926535836625
192 3.1415926585943867 3.1415926648502488 3.1415926535896972
384 3.1415926539025598 3.1415926542935209 3.1415926535897909
768 3.1415926536093401 3.1415926536337748 3.1415926535897922
1536 3.1415926535910144 3.1415926535925416 3.1415926535897927
=======================================
Update:
9*a*c 7*(a  4*c)
10*pi ~   
4*a  c 3
where
a = lower bound
c = next lower bound
avoids the cube root exraction and gives out the same accuracy (1.8 decimal places per iteration).
=======================================
Update:
Found a relationship between the previous formulas and had WA solve the equation (solve (p(27*(a^2*b)^(1/3)4*(2*a+b))/15 )/((9*a*c/(4*ac)7/3*(a4*c))/10p)==128/19 for p) and simplify the result:
p ~ (513*(a^2*b)^(1/3) + 24*a*((288*a)/(c  4*a)  97)  76*b + 1792*c)/2205
where
a = lower bound
b = upper bound
c = the next lower bound
This gives 2.4 decimal digits per iteration.
Archimedes' 96gon is enough for 15 correct digits:
3.14159265358979259
The following procedure can be used to check the formula without running all the iterations:
n:=50;
a:=3*2^n*sin(pi/(3*2^n));
b:=3*2^n*tan(pi/(3*2^n));
c:=3*2^(n+1)*sin(pi/(3*2^(n+1)));
p:=(513*(a^2*b)^(1/3)+24*a*((288*a)/(c4*a)97)76*b+1792*c)/2205;
More details here.
Edited: 22 July 2011, 4:58 p.m. after one or more responses were posted
