简单方法算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^12公里，在这个球中装满了微生物，假定球的每1立方毫米中有10^10个微生物，然后把所有微生物排列在一条线上，使每两个相邻微生物的间距重新等于地球到天狼星的距离，那么，拿这个幻想长度来作为圆的直径，取pi的值们确到小数点后100位，可以算出这个巨圆的周长们确到1/1000000毫米以下．法国天文学家阿拉哥曾说过“无休止地追求的精确值，没有丝毫精确意义”．”   工程上的计算，现实主义上的计算，都不会要求这样高精度的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值部分结果] 正在计算,请等待 ____________________ >>>>>>>>>>>>>>>>>>>  pi=3. 14159265358979323846264338327950288419716939937510582097494459230781640628620899 86280348253421170679821480865132823066470938446095505822317253594081284811174502 84102701938521105559644622948954930381964428810975665933446128475648233786783165 27120190914564856692346034861045432664821339360726024914127372458700660631558817 48815209209628292540917153643678925903600113305305488204665213841469519415116094 33057270365759591953092186117381932611793105118548074462379962749567351885752724 89122793818301194912983367336244065664308602139494639522473719070217986094370277 05392171762931767523846748184676694051320005681271452635608277857713427577896091 73637178721468440901224953430146549585371050792279689258923542019956112129021960 86403441815981362977477130996051870721134999999837297804995105973173281609631859 50244594553469083026425223082533446850352619311881710100031378387528865875332083 81420617177669147303598253490428755468731159562863882353787593751957781857780532 17122680661300192787661119590921642019893809525720106548586327886593615338182796 82303019520353018529689957736225994138912497217752834791315155748572424541506959 50829533116861727855889075098381754637464939319255060400927701671139009848824012 85836160356370766010471018194295559619894676783744944825537977472684710404753464 62080466842590694912933136770289891521047521620569660240580381501935112533824300 35587640247496473263914199272604269922796782354781636009341721641219924586315030 28618297455570674983850549458858692699569092721079750930295532116534498720275596 02364806654991198818347977535663698074265425278625518184175746728909777727938000 81647060016145249192173217214772350141441973568548161361157352552133475741849468……  

评论