[基础编程]求π10000位(最简单算法)
练习编程能力~~


以前总是想那帮日本人为啥能算到π后面几万亿位的~


现在发现不是不可能的。。。


算法非常经典~——沙-波法计算π
具体过程为:
1.令a0=1,b0=1/√2

t0=1/4
2.按下列顺序迭代:
①an+1=(an+bn)/2
②bn+1=√(an·bn)
③tn+1=tn-2n-1·(an-an+1)2


并且这个算法是平方收敛的(每次迭代有效位数增加一倍),
这样要求到Q位有效数字,就至少要log2Q+1次迭代

显然,这种计算是要高精度表示的,
我所用的算法就是普通的高精度加减乘除
加减是O(n)的,时间花费很少。
乘除是O(n2)的,主要的时间花费都在这上面。

而由这个算法中可以看到,bn的计算需要开根号。
这里有两种开根号的办法:
①用初中教我们的方法,每次开两位根号,用除法减掉,
这样所需时间为n2·n/2=O(n3),不可忍受!
②用牛顿迭代法。

x0=任意正数,
xn+1=(xn+a/xn)/2
同样,这个算法也是平方收敛的。
其中除法是O(n2)的
故还是需要O(n2log[sub]2n)的时间

显然方法②要好的多

现在来盘点一下:
计算an需要一次加法,一次折半(是O(n)的)
计算bn需要一次乘法(O(n2)),一次根号(O(n2log2n))
计算t n需要一次减法(O(n)),一次乘法(O(n2)),一次移位(计算2n-1,O(1)),一次减法
需要迭代log2n次

总时间复杂度:
(O(n)+O(n2)+O(n2log2n)+2*O(n)+O(n2)+O(1))*log2n
=O(n2log22n)

或者更精确地,计算次数的上界:
(n+n*(n+1)/2+n*(n+1)/2*log2n+2*n+n*(n+1)/2+1)*log2n
=(3n+n^2+n+n*(n+1)/2*log2n+1)*log2n
将n=10000代入,
大概计算次数(循环总次数)为
(30000+10000*10000+10000*10001/2*15+1)
<150亿次循环
假设按照1亿次加法需要0.2s来计算:
则每次循环需要约50次加法的时间(算上求余、乘法,以20倍加法时间计)
即一亿次循环需要10s
这样总时间大概需要10s*150=1500s≈25分钟
以上都是运算时间的推导

而这是最慢、最朴素的算法!
有快速离散傅里叶算法(FFT),使乘法变成O(nlog2n)
相应除法也一样。
如此一来,
总的时间复杂度立即下降到O(nlog32n)
接近线性。
这样一来,再用上好点的计算机(超高速的),算出几万亿位也不是不可能了。

程序不长,200行上下,
但求出来的10000位数据经过mathematica检验,精度达到了9992位,达标
由于初始化时需要用到根号二,所以我首先在一个文件(sqrt.dat,在附件)里放上根号0.5的10万位以供直接读取,省时间

可能有人鄙视用C编的程序
但是C++由于里面的一大堆东西显然会降低速度,
原来用Pascal编的,除法一次都要15秒。。。用C编只用7~8秒。(精度一万位)
贴上程序:

/*
* Extended.c
*
*  Created on: 2010-11-13
*      Author: Administrator
*/

#include&lt;string.h&gt;
#include&lt;stdio.h&gt;
#include&lt;stdlib.h&gt;
#include&lt;math.h&gt;
#include&lt;time.h&gt;

#define min(a,b) (a)&lt;(b)?a:b
#define Q 10000
#define M 10

typedef int Extended[Q + 5];
void set(Extended a, int n) {
    memset(a, 0, 4 * (Q + 5));
}

void add(Extended a, Extended b, int pcs) {
    int i;
    for (i = pcs; i &gt;= 1; i--) {
        a += b;
        a[i - 1] += a / M;
        a %= M;
    }
    a[0] += b[0];
}
void minus(Extended a, Extended b, int pcs) {
    int i;
    for (i = pcs; i &gt;= 0; i--) {
        a -= b;
        if (a &lt; 0)
            a[i - 1]--, a += M;
    }
}
void mulnum(Extended a, int n, int pcs) {
    int i, p = 0;
    for (i = pcs; i &gt;= 0; i--) {
        a *= n;
        a += p;
        p = a / M;
        a %= M;
    }
}
void divnum(Extended a, int n, int pcs) {
    int i, p = 0;
    for (i = 0; i &lt;= pcs; i++) {
        a += p * M;
        p = a % n;
        a /= n;
    }
}

