[基础编程]求π10000位(最简单算法)
caoyuan96422010/11/15软件综合 IP:安徽
练习编程能力~~


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


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


算法非常经典~——沙-波法计算π
具体过程为:
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<string.h>
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>

#define min(a,b) (a)<(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 >= 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 >= 0; i--) {
        a -= b;
        if (a < 0)
            a[i - 1]--, a += M;
    }
}
void mulnum(Extended a, int n, int pcs) {
    int i, p = 0;
    for (i = pcs; i >= 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 <= pcs; i++) {
        a += p * M;
        p = a % n;
        a /= n;
    }
}

int ge(Extended a, Extended b) {
    int* p;
    for (p = a; a <= p + Q; a++, b++)
        if (*a > *b)
            return 1;
        else if (*a < *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 <= Q && !b[s])
        s++;
    if (s == Q)
        return;//divby0
    for (i = s; i <= pcs; i++) {
        while (ge(a, b)) {
            minus(a, b, pcs);
            c[i - s]++;
        }
        for (j = pcs; j >= 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("  multiplying ...\n");
    fflush(stdout);
    for (i = 0; i <= pcs; i++) {
        for (j = 0; j <= pcs; j++)
            if (i + j <= pcs && a && 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 < log2(Q) + 1; i++) {
        memcpy(d, c, 4 * (Q + 5));
        divide(d, a, pcs);
        add(a, d, pcs);
        divnum(a, 2, pcs);
        printf("  sqrt %d finished\n", i + 1);
        fflush(stdout);
    }
    printf("  sqrt finished\n");
    fflush(stdout);
}

int main() {
    Extended a, b, t, a0, t1, pi;
    long long p;
    FILE *sqrt2 = fopen("sqrt2.dat", "r");
    FILE *f = fopen("pi.txt", "w");
    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, "0.");
    for (i = 1; i <= Q; i++){
        b=fgetc(sqrt2);
        b-=48;
    }

    for (i = 1; i <= log2(Q)+2; i++) {
        printf("iterating No.%d...\n", i);
        fflush(stdout);
        printf("calculating a %d...\n", i);
        fflush(stdout);
        memcpy(a0, a, 4 * (Q + 5));
        add(a, b, Q);
        divnum(a, 2, Q);
        printf("calculating b %d...\n", i);
        fflush(stdout);
        mul(b, a0, Q);
        sqroot(b, Q);
        printf("calculating t %d...\n", 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 <<= 1;

        printf("finish iterate %d\n", 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,"result of iterate %d\n",i+1);
        fprintf(f, "%d.", pi[0]);
        for (j = 1; j <= Q; j++) {
            fprintf(f, "%d", pi[j]);
            //if (i % 50 == 0)
            //printf("\n");
        }
        fprintf(f,"\n\n");
    }
    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,"\n\nfinal result");
    fprintf(f, "%d.", pi[0]);
    for (i = 1; i <= Q; i++) {
        fprintf(f, "%d", pi);
        //if (i % 50 == 0)
        //printf("\n");
    }
    printf("time ellapsed...%d s",time(0)-timen);
    fclose(f);
    return 0;
}



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


这只是一个有趣的尝试,但其中涉及的沙-波法应该是求π非常好的一个算法了
还有如果谁能解释下怎么用FFT来做快速乘法和除法的话,
本人不胜感激
attachment icon PiCalculator.rar 44.01KB RAR 76次下载
+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 好帖!
+1  学术分    科创网    2010/11/15 赞扬这样的作品。
来自:计算机科学 / 软件综合
26
 
已屏蔽 原因:{{ notice.reason }}已屏蔽
{{notice.noticeContent}}
~~空空如也
caoyuan9642 作者
13年6个月前 IP:未同步
265658
编程过程中遇到的困难:
①除法的处理极易出错,因为要减、移位、处理小数点位置等等
②由于在C中数组的传递使用指针,速度固然很快,但是由于本人的误操作,造成指针乱指。。电脑多次崩溃。
所以我这才明白为啥以前学C的时候总是警告我们要小心使用指针。。
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
caoyuan9642作者
13年6个月前 IP:未同步
265695
喂还是发表点看法呐~
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
phpskycn
13年6个月前 IP:未同步
265697
赞扬!望LZ解释为什么精度没有到100%。另外建议LZ在fopen时还是try下吧。
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
joyeep
13年6个月前 IP:未同步
265700
这个算法的核心部位如果内嵌汇编,我相信还要快

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

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

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

LZ 寻求更高效的算法,可以参考《算法导论》 《算法设计》 《计算机程序设计艺术》
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
虎哥
13年6个月前 IP:未同步
265703
本人对程序和算法设计是外行,离散数学和数值计算方法。。都还给老师了。记得FFT的算法都是在数字信号处理里面学的,当时老师讲得唾沫四溅,但是下面似乎除了前四排都在睡觉。[s:94]遗憾的是我坐在第五排。。

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

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

谢谢LS几位的指导!


其实我看过FFT的算法,但是那个要涉及到浮点数(复数)的运算,我怀疑速度会大大减慢,[s:246],怎么处理?
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
小俊
13年6个月前 IP:未同步
266011
大数乘法用FFT算的话效率肯定可以比常规的方法要高很多,尤其是位数非常高的时候。原理我可以简单说说:

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

详细原理推导和代码网上应该有很多。如果乘法位数非常高的话甚至还可以考虑用多线程或CUDA加速。
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
fw。
13年5个月前 IP:未同步
269612
手写高精度除法的男人
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
qharryq
13年5个月前 IP:未同步
270319
话说π笔算到100位是不是要用1个月?
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
talenth1
13年4个月前 IP:未同步
276760
楼主神牛!!
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
新奕秀颖
13年0个月前 IP:未同步
296692
[s:222]楼主大牛!
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
ltl
13年0个月前 IP:未同步
297336
引用第8楼fw。于2010-12-17 23:01发表的  :
手写高精度除法的男人

高精度除法难道还可以不用手写?
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
caoyuan9642作者
13年0个月前 IP:未同步
297511
话说我现在知道肿么用FFT做乘法了。。
不过自从有了Java...[s:233]
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
ltl
12年11个月前 IP:未同步
298452
靠……常数会大到死啊………………Java………………

高精度还是手写好玩……
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
nihao_x
11年6个月前 IP:未同步
470448
哇哈哈, 这个我要试试
-5
科创币
justinpiggy
2012-11-14
挖坟!
+5
科创币
孤独的酒精灯
2012-11-15
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
slof
11年6个月前 IP:未同步
470453
想当初我在极其无聊的时候,曾经将π值背到过600位,现在已经基本忘光了
+1
科创币
phpskycn
2012-11-15
无聊到极点了……
+1
科创币
找自己
2013-04-06
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
yaozhiyu
11年6个月前 IP:未同步
470647
呵呵,让我想起20年前,刚毕业,实在无聊,在XT上算π,10000位算了好像2个小时,全汇编写的,
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
相对论万岁
11年6个月前 IP:未同步
470664
表示我背到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;,变量名);吧
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
xiying12571
11年6个月前 IP:未同步
470666
代码好复杂,慢慢研究下!
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
frival
11年4个月前 IP:未同步
489655
楼主用的啥编译器啊?我用VC6.0,“long long”类型不支持,改成long,再加了几个头文件可以编译通过,但运行了十几分钟还没出结果,不知道是怎么回事?
+1
科创币
acmilan
2013-01-27
int64,,,
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
zx-16533
11年4个月前 IP:未同步
489677
用ASM应该快不少
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
相对论万岁
11年4个月前 IP:未同步
489734
回 20楼(frival) 的帖子
C Free编译成功
你用gcc也可以
这个程序运行本来就要很久……
long long改成long估计精度上差很多……
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
冰虹
11年4个月前 IP:未同步
490281
上百度搜一搜"初音圆周率10000位洗脑歌"...
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
justinpiggy
11年4个月前 IP:未同步
490284
默认都是3位吧,如果要保留1位应该要printf("%.1",变量名);吧
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
量子隧道
11年4个月前 IP:未同步
490340
引用第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回时域就是结果了。
不知此猜测对否?
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论

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

所属专业
上级专业
同级专业
caoyuan9642
学者 笔友
文章
49
回复
810
学术分
4
2009/05/06注册,7年10个月前活动
暂无简介
主体类型:个人
所属领域:无
认证方式:邮箱
IP归属地:未同步
文件下载
加载中...
{{errorInfo}}
{{downloadWarning}}
你在 {{downloadTime}} 下载过当前文件。
文件名称:{{resource.defaultFile.name}}
下载次数:{{resource.hits}}
上传用户:{{uploader.username}}
所需积分:{{costScores}},{{holdScores}}下载当前附件免费{{description}}
积分不足,去充值
文件已丢失

当前账号的附件下载数量限制如下:
时段 个数
{{f.startingTime}}点 - {{f.endTime}}点 {{f.fileCount}}
视频暂不能访问,请登录试试
仅供内部学术交流或培训使用,请先保存到本地。本内容不代表科创观点,未经原作者同意,请勿转载。
音频暂不能访问,请登录试试
支持的图片格式:jpg, jpeg, png
插入公式
评论控制
加载中...
文号:{{pid}}
投诉或举报
加载中...
{{tip}}
请选择违规类型:
{{reason.type}}

空空如也

加载中...
详情
详情
推送到专栏从专栏移除
设为匿名取消匿名
查看作者
回复
只看作者
加入收藏取消收藏
收藏
取消收藏
折叠回复
置顶取消置顶
评学术分
鼓励
设为精选取消精选
管理提醒
编辑
通过审核
评论控制
退修或删除
历史版本
违规记录
投诉或举报
加入黑名单移除黑名单
查看IP
{{format('YYYY/MM/DD HH:mm:ss', toc)}}