简单方法算pi和e，计算到小数点后面一万位
1， 计算pi(圆周率)
算pi的公式非常之多，下面是几个有趣的：
pi/4=1+(1*1/(2+(3*3/(2+(5*5/(2+…)))))) [布朗克连分式]
2/pi=(sqrt(2)/2)*(sqrt(2+sqrt(2))/2)*(sqrt(2+sqrt(2+sqrt(2)))/2)*… [韦达恒等式]
pi/4=(2/1)*(2/3)*(4/3)*(4/5)*(6/5)*(6/7)*(8/7)*(8/9)*… [华里达表达式]
pi*pi/6=1/1*1+1/2*2+1/3*3+… [欧拉等式]
虽然这些连分式、无穷级数、无穷乘积都收敛到pi或和pi相关的数值，但是收敛速度却各不相同，收敛和收敛速度是两个完全不同的概念。
现在PC机上很有名的SuperPi，采用了Gauss-Legendre算法，这个公式每迭代一次将得到双倍的十进制精度，比如要计算100万位，迭代20次就够了。虽然收敛速度是超线性的，但是随着展开的数据位数越来越高，迭代一次的代价会越来越大，这关系到高精度算法，实现起来相当复杂。高精度，这是个值得思考的问题。
在这里采用一个简单的公式计算，Machin公式：pi=16arctg(1/5)-4arctg(1/239) （可以把这个表达式输入到google上，然后google会告诉你结果是3.1415926）
下面从arctg(x)的展开出发，简单讨论一下这个公式。
可以从等比数列的极限开始：设|x|<1，1+x+x^2+x^3+…=1/(1-x)，等式左边是无限项，如果只取有限的n项，结果会是什么呢？这关系到一个余项的的问题，对右边的式子进行泰勒展开，有泰勒余项，于是：1/(1-x)=1+x+x^2+…+x^n+ x^(n+1)/(1-eps)^(n+1) （1），eps在0到x之间，最后的那个余项是进行误差分析的基础，也是决定收敛速度的式子。（1）只展开到有限项，于是对它两边进行代换得到，
1/(1+x^2)=1-x^2+x^4-…+(-1)^n*x^(2n)+ x^(2n+2)/(1+eps*eps)^(n+1)
再逐项积分得到，
arctg(x)=x-x^3/3+x^5/5-…+(-1)^n*x^(2n+1)/(2n+1)+ x^(2n+3)/[(2n+3)(1+eps*eps)^(n+1)] (2)
（2）的余项取绝对值，再取eps=0（保守估计），得到R(x)<=x^(2n+3)/(2n+3)<10^(-M) （3），10^(-M)就是我们需要达到的误差限，即十进制小数点后面M位。
如果选择pi=4arctg(1)来算pi，那R(x)<=1/(2n+3)，速度太慢了，比如，如果要算到M=4位，那保守估计n需要取n=4999，实验一下，如下程序：
int main(int argc, char* argv[])
{
int tag=0;
double s=0;
for(double i=4998.0;i>=0.0;i-=1.0)
{
if(tag)
{
s-=4.0/(i*2+1.0);
tag=0;
}
else
{
s+=4.0/(i*2+1.0);
tag=1;
}
printf("%5lf : %lf\n",i,s);
}
return 0;
}
最后输出的最后一部分：
…
10.000000 : 0.099953
9.000000 : -0.110573
8.000000 : 0.124721
7.000000 : -0.141946
6.000000 : 0.165747
5.000000 : -0.197890
4.000000 : 0.246555
3.000000 : -0.324874
2.000000 : 0.475126
1.000000 : -0.858207
0.000000 : 3.141793
可见收敛速度极慢，为什么第4位不是5？实际上第4位都不准确，这是因为公式是摆动的，最后的数据会在3.1413到3.1417之间，如果你要精确到5位，那就要扩大10倍的运算量，就是说这样的摆动会从5千次一直持续到5万次，然后才会稳定下来，如果要算更高位数，完全无法想像。甚至会觉得这样算比割圆还慢，祖冲之不是算到了第7位是6~7之间吗，他给出的密率是相当准确的，为我国的古代数学家感到自豪。
世界上有段时间是以各个国家计算圆周率的位数标志数学水平的，到了现代有人开始怀疑做这样的事是否有意义？现在的SuperPi可以用来测PC机上浮点数运算的能力，因为它会运用大量的FFT算法，浮点数运算量巨大。而目前的世界记录还是在日本，因为他有世界上性能最强的计算机，记录好像是1.2兆（亿上面一个单位）。以下一段摘自互联网：
“原苏联数学家格拉维夫斯基证明了pi的值即使算到100位已完全没有必要了．他算出，假设有一个球体，它的半径等于地球到天狼星的距离1.32*10^
工程上的计算，现实主义上的计算，都不会要求这样高精度的pi，那计算pi就成了数学家的游戏？成了毫无意义的事吗？不完全是，计算pi的值本身虽然没有意义，但这个过程很有意义，算法的利用，计算机上的实现，展现计算机的性能等。
（3）是余项式，如果arctg(x)中的x表现成1/m的形式，则余项将以线性收敛，这样的改进对效率有大幅的提高，最有名的一个公式是Machin公式，如前面所示，两项的和，对于整体的精度，可以取较小的一个，也就是5，保守误差是1/(5^(2n+3))(2n+3)<1/10^M，也就是25^n约为10^M，即算1项，将得到log(25)= 1.398（以10为底）位的十进制精度，于是要计算到M=10000位，那N要取7154。
然后再思考高精度的问题。
如何做到高精度，计算机内的数据一般是定长位的，比如32位表示一个long整数，如果超出了最大范围怎么办？那就用64位的long long，如果是32位的计算机，那实际上也还是将64位拆成高32位和低32位的形式处理，计算机只能处理定长，如果要高精就必须牺牲速度，很简单的道理，算一个32位整数的加法和算一个3万位整数的加法显然是代价不同的，于是就只需要定义数据结构，和实现基本运算了。
考虑级数表达式，最好的情况是只有加减乘除，数据只有整数，没错Machin公式就是如此。加法减法和乘法不会丢失精度，而除法会丢失精度，比如7/3=2.333…它不可能用一个整数表示，就是说小数点怎么移都不会是整数，而你若取2.33或者2.34都不准确，因为它始终存在余数，余数是丢失精度的原因，所以要保存精度就是保存余数，在这里的除法将产生两个整数（相当于有理数），一个是商一个是余数，形如：
a / b = (d,r)
如果以10为基，那高精度无非就是对这样的除法继续做下去，将r*10恢复精度，再除，直到满足精度为止。想像这里的a和b也是(d,r)的形式，再考查连除：
a/b/c=(d1,r1)/c=d1/c=(d2,r2)
恢复一次：r1*10/b=(d1,r1) -> (r2*10+d1)/c=(d2,r2)
可以证明是正确的，于是可以定义为：
(a,0)/b/c=(d1,r1)/c=(d2,r2)
对于每个数据用(d,r)保存，定义(d1,r1)/(b,0)=(r2*10+d1)/b=(d2,r2)，(d,r)中的r表示前一次除得到的余数，所以在当前的除法中要先恢复余数，再加上被除数。
另外再用一个连除的技巧，如下一序列：
1/5 1/5^3 1/5^5 … (4)
是个等比数列，要计算得到结果，要进行多少次除法呢？如果每项都直接做除法，那就需要n*n量级的除法。连除就是把它们连起来除，前一个的结果除25得到后一个结果，总共只需要n次除法。
对比Machin公式，具体计算过程就是先算一遍(4)，再将得到的序列分别除上1,3,5,…，再将结果正负交替地累加起来，如果取基是10000，那算一遍将可以得到十进制4位有效数，结果存到一个数组，全部算完后再整理一遍，因为有可能有的项为负，整理完就得到结果。
源程序是：
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;
}
2、计算e（自然对数的底）
算e有个形式简单，又收敛得比较快的级数，原形是由e^x的泰勒展开，
e=1+1+1/2+1/3!+…+1/n!+… (5)
由于阶乘跑得比指数快，所以会以超线性收敛，而且形式简单，便于编程。
在计算之前也做个技巧，就是把它们叠乘起来，将(5)改写成：
e=2+1/2(1+1/3(1+1/4(1+…1/(n-1)(1+1/n)…)))
从最里面的括号往外算，共做n次除法和加法得一段结果，运算效率也是O(N*M)，但是由于收敛速度快些，所以N项节省一些，所以算一万位的e比上面的程序快很多。
源程序是：
#ifndef UINT
typedef unsigned int UINT;
#endif
int GetLength(int m)
{
int n=1;
double x=2.0/2.72;
while(m)
{
while(x<1.0)
{
x*=(double)(++n+1);
}
while(x>=1.0 && m)
{
x/=10.0;
--m;
}
}
while(x<1.0)
{
x*=(double)(++n+1);
}
return n;
}
int main(int argc, char* argv[])
{
const UINT base=1000000;
int i,k,m,n;
printf("输入计算e的有效位数n=");
scanf("%d",&n);
m=GetLength(n);
UINT *r=new UINT[m+1];
assert(r!=NULL);
for(i=0;i<=m;++i)
r[i]=1;
//memset(r,1,sizeof(unsigned int)*(m+1));
printf("e=2.\n");
for(k=n;k>0;k-=6)
{
UINT y=0;
for(i=m;i>=2;--i)
{
y=y+r[i]*base;
r[i]=y%i;
y/=i;
}
if(k<6)
{
char str[128];
sprintf(str,"%06d",y%base);
str[k]='\0';
printf("%s\t\n",str);
}
else
{
if(y<base)
printf("%06d\t",y);
else
{
printf("XXXXXX\t");
if(r)
delete[]r;
return 1;
}
}
}
if(r)
delete[]r;
return 0;
}
这里由于和式每一项的数据控制在较低精度，再由余数恢复的方法保存精度，再逐段展开，就可以不用专门设计高精运算库，实现较高精度的计算。如果对性能要求更高，有必要好好学一下高精算法，再用高级的收敛速度更快的公式运算，希望能获得更多的惊喜！
正在计算,请等待
____________________
>>>>>>>>>>>>>>>>>>>
pi=3.
14159265358979323846264338327950288419716939937510582097494459230781640628620899
86280348253421170679821480865132823066470938446095505822317253594081284811174502
84102701938521105559644622948954930381964428810975665933446128475648233786783165
27120190914564856692346034861045432664821339360726024914127372458700660631558817
48815209209628292540917153643678925903600113305305488204665213841469519415116094
33057270365759591953092186117381932611793105118548074462379962749567351885752724
89122793818301194912983367336244065664308602139494639522473719070217986094370277
05392171762931767523846748184676694051320005681271452635608277857713427577896091
73637178721468440901224953430146549585371050792279689258923542019956112129021960
86403441815981362977477130996051870721134999999837297804995105973173281609631859
50244594553469083026425223082533446850352619311881710100031378387528865875332083
81420617177669147303598253490428755468731159562863882353787593751957781857780532
17122680661300192787661119590921642019893809525720106548586327886593615338182796
82303019520353018529689957736225994138912497217752834791315155748572424541506959
50829533116861727855889075098381754637464939319255060400927701671139009848824012
85836160356370766010471018194295559619894676783744944825537977472684710404753464
62080466842590694912933136770289891521047521620569660240580381501935112533824300
35587640247496473263914199272604269922796782354781636009341721641219924586315030
28618297455570674983850549458858692699569092721079750930295532116534498720275596
02364806654991198818347977535663698074265425278625518184175746728909777727938000
81647060016145249192173217214772350141441973568548161361157352552133475741849468……
评论