int ge(Extended a, Extended b) {
    int* p;
    for (p = a; a &lt;= p + Q; a++, b++)
        if (*a &gt; *b)
            return 1;
        else if (*a &lt; *b)
            return 0;
    return 1;
}
void divide(Extended a, Extended d, int pcs) {
    Extended c, b;
    int s = 0, i, j;
    memcpy(b, d, 4 * (Q + 5));
    set(c, 0);
    while (s &lt;= Q &amp;&amp; !b[s])
        s++;
    if (s == Q)
        return;//divby0
    for (i = s; i &lt;= pcs; i++) {
        while (ge(a, b)) {
            minus(a, b, pcs);
            c[i - s]++;
        }
        for (j = pcs; j &gt;= 1; j--)
            b[j + 1] = b[j];
        b[1] = b[0] % M;
        b[0] /= M;
    }
    memcpy(a, c, 4 * pcs);
}
void mul(Extended a, Extended b, int pcs) {
    int i, j;
    Extended c;
    set(c, 0);
    printf(&quot;  multiplying ...\n&quot;);
    fflush(stdout);
    for (i = 0; i &lt;= pcs; i++) {
        for (j = 0; j &lt;= pcs; j++)
            if (i + j &lt;= pcs &amp;&amp; a &amp;&amp; b[j]) {
                c[i + j] += a * b[j];
                c[i + j - 1] += c[i + j] / M;
                c[i + j] %= M;
            }

    }
    memcpy(a, c, 4 * pcs);
}
void sqroot(Extended a, int pcs) {
    int i;
    Extended c, d;
    memcpy(c, a, 4 * (Q + 5));
    for (i = 0; i &lt; log2(Q) + 1; i++) {
        memcpy(d, c, 4 * (Q + 5));
        divide(d, a, pcs);
        add(a, d, pcs);
        divnum(a, 2, pcs);
        printf(&quot;  sqrt %d finished\n&quot;, i + 1);
        fflush(stdout);
    }
    printf(&quot;  sqrt finished\n&quot;);
    fflush(stdout);
}

int main() {
    Extended a, b, t, a0, t1, pi;
    long long p;
    FILE *sqrt2 = fopen(&quot;sqrt2.dat&quot;, &quot;r&quot;);
    FILE *f = fopen(&quot;pi.txt&quot;, &quot;w&quot;);
    int i,j;
    int timen=time(0);
    set(a, 0);
    a[0] = 1;
    set(t, 0);
    t[1] = 2;
    t[2] = 5;
    p = 1;
    fscanf(sqrt2, &quot;0.&quot;);
    for (i = 1; i &lt;= Q; i++){
        b=fgetc(sqrt2);
        b-=48;
    }

    for (i = 1; i &lt;= log2(Q)+2; i++) {
        printf(&quot;iterating No.%d...\n&quot;, i);
        fflush(stdout);
        printf(&quot;calculating a %d...\n&quot;, i);
        fflush(stdout);
        memcpy(a0, a, 4 * (Q + 5));
        add(a, b, Q);
        divnum(a, 2, Q);
        printf(&quot;calculating b %d...\n&quot;, i);
        fflush(stdout);
        mul(b, a0, Q);
        sqroot(b, Q);
        printf(&quot;calculating t %d...\n&quot;, i);
        fflush(stdout);
        memcpy(t1, a0, 4 * (Q + 5));
        minus(t1, a, Q);
        mul(t1, t1, Q);
        mulnum(t1, p, Q);
        minus(t, t1, Q);
        p &lt;&lt;= 1;

        printf(&quot;finish iterate %d\n&quot;, i);
        fflush(stdout);

        memcpy(t1, a, 4 * (Q + 5));
        add(t1, b, Q);
        mul(t1, t1, Q);
        divnum(t1, 4, Q);
        divide(t1, t, Q);
        memcpy(pi, t1, 4 * (Q + 5));
        fprintf(f,&quot;result of iterate %d\n&quot;,i+1);
        fprintf(f, &quot;%d.&quot;, pi[0]);
        for (j = 1; j &lt;= Q; j++) {
            fprintf(f, &quot;%d&quot;, pi[j]);
            //if (i % 50 == 0)
            //printf(&quot;\n&quot;);
        }
        fprintf(f,&quot;\n\n&quot;);
    }
    memcpy(t1, a, 4 * (Q + 5));
    add(t1, b, Q);
    mul(t1, t1, Q);
    divnum(t1, 4, Q);
    divide(t1, t, Q);
    memcpy(pi, t1, 4 * (Q + 5));

    fprintf(f,&quot;\n\nfinal result&quot;);
    fprintf(f, &quot;%d.&quot;, pi[0]);
    for (i = 1; i &lt;= Q; i++) {
        fprintf(f, &quot;%d&quot;, pi);
        //if (i % 50 == 0)
        //printf(&quot;\n&quot;);
    }
    printf(&quot;time ellapsed...%d s&quot;,time(0)-timen);
    fclose(f);
    return 0;
}



