主题:[原创]写了个Machin公式算pi的程序,速度不快,算到小数点后面1万位
Rick0ne
[专家分:1490] 发布于 2007-03-28 21:35:00
公式是
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 楼
Rick0ne [专家分:1490] 发布于 2007-03-29 19:40:00
改了一下,
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 楼
Rick0ne [专家分:1490] 发布于 2007-03-29 20:04:00
一个数学达人告诉我的,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 楼
knate [专家分:570] 发布于 2007-03-30 18:18:00
这个是计算pi的快速收敛法之1,
我见过更快的,
估计算法最快1次可以达到10位左右??!
14 楼
Rick0ne [专家分:1490] 发布于 2007-03-31 11:40:00
哈哈,Gauss公式确实快些,它是计算3项得到2.5个十进制位精度,Machin公式是计算2项得到1.39个十进制位精度。
计算一项得到10位的,也是线性收敛,但并不是O(n),因为表示的数据会越来越长,所以超高精度的计算还是相当地慢~~
15 楼
liangbch [专家分:1270] 发布于 2007-04-02 20:06:00
楼主:
我更喜欢结构化的程序,这个算法并不复杂,但你的代码对大多数人来讲,阅读起来似乎费劲了一点。我无意贬低楼主,如果我和楼主有同样年限的编程经验,我的代码可能更糟。
另外,数组的个数可以更少一点,楼主用了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 楼
liangbch [专家分:1270] 发布于 2007-04-02 20:19:00
我的代码量确实很大,许多函数仍可以优化,使其代码更短。我想知道,我的代码速度会不会很慢,楼主可否测试一下?
17 楼
liangbch [专家分:1270] 发布于 2007-04-02 20:50:00
这里给出一些关于圆周率计算的web 站点,前者是我最早看到关于圆周率计算最详尽的一个web站点,后者集中了几乎所有的计算圆率的公式,你看了以后可能不会再编普通的计算圆周率的程序了。值的一提的是,后者的网页出自一个中学生之手。
http://www.jason314.com
http://enjoy-math.equn.com/Math/PiFormula.html
18 楼
Rick0ne [专家分:1490] 发布于 2007-04-02 22:37:00
这些公式都没给余项~没法编在我看来
后来我改了,只要4个数组,呵
19 楼
liangbch [专家分:1270] 发布于 2007-04-03 12:08:00
我看到一片文章,其中罗列了当前世界上运行在PC上的几个最快的计算周率的程序,如果你感兴趣,可以去看看。在这些程序中,有相当一部分使用了AGM算法。这个算法较简单,迭代速度快,需要用到大数乘法运算。
http://myownlittleworld.com/miscellaneous/computers/pilargetable.html
20 楼
Rick0ne [专家分:1490] 发布于 2007-04-03 13:01:00
感谢楼上的,我正在写大数乘法,那个AGM算法是几何算术平均值法吗~~
我来回复