祖冲之的割圆术与圆周率

lihpb00

收敛太慢了

nyy

lihpb00 发表于 2023-12-23 23:09
收敛太慢了

看看马青公式，计算前5项，得到祖冲之的精度。
\[x-\frac{x^3}{3}+\frac{x^5}{5}-\frac{x^7}{7}+\frac{x^9}{9}\]

1. Clear["Global`*"];(*清除所有变量*)
   
2. $MaxExtraPrecision=500;(*最大额外精度*)
   
3. aa=Sum[(-1)^(k+1)*x^(2k-1)/(2k-1),{k,1,4}];(*计算前n项和*)
   
4. pi=(16*aa/.x->1/5)-(4*aa/.x->1/239)(*用马青公式计算圆周率*)
   
5. N[pi-Pi,100](*计算差值*)
   
6. N[pi,100](*计算近似值*)
   

复制代码

得到
3.141592682404399517240259836073575860489757996720135506301754703591490348774810995008612800369424553

如果计算72项，大概就能得到圆周率的100位
\[x-\frac{x^3}{3}+\frac{x^5}{5}-\frac{x^7}{7}+\frac{x^9}{9}-\frac{x^{11}}{11}+\frac{x^{13}}{13}-\frac{x^{15}}{15}+\frac{x^{17}}{17}-\frac{x^{19}}{19}+\frac{x^{21}}{21}-\frac{x^{23}}{23}+\frac{x^{25}}{25}-\frac{x^{27}}{27}+\frac{x^{29}}{29}-\frac{x^{31}}{31}+\frac{x^{33}}{33}-\frac{x^{35}}{35}+\frac{x^{37}}{37}-\frac{x^{39}}{39}+\frac{x^{41}}{41}-\frac{x^{43}}{43}+\frac{x^{45}}{45}-\frac{x^{47}}{47}+\frac{x^{49}}{49}-\frac{x^{51}}{51}+\frac{x^{53}}{53}-\frac{x^{55}}{55}+\frac{x^{57}}{57}-\frac{x^{59}}{59}+\frac{x^{61}}{61}-\frac{x^{63}}{63}+\frac{x^{65}}{65}-\frac{x^{67}}{67}+\frac{x^{69}}{69}-\frac{x^{71}}{71}+\frac{x^{73}}{73}-\frac{x^{75}}{75}+\frac{x^{77}}{77}-\frac{x^{79}}{79}+\frac{x^{81}}{81}-\frac{x^{83}}{83}+\frac{x^{85}}{85}-\frac{x^{87}}{87}+\frac{x^{89}}{89}-\frac{x^{91}}{91}+\frac{x^{93}}{93}-\frac{x^{95}}{95}+\frac{x^{97}}{97}-\frac{x^{99}}{99}+\frac{x^{101}}{101}-\frac{x^{103}}{103}+\frac{x^{105}}{105}-\frac{x^{107}}{107}+\frac{x^{109}}{109}-\frac{x^{111}}{111}+\frac{x^{113}}{113}-\frac{x^{115}}{115}+\frac{x^{117}}{117}-\frac{x^{119}}{119}+\frac{x^{121}}{121}-\frac{x^{123}}{123}+\frac{x^{125}}{125}-\frac{x^{127}}{127}+\frac{x^{129}}{129}-\frac{x^{131}}{131}+\frac{x^{133}}{133}-\frac{x^{135}}{135}+\frac{x^{137}}{137}-\frac{x^{139}}{139}+\frac{x^{141}}{141}-\frac{x^{143}}{143}\]
差值
-4.7347325070038906601393332841512003488000614441363296637426861454146\
16030131177859665455168151268452*10^-103

估计71项达不到精度，但是我没仔细研究！

nyy

