我想的圆周率计算法，１万位十进制精度计算时间在数秒至数十秒左右

wayne

不过，这些信息其实就相当于什么都没说，

abiao

29# liangbch sin的效率我目前只能死算，一次次叠加，不过还好收敛的速度可以接受 叠加一次需要做一次全精度的除法（这个是大头），和几次乘法（计算位数不多），一次加法（可以忽略）

无心人

sin有AGM算法吧

LZC_314

以前看过一个程序，计算pi的，和大家分享一下 C++的，计算的挺快的，就是输出费时间

1. #include < stdio.h>
2. long b, c=2800, d, e, f[2801];
3. void main(){
4. while(b-c!=0){
5. f[b]=2000;
6. b++;
7. }
9. while(c!=0){
10. b=c;
11. d=f[b]*10000;
12. f[b]=d%(b*2-1);
13. d=d/(b*2-1);
14. --b;
15. while(b!=0){
16. d=d*b+f[b]*10000;
17. f[b]=d%(b*2-1);
18. d=d/(b*2-1);
19. --b;
20. }
21. c-=14;
22. printf("%.4d",e+d/10000);
23. e=d%10000;
24. }
25. }

复制代码

liangbch

明白了，楼主的方式本质上是牛顿迭代法，首先求的一个初值，然后逐步求精。仔细推算了下发现，如果第n次迭代可以精确到d位有效数字，第n+1次迭代则可以精确到3d位有效数数字。即若第n次迭代的误差是e,则第n+1次迭代则的误差为$1/6*e^3$. 对于直接使用泰勒公式求sin(x),当x在pi附近，收敛的仍然比较慢。可考虑先求sin(x/9), 然后利用正弦倍角公式$sin(9x)=9sin(x)-120sin^3(x)+432sin^5(x)-576sin^7(x)+256sin^9(x)$ 求得sin(x)的值。

liangbch

9# abiao 无理论是理论上还是时间上，我认为你的算法不会比Martin公式更快。 从理论上讲， 使用泰勒展开，计算sin(x）的复杂度约为O(d^3), d是要精确的位数。而Martin公式的复杂是O(d^2),而且编程非常简单。 从实践上，你说你的程序精确到1万位需要几秒到几十秒。而我在2000年编写的超级计算器中就包含了计算圆周率的功能，在我的这台老爷机(PM1.7,2005年购买）使用这个计算器计算1万位圆周率也只需0.56秒。 注：我的那个超级计算器可从http://download.csdn.net/download/liangbch/153669 下载，由于csdn换版，错误的把上载日期改为2006年，实际上我是在2001年上载的。

liangbch

这个算法的优雅之处在于简单（原理简单，易于理解），和不算差的性能，比起Martin公式，Ramanujan公式和其他乌七八糟的公式比起来简单得多了abiao 发表于 2010-8-17 21:04

本着科学认真的精神，还是评议一下吧，楼主不要见怪。 “不算差的性能”，我不敢苟同，单是计算sin(x)，你的算法的复杂度就达到O(d^3),而目前最好的算法几乎是O(d). "比起Martin公式，Ramanujan公式和其他乌七八糟的公式比起来简单得多了" 公式是简单了，但计算过程一点也不简单，Martin公式的计算非常简单，在100行C代码完全可以实现。Ramanujan是个数学天才，他的那些公式极富创意，非普通人可发现，而且收敛速度极快，目前几个计算圆周率的程序，其中好多用的就是Ramanujan的公式，另外一些用的是AGM算法。 而你的这个公式却非常简单，发现和证明这个公式只需用到泰勒公式。让我来分析一下。 若$pi=x_0+d$ 则$x_1=x_0+sin(x_0)$ $=x_0 + sin(pi-d)$ $=x_0+d - 1/6 * d ^3+...$ 故x1比x0更接近于pi，其误差约为 $x_1-pi= x_0+d - 1/6 * d ^3+... -(x_0+d)= 1/6 * d^3$

Ian

樓主的方法其實就是迭代法解方程式 欲解 f(x)=0, 可以考慮$ x_{n+1}=x_n + f(x_n).$ 若${x_n}$收斂至 a, 則 a=a+f(a) ==> f(a)=0. 令 f(x)=sin(x), 則 ${x_n}$收斂至 a=奇數*Pi. 如果 $x_0=3$, 則$ lim x_k = \pi$ 如果 $x_0=7$, 則$ lim x_k = 3\pi$. 類似的, 也可以考慮 $x_{n+1}=x_n + cos(x_n)$, or $x_{n+1}=x_n-tan(x)$ ...etc

荼靡破晓

下好了。

mathematica

N[Pi,10000]