由于我的本子相当慢。。(512M,1.7GHz),所以运行了将近30分钟才算完。
算出来的数字都输出在pi.txt里面


这只是一个有趣的尝试,但其中涉及的沙-波法应该是求π非常好的一个算法了
还有如果谁能解释下怎么用FFT来做快速乘法和除法的话,
本人不胜感激
+1  学术分    科创论坛   2010-11-15   赞扬这样的作品。
+100  科创币    joyeep   2010-11-15   优秀文章
+100  科创币    小俊   2010-11-16   GuLi YuanChuang
+150  科创币    novakon   2010-11-17   当场落泪
+25  科创币    frival   2013-01-19   佩服,果然是编程高手!
+25  科创币    pl_014   2013-01-19   高质量发帖
+50  科创币    量子隧道   2013-01-21   好帖!
来自 软件综合
 
2010-11-15 18:25:35
caoyuan9642(作者)
1楼
编程过程中遇到的困难:
①除法的处理极易出错,因为要减、移位、处理小数点位置等等
②由于在C中数组的传递使用指针,速度固然很快,但是由于本人的误操作,造成指针乱指。。电脑多次崩溃。
所以我这才明白为啥以前学C的时候总是警告我们要小心使用指针。。
折叠评论
加载评论中,请稍候...
折叠评论
caoyuan9642(作者)
2楼
喂还是发表点看法呐~
折叠评论
加载评论中,请稍候...
折叠评论
3楼
赞扬!望LZ解释为什么精度没有到100%。另外建议LZ在fopen时还是try下吧。
折叠评论
加载评论中,请稍候...
折叠评论
4楼
这个算法的核心部位如果内嵌汇编,我相信还要快

LZ 把计算结果都放在内存里面,而是在最后写到文件中,计算过程估计会快3~5倍吧。

现在都是双核CPU,开辟一个线程来读写文件。会比较高效,不会拖累算法的时间。

看到很多for for 及while嵌套,这些是算法不良习惯,要是能够优化更好。

LZ 寻求更高效的算法,可以参考《算法导论》 《算法设计》 《计算机程序设计艺术》
折叠评论
加载评论中,请稍候...
折叠评论
2010-11-16 00:03:22
5楼
本人对程序和算法设计是外行,离散数学和数值计算方法。。都还给老师了。记得FFT的算法都是在数字信号处理里面学的,当时老师讲得唾沫四溅,但是下面似乎除了前四排都在睡觉。[s:94]遗憾的是我坐在第五排。。

我只是从文章的结构和讲述的思路来看待这篇帖子,觉得写得很好,就怂恿同志们加了分。文章以常用的数学语言,通俗易懂的介绍了一个很专业的话题。并且,通过简短的介绍和设计思路的回顾,让读者能够轻松的了解不同算法是怎样节省机器资源来尽快的达到目的的。文章阐述了设计思路,经过假设、推论、实验设计,最终给出了实验结论,最后还对实验进行了讨论,指明了探索的方向。作者思路清晰,逻辑严谨,理据充分,表达准确,表现出较高的科学素质,因此,这是一篇可以作为范例的好文章。
折叠评论
加载评论中,请稍候...
折叠评论
caoyuan9642(作者)
6楼
额。。for循环嵌套有什么好的解决办法吗?
我这个算法读写文件不是很经常的,
只是每次迭代结束才写一次文件,是调试时加的

精度没达到100%是末尾几位有计算误差积累,不可避免的。

谢谢LS几位的指导!


其实我看过FFT的算法,但是那个要涉及到浮点数(复数)的运算,我怀疑速度会大大减慢,[s:246],怎么处理?
折叠评论
加载评论中,请稍候...
折叠评论
2010-11-18 17:26:34
2010-11-18 17:26:34
7楼
大数乘法用FFT算的话效率肯定可以比常规的方法要高很多,尤其是位数非常高的时候。原理我可以简单说说:

1、两个数相乘,如果对两个乘数调整一下的话,其积可以写成原乘数各位数字序列的卷积形式。
2、时域上的卷积,就是频域上的点乘。
3、剩下的自己发挥。

详细原理推导和代码网上应该有很多。如果乘法位数非常高的话甚至还可以考虑用多线程或CUDA加速。
折叠评论
加载评论中,请稍候...
折叠评论
2010-12-17 23:01:51
2010-12-17 23:01:51
8楼
手写高精度除法的男人
折叠评论
加载评论中,请稍候...
折叠评论
2010-12-23 13:03:28
2010-12-23 13:03:28
9楼
话说π笔算到100位是不是要用1个月?
折叠评论
加载评论中,请稍候...
折叠评论
2011-1-26 20:38:44
2011-1-26 20:38:44
10楼
楼主神牛!!
折叠评论
加载评论中,请稍候...
折叠评论
2011-5-26 01:00:48
2011-5-26 01:00:48
11楼
[s:222]楼主大牛!
折叠评论
加载评论中,请稍候...
折叠评论
2011-5-30 07:46:01
2011-5-30 07:46:01
12楼
引用第8楼fw。于2010-12-17 23:01发表的  :
手写高精度除法的男人

高精度除法难道还可以不用手写?
折叠评论
加载评论中,请稍候...
折叠评论
2011-5-31 19:48:44
2011-5-31 19:48:44
caoyuan9642(作者)
13楼
话说我现在知道肿么用FFT做乘法了。。
不过自从有了Java...[s:233]
折叠评论
加载评论中,请稍候...
折叠评论
2011-6-5 21:59:00
2011-6-5 21:59:00
14楼
靠……常数会大到死啊………………Java………………

高精度还是手写好玩……
折叠评论
加载评论中,请稍候...
折叠评论
2012-11-14 22:10:29
2012-11-14 22:10:29
15楼
哇哈哈, 这个我要试试
-5  科创币    justinpiggy   2012-11-14   挖坟!
+5  科创币    孤独的酒精灯   2012-11-15   
折叠评论
加载评论中,请稍候...
折叠评论
16楼
想当初我在极其无聊的时候,曾经将π值背到过600位,现在已经基本忘光了
+1  科创币    phpskycn   2012-11-15   无聊到极点了……
+1  科创币    找自己   2013-04-06   
折叠评论
加载评论中,请稍候...
折叠评论
2012-11-15 21:10:32
17楼
呵呵,让我想起20年前,刚毕业,实在无聊,在XT上算π,10000位算了好像2个小时,全汇编写的,
折叠评论
加载评论中,请稍候...
折叠评论
18楼
表示我背到150位,现在忘了40位……
曾试过用pi/4=1-1/3 1/5-1/7 ……算pi
用C语言编的,代码不到15行……
效率很低,经过统计,计算到6位,共循环50w次
计算位时到了5000w次……
8位开了1min,显示50亿次……
话说为什么我定义double time=0.0
time=time 0.001
运行后显示time3位后还有一堆数字?
+1  科创币    justinpiggy   2013-01-21   默认都是3位吧,如果要保留1位应该要printf(&quot;%.1&quot;,变量名);吧
折叠评论
加载评论中,请稍候...
折叠评论
19楼
代码好复杂,慢慢研究下!
折叠评论
加载评论中,请稍候...
折叠评论
2013-1-19 15:09:35
2013-1-19 15:09:35
20楼
楼主用的啥编译器啊?我用VC6.0,“long long”类型不支持,改成long,再加了几个头文件可以编译通过,但运行了十几分钟还没出结果,不知道是怎么回事?
+1  科创币    acmilan   2013-01-27   int64,,,
折叠评论
加载评论中,请稍候...
折叠评论
21楼
用ASM应该快不少
折叠评论
加载评论中,请稍候...
折叠评论
22楼
回 20楼(frival) 的帖子
C Free编译成功
你用gcc也可以
这个程序运行本来就要很久……
long long改成long估计精度上差很多……
折叠评论
加载评论中,请稍候...
折叠评论
2013-1-21 06:55:21
2013-1-21 06:55:21
23楼
上百度搜一搜"初音圆周率10000位洗脑歌"...
折叠评论
加载评论中,请稍候...
折叠评论
24楼
默认都是3位吧,如果要保留1位应该要printf("%.1",变量名);吧
折叠评论
加载评论中,请稍候...
折叠评论
25楼
引用第7楼小俊于2010-11-18 17:26发表的  :

1、两个数相乘,如果对两个乘数调整一下的话,其积可以写成原乘数各位数字序列的卷积形式。
.......


我猜是例如12345.678...之类的数可以写成1×10^4+2×10^3+3×10^2+4×10^1+5×10^0+6×10^-1+7×10^-2+8×10^-3+...之类的形式。
然后两个数相乘,可以根据上诉表达,写成卷积式。
如12×34=(1×10+2×1)×(3×10+4×1)=1×10×3×10+2×1×3×10+1×10×4×1+2×1×4×1.这是卷积表示的乘法。
然后根据时域卷积等效于频域相乘,则可以把各个位上的数作为时域数列,FFT变换到时域做点乘。然后把结果再iFFT回时域就是结果了。
不知此猜测对否?
折叠评论
加载评论中,请稍候...
折叠评论

想参与大家的讨论?现在就 登录 或者 注册

{{submitted?"":"投诉"}}
请选择违规类型:
{{reason.description}}
支持的图片格式:jpg, jpeg, png