nyy 发表于 2023-12-25 12:12
看看马青公式，计算前5项，得到祖冲之的精度。
\[x-\frac{x^3}{3}+\frac{x^5}{5}-\frac{x^7}{7}+\frac{x^ ...

如果要计算到小数点后100位，那么需要2^168=3.741444191*10^(50)正多边形，代码如下

1. Clear["Global`*"];(*清除所有变量*)
   
2. a=Sqrt[2];(*正方形的初始边长*)
   
3. n=4;(*正方形边数*)
   
4. Do[a=Sqrt[2-Sqrt[4-a^2]];(*由n边形的边长,计算出2n边形的边长*)
   
5.    n=n*2;(*边数从n变成2n*)
   
6.    pi=a*n/2,(*周长除以直径=圆周率*)
   
7. {k,1,166}]
   
8. N[pi-Pi,110]
   
9. aa=N[pi,103](*鲁道夫的圆周率*)
   
10. bb=N[Pi,103](*真实的圆周率*)
    

复制代码

结果
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067945
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982

需要连续开167次方！

nyy

nyy 发表于 2023-12-25 12:12
看看马青公式，计算前5项，得到祖冲之的精度。
\[x-\frac{x^3}{3}+\frac{x^5}{5}-\frac{x^7}{7}+\frac{x^ ...

利用马青公式\(48\arctan\frac{1}{18}+32\arctan\frac{1}{57}-20\arctan\frac{1}{239}\)
只需要计算前40项就可以了！

代码如下：

1. Clear["Global`*"];(*清除所有变量*)
   
2. $MaxExtraPrecision=100;(*最大额外精度*)
   
3. aa=Sum[(-1)^(k+1)*x^(2k-1)/(2k-1),{k,1,50}];(*计算前n项和*)
   
4. pi=(48*aa/.x->1/18)+(32*aa/.x->1/57)-(20*aa/.x->1/239)(*用马青公式计算圆周率*)
   
5. N[pi-Pi,103](*计算差值*)
   
6. N[pi,104](*计算近似值*)
   
7. N[Pi,104](*计算近似值*)
   

复制代码

结果如下
3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679809
3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821
可以看到，实际上精确到小数点后101位。

nyy

nyy 发表于 2023-12-29 11:47
利用马青公式\(48\arctan\frac{1}{18}+32\arctan\frac{1}{57}-20\arctan\frac{1}{239}\)
只需要计算前40 ...

利用马青公式
\(\frac{\pi}{4}=44\arctan\frac{1}{57}+7\arctan\frac{1}{239}-12\arctan\frac{1}{682}+24\arctan\frac{1}{12943}\)
只需要计算29项

1. Clear["Global`*"];(*清除所有变量*)
   
2. $MaxExtraPrecision=250;(*最大额外精度*)
   
3. aa=Sum[(-1)^(k+1)*x^(2k-1)/(2k-1),{k,1,29}];(*计算前n项和*)
   
4. (*用马青公式计算圆周率*)
   
5. pi=4*((44*aa/.x->1/57)+(7*aa/.x->1/239)-(12*aa/.x->1/682)+(24*aa/.x->1/12943));
   
6. N[pi-Pi,103](*计算差值*)
   
7. N[pi,104](*计算近似值*)
   
8. N[Pi,104](*计算近似值*)
   

复制代码

计算结果
3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679822
3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821

tdwhh

我偶然发现的一个圆周率公式,大家看下：一个简洁的圆周率算法 - happytdw的文章 - 知乎
https://zhuanlan.zhihu.com/p/675031866

nyy

tdwhh 发表于 2023-12-31 11:58
我偶然发现的一个圆周率公式,大家看下：一个简洁的圆周率算法 - happytdw的文章 - 知乎
https://zhuanlan.z ...

1. (*一个简洁的圆周率算法 https://zhuanlan.zhihu.com/p/675031866*)
   
2. Clear["Global`*"];(*清除所有变量*)
   
3. a0=1;a1=2;(*a数列的初始值*)
   
4. b0=1;b1=2;(*b数列的初始值*)
   
5. (*利用循环进行递推*)
   
6. Do[a2=a1+n*a0;(*计算递推结果*)
   
7.    b2=b1+n*b0+1;
   
8.    {a0,a1}={a1,a2};(*重新赋值*)
   
9.    {b0,b1}={b1,b2},
   
10. {n,2,400}]
    
11. pi=2*(b2/a2)^2//Simplify;
    
12. N[pi-Pi,100]
    

复制代码

pi的计算结果
62025000939678909516687735193323735679065222430471033763381665839374761995927718792304418729658751210049860488860911551680939456432144359841211072335895952147921464697513968242090392315718012516015590713627465244700797183151760707771499894364547214710591597600846408933165356599304275471466095871602326246357731878903070418372922554993470258688672327713486649550161027062315402388427258430020421815659675995218714419996458444204554859904759312408758535748607611738041254470374876593240980407603750246924235349756903590684663913440602428379474402975747442120724740808012952797639745546189266638840413443786459171826980252856226982197660944851286153199683302766825452008610849213435634382955651685398647960622185880437762598492223473579124892773167361353256715117973875293103124963822411256705801871887131620488818568592134867440992370758611285265299055526977458774244746074556800428542081/19743170989658702829708591285591630717562365062427886961099403631133476119547148422881762950355106763464658566139692506973077016062198092073282093446192327692126784166349645485356214315502798552923766493857569813524669593818230628367479232839474331330314390642177373411824165311786333767561565008246688177567595388974795002389384579994311747141634363559263064669613868071103416913625011554865450732523170394255632143364262004632226712388029603885440999693919866053401215899783057464554443190561413222910839475667718942779276999809824558685043390141945559782178717329802237035412583201025565795456311920839240762285416938939318009154045190203091851759480430913872416213859324276124729267616043236366753593124787498572461727221865110910684771119362709648153306017044421356246752838000108858440613003052598961859241908248055209419156532740031605893034256340075591729806631243285637929893888

误差
4.932603769219693390126654336135238877664753636994259082389243155687944028218791267849224234479444015*10^-17

根据我的计算结果，收敛速度好像并不快，不知道原理是什么

nyy

tdwhh 发表于 2023-12-31 11:58
我偶然发现的一个圆周率公式,大家看下：一个简洁的圆周率算法 - happytdw的文章 - 知乎
https://zhuanlan.z ...

利用这儿的代码
https://bbs.emath.ac.cn/forum.ph ... 17586&pid=98458

1. (*拉马努金计算圆周率*)
   
2. Clear["Global`*"];(*清除所有变量*)
   
3. aa=Sum[(4k)!*(1103+26390k)/((k!)^4*396^(4k)),{k,0,2}]//Simplify
   
4. pi=9801/(2*Sqrt[2])/aa
   
5. N[pi-Pi,20]
   

复制代码

圆周率的近似表达式=
\[\frac{2286635172367940241408 \sqrt{2}}{1029347477390786609545}\]

与精确值的误差
5.68242325601395950808108722896*10^-24

这个表达式，比上面的表达式计算要小很多，且精度要搞几个有效数字。

所以我推测你的算法不仅仅计算量大且精度小。
这个是我的论据。

nyy

nyy 发表于 2024-1-3 13:26
利用这儿的代码
https://bbs.emath.ac.cn/forum.php?mod=redirect&goto=findpost&ptid=17586&pid=98458

计算前五项的值

1. (*拉马努金计算圆周率*)
   
2. Clear["Global`*"];(*清除所有变量*)
   
3. aa=Table[9801/(2*Sqrt[2])/Sum[(4k)!*(1103+26390k)/((k!)^4*396^(4k)),{k,0,n}],{n,0,5}]
   

复制代码

计算结果
\[\begin{array}{c}
\frac{9801}{2206 \sqrt{2}} \\
\frac{2510613731736 \sqrt{2}}{1130173253125} \\
\frac{2286635172367940241408 \sqrt{2}}{1029347477390786609545} \\
\frac{17252765328978109815564789153792 \sqrt{2}}{7766473062254307011793347201855} \\
\frac{1131379202490552979877435552947122965839872 \sqrt{2}}{509299577881529611662930757403081523769055} \\
\frac{128805730098892711723125911845114081418091536842752 \sqrt{2}}{57982950211280781944919792648021104999982386829481} \\
\end{array}\]