回 帖 发 新 帖 刷新版面

主题:[原创]写了个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个回复)

11 楼

改了一下,

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,d2;
    int r3[N]={0},r4[N]={0},d3,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=t/0x5;
        r1[0]=t%0x5;
        //
        t=r3[0]*B;
        d3=t/0xEF;
        r3[0]=t%0xEF;
        s[k]=d1-d3;
        int tag=0;
        for(i=1;i<N;++i)
        {
            t=r1[i]*B+d1;
            d1=t/0x19;
            r1[i]=t%0x19;
            t=r2[i]*B+d1;
            d2=t/(2*i+1);
            r2[i]=t%(2*i+1);
            //
            t=r3[i]*B+d3;
            d3=t/0xDF21;
            r3[i]=t%0xDF21;
            t=r4[i]*B+d3;
            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\t",s[i]);
    return 0;
}


正在改另外一个公式的,速度会快一些。。。

12 楼

一个数学达人告诉我的,Gauss 公式:

pi=48arctan(1/18)+32arctan(1/57)-20arctan(1/239)

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

    r1[0]=12;
    r1[1]=2;
    r3[0]=32;
    r5[0]=20;
    s[0]=-10000;

    printf("正在计算,请等待\n___________________\n");

    for(k=0;k<M/L;++k)
    {
        t=r1[0]*B;
        d1=t/18;
        r1[0]=t%18;
        //
        t=r3[0]*B;
        d3=t/57;
        r3[0]=t%57;
        //
        t=r5[0]*B;
        d5=t/239;
        r5[0]=t%239;

        s[k]+=d1+d3-d5;
        int tag=0;
        for(i=1;i<N;++i)
        {
            t=r1[i]*B+d1;
            d1=t/324;
            r1[i]=t%324;
            t=r2[i]*B+d1;
            d2=t/(2*i+1);
            r2[i]=t%(2*i+1);
            //
            t=r3[i]*B+d3;
            d3=t/3249;
            r3[i]=t%3249;
            t=r4[i]*B+d3;
            d4=t/(2*i+1);
            r4[i]=t%(2*i+1);
            //
            t=r5[i]*B+d5;
            d5=t/57121;
            r5[i]=t%57121;
            t=r6[i]*B+d5;
            d6=t/(2*i+1);
            r6[i]=t%(2*i+1);
            if(tag)
            {
                s[k]+=(d2+d4-d6);
                tag=0;
            }
            else
            {
                s[k]+=(d6-d4-d2);
                tag=1;
            }
        }
        if(p==mp)
        {
            printf(">");
            p=0;
        }
        else
            p++;
    }
    for(i=M/L-1;i>=0;i--)
    {
        if(s[i]>=B)
        {
            s[i-1]+=(s[i]/B);
            s[i]%=B;
        }
        else if(s[i]<0)
        {
            s[i-1]-=(1+(-s[i])/B);
            s[i]=(B-(-s[i])%B);
        }
    }
    printf("\npi=3.\n");
    for(k=0,i=0;i<M/L;++i,++k)
    {
        if(k==250)
        {
            k=0;
            printf("\n");
        }
        printf("%04d\t",s[i]);
    }
    return 0;
}

理论上会快20%左右,不过好像不明显~~

13 楼

这个是计算pi的快速收敛法之1,
我见过更快的,
估计算法最快1次可以达到10位左右??!

14 楼

哈哈,Gauss公式确实快些,它是计算3项得到2.5个十进制位精度,Machin公式是计算2项得到1.39个十进制位精度。

计算一项得到10位的,也是线性收敛,但并不是O(n),因为表示的数据会越来越长,所以超高精度的计算还是相当地慢~~

15 楼

楼主:
  我更喜欢结构化的程序,这个算法并不复杂,但你的代码对大多数人来讲,阅读起来似乎费劲了一点。我无意贬低楼主,如果我和楼主有同样年限的编程经验,我的代码可能更糟。

  另外,数组的个数可以更少一点,楼主用了7个数组,我的算法只用了4个数组,实现了同样的功能。下面贴出我刚写的代码。

  声明,我的观点和我推荐的代码风格仅代表我个人观点,我见过好多写的很复杂(难以阅读)的代码,但这些代码对聪明的人来说也许并不坏。谢绝争论编程风格等问题。

#include "stdlib.h"
#include "stdio.h"
#include "memory.h"

// 2个4位数相乘,积的结果小于 10^8 < 2^32,故定义RADIX =10000
#define RADIX 10000
typedef unsigned long DWORD;

void print(DWORD arr[],int len)
{
    int i;
    printf("%d.",arr[1]);
    for (i=2;i<len;i++)
        printf("%04d",arr[i]);
    printf("\n");
}

//arr1= arr1 +arr2, arr1[i]=arr1[i]+arr2[i],并处理进位,使得arr1[i]<RADIX
void Array_Add(DWORD arr1[],DWORD arr2[],int len)
{
    DWORD *p1=arr1+len-1;
    DWORD *p2=arr2+len-1;
    DWORD carry=0;
    DWORD sum;
    while (p1>=arr1)
    {
        sum=carry + *p1 + *p2;
        
        if (sum>=RADIX)
        {
            carry=1;
            sum -= RADIX;
        }
        else
            carry=0;
        *p1=sum;
        p1--;    p2--;

    }
}

//arr1 =arr1 - arr2;
//arr1= arr1 - arr2, arr1[i]=arr1[i]-arr2[i],并处理借位,使得arr1[i]<RADIX
void Array_Sub(DWORD arr1[],DWORD arr2[],int len)
{
    DWORD *p1=arr1+len-1;
    DWORD *p2=arr2+len-1;
    DWORD carry=0;
    while (p1>=arr1)
    {
        DWORD t=carry +  *p2;
        if ( *p1 < t)  //不够减,
        {
            *p1 = *p1 + RADIX -t; 
            carry=1;//置借位标志
        }
        else
        {
            *p1 -= t;
            carry=0;
        }
    
        p1--;    p2--;
    }
}

//arr = arr * n
void Array_MulDWORD(DWORD arr[],int len,DWORD n)
{
    DWORD *p=arr+len-1;
    DWORD carry,prod;
    carry=0;
    while (p>=arr)
    {
        prod= (*p) * n + carry;
        *p=prod % RADIX;
        carry=prod / RADIX;
        p--;
    }

}

void Array_DivDWORD(DWORD arr[],int len,DWORD n)
{
    DWORD *p=arr;
    DWORD mod,prod;
    mod=0;
    while (p<arr+len)
    {
        prod=mod * RADIX + *p;
        *p= prod / n;
        mod=prod % n;
        p++;
    }
}


// 计算arctan(1/k)
// arr[0],arr[1]:表示整数部分,arr[2],arr[3]...表示小数部分
// 每个元素表示4位10进制数字
// arctan(x) = x - x^3/3 + x^5/5 - x^7/7
void arctgx(DWORD arr1[],DWORD arr2[],DWORD tmpArr[],int len,int k)
{
    DWORD *p;
    DWORD n;
    int i;
    
    memset(arr1,0,sizeof(DWORD)*len);
    arr1[1]=1;
    
    Array_DivDWORD(arr1,len,k);            // arr1= 1/k
    memcpy(arr2,arr1,sizeof(DWORD)*len);// arr2= 1/k;

    i=0;
    p=arr1;
    n=3;

    while (true )
    {
        Array_DivDWORD(p+i,len-i,k*k); 
        
        while (p[i]==0 && p+i <arr1+len)
            i++;
        
        if (i>=len)
            break;

        memcpy(tmpArr,arr1,sizeof(DWORD)*len);
    
        Array_DivDWORD(tmpArr+i,len-i,n); 

        if (n % 4==1)
            Array_Add(arr2,tmpArr,len);
        else
            Array_Sub(arr2,tmpArr,len);

        n+=2;
    }    
}

// pi=16arctg(1/5)-4arctg(1/239)
void calcPI(DWORD arr1[],DWORD arr2[],DWORD tmpArr[],DWORD arr3[],int len)
{
    memset(arr3,0,sizeof(DWORD)*len);
    
    arctgx(arr1,arr2, tmpArr,len,5);    //arr2=arctg(1/5)
    Array_MulDWORD(arr2,len,16);        //arr2= 16* arctg(1/5);

    memcpy(arr3,arr2,sizeof(DWORD)*len); //arr3= 16* arctg(1/5);

    arctgx(arr1,arr2,tmpArr, len,239);    //arr2 = arctg(1/239)
    Array_MulDWORD(arr2,len,4);            //arr2 = 4* arctg(1/239)

    Array_Sub(arr3,arr2,len);   //arr3= 16* arctg(1/5) -4*arctg(1/239),arr3 即为PI
}


int main(int argc, char* argv[])
{
    int digLen=10000; //计算10000位圆周率

    int len= digLen/4 +3;

    DWORD *arr1=new DWORD[ len];
    DWORD *arr2=new DWORD[ len];
    DWORD *arr3=new DWORD[ len];
    DWORD *tArr=new DWORD[ len];

    calcPI(arr1,arr2,tArr,arr3,len);
    print(arr3,len-1);    //最后一个元素是错的,不打印
    
    delete[] arr1;
    delete[] arr2;
    delete[] arr3;
    delete[] tArr;

    return 0;
}


16 楼

我的代码量确实很大,许多函数仍可以优化,使其代码更短。我想知道,我的代码速度会不会很慢,楼主可否测试一下?

17 楼

这里给出一些关于圆周率计算的web 站点,前者是我最早看到关于圆周率计算最详尽的一个web站点,后者集中了几乎所有的计算圆率的公式,你看了以后可能不会再编普通的计算圆周率的程序了。值的一提的是,后者的网页出自一个中学生之手。

http://www.jason314.com
http://enjoy-math.equn.com/Math/PiFormula.html

18 楼

这些公式都没给余项~没法编在我看来

后来我改了,只要4个数组,呵

19 楼

我看到一片文章,其中罗列了当前世界上运行在PC上的几个最快的计算周率的程序,如果你感兴趣,可以去看看。在这些程序中,有相当一部分使用了AGM算法。这个算法较简单,迭代速度快,需要用到大数乘法运算。

http://myownlittleworld.com/miscellaneous/computers/pilargetable.html

20 楼

感谢楼上的,我正在写大数乘法,那个AGM算法是几何算术平均值法吗~~

我来回复

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