练习编程能力~~


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


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


算法非常经典~——沙-波法计算π
具体过程为：
1.令a₀=1,b₀=1/√2

t₀=1/4
2.按下列顺序迭代：
①a_n+1=(a_n+b_n)/2
②b_n+1=√(a_n·b_n)
③t_n+1=t_n-2^n-1·(a_n-a_n+1)²


并且这个算法是平方收敛的（每次迭代有效位数增加一倍），
这样要求到Q位有效数字，就至少要log₂Q+1次迭代

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

而由这个算法中可以看到，b_n的计算需要开根号。
这里有两种开根号的办法：
①用初中教我们的方法，每次开两位根号，用除法减掉，
这样所需时间为n²·n/2=O(n³),不可忍受！
②用牛顿迭代法。
即
x₀=任意正数，
x_n+1=(x_n+a/x_n)/2
同样，这个算法也是平方收敛的。
其中除法是O(n²)的
故还是需要O(n²log[sub]2n)的时间

显然方法②要好的多

现在来盘点一下：
计算a_n需要一次加法，一次折半（是O(n)的）
计算b_n需要一次乘法(O(n²))，一次根号(O(n²log₂n))
计算t _n需要一次减法(O(n)),一次乘法(O(n²)),一次移位（计算2^n-1,O(1)),一次减法
需要迭代log₂n次

总时间复杂度：
(O(n)+O(n²)+O(n²log₂n)+2*O(n)+O(n²)+O(1))*log₂n
=O(n²log²₂n)

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

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

程序不长，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来做快速乘法和除法的话，
本人不胜感激
PiCalculator.rar 44.01KB RAR 83次下载