回 帖 发 新 帖 刷新版面

主题:[原创]写了个Machin公式算pi的程序,速度不快,算到小数点后面1万位

公式是

pi=16arctg(1/5)-4arctg(1/239)

int main(int argc, char* argv[])
{
    int const N=7200;
    int const M=10000;
    int const B=10000;
    int const L=4;
    // Machin公式 计算pi到一万位
    int s[M/L];
    int r1[N]={0},r2[N]={0},d1[N]={0},d2;
    int r3[N]={0},r4[N]={0},d3[N]={0},d4;
    int i,k,t,p=0,mp=M/L/20;

    r1[0]=1;
    r1[1]=3;
    r3[0]=4;

    printf("正在计算,请等待\n____________________\n");
    
    for(k=0;k<M/L;++k)
    {
        t=r1[0]*B;
        d1[0]=t/0x5;
        r1[0]=t%0x5;
        //
        t=r3[0]*B;
        d3[0]=t/0xEF;
        r3[0]=t%0xEF;
        s[k]=d1[0]-d3[0];
        int tag=0;
        for(i=1;i<N;++i)
        {
            t=r1[i]*B+d1[i-1];
            d1[i]=t/0x19;
            r1[i]=t%0x19;
            t=r2[i]*B+d1[i];
            d2=t/(2*i+1);
            r2[i]=t%(2*i+1);
            //
            t=r3[i]*B+d3[i-1];
            d3[i]=t/0xDF21;
            r3[i]=t%0xDF21;
            t=r4[i]*B+d3[i];
            d4=t/(2*i+1);
            r4[i]=t%(2*i+1);
            if(tag)
            {
                s[k]+=(d2-d4);
                tag=0;
            }
            else
            {
                s[k]+=(d4-d2);
                tag=1;
            }
        }
        if(p==mp)
        {
            printf(">");
            p=0;
        }
        else
            p++;
    }
    for(i=M/L-1;i>=0;i--)
    {
        while(s[i]>=B)
        {
            s[i-1]++;
            s[i]-=B;
        }
        while(s[i]<0)
        {
            s[i-1]--;
            s[i]+=B;
        }
    }
    printf("\npi=3.\n");
    for(i=0;i<M/L;++i)
        printf("%04d",s[i]);
    return 0;
}

复杂度是O(M*N),N是计算的项数,M是要精确到的位数,它们的关系要解个不等式。在Release下编译速度快一点,可能的优化是/和%运算的重复。

结果就自己试试,不帖出来了,经检查,发现最后两位不准确,公式本身决定的,前面都抽样检查过,应该没问题。

回复列表 (共32个回复)

31 楼

boxertony :
   你手头有分治乘法的代码吗,能弄一份给我吗?thx!


还看到一个BBP算法的公式,说可以分布式计算,也从公式本身看不明白,去搜源论文看估计我也难看懂,技术细节还是挺重要的~~

32 楼

牛人啊!拿回去研究一下。

我来回复

您尚未登录,请登录后再回复。点此登录或注册