[基础编程]求π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来做快速乘法和除法的话,
本人不胜感激
PiCalculator.rar
45.0k
RAR
15次下载

+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楼
喂还是发表点看法呐~
折叠评论
加载评论中,请稍候...
折叠评论
2010-11-16 11:45:50
caoyuan9642(作者)
3楼
额。。for循环嵌套有什么好的解决办法吗?
我这个算法读写文件不是很经常的,
只是每次迭代结束才写一次文件,是调试时加的

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

谢谢LS几位的指导!


其实我看过FFT的算法,但是那个要涉及到浮点数(复数)的运算,我怀疑速度会大大减慢,[s:246],怎么处理?
折叠评论
加载评论中,请稍候...
折叠评论
2011-05-31 19:48:44
2011-5-31 19:48:44
caoyuan9642(作者)
4楼
话说我现在知道肿么用FFT做乘法了。。
不过自从有了Java...[s:233]
折叠评论
加载评论中,请稍候...
折叠评论

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

caoyuan9642
学者 笔友
文章
49
回复
810
学术分
4
2009/05/06注册,4 年前活动
暂无简介
插入资源
全部
图片
视频
音频
附件
全部
未使用
已使用
正在上传
空空如也~
上传中..{{f.progress}}%
处理中..
上传失败,点击重试
{{f.name}}
空空如也~
(视频){{r.oname}}
{{selectedResourcesId.indexOf(r.rid) + 1}}
插入表情
我的表情
共享表情
Emoji
上传
注意事项
最大尺寸100px,超过会被压缩。为保证效果,建议上传前自行处理。
建议上传自己DIY的表情,严禁上传侵权内容。
点击重试等待上传{{s.progress}}%处理中...已上传
空空如也~
草稿箱
加载中...
此处只插入正文,如果要使用草稿中的其余内容,请点击继续创作。
{{fromNow(d.toc)}}
{{getDraftInfo(d)}}
标题:{{d.t}}
内容:{{d.c}}
继续创作
删除插入插入
{{forum.displayName}}
{{forum.countThreads}}
篇文章,
{{forum.countPosts}}
条回复
{{forum.description || "暂无简介"}}
ID: {{user.uid}}
{{submitted?"":"投诉"}}
请选择违规类型:
{{reason.description}}
支持的图片格式:jpg, jpeg, png