PI的计算算法

上机7：文件操作一、问题：古人计算pi （3.1415926535 89793 23846……），一般是用割圆法。

即用圆的内接或外切正多边形来逼近圆的周长。

大约公元前1200年，中国人计算到小数点后1位；Archimedes 在公元前250年用正96边形得到pi 小数点后3位的精度；公元263年，刘徽用正3072边形计算到小数点后5位 ；公元480年，祖冲之计算到小数点后7位；Ludolph Van Ceulen 在1609年用正262边形得到了35位精度；1706年，John Machin 计算到小数点后100位。

1874年，William Shanks 计算到小数点后707位（527位正确）。

这种基于几何的算法计算量大，速度慢，吃力不讨好。

1973年Guilloud & Bouyer 用 CDC 7600计算机计算到1,001,250位；1989年Kanada & Tamurar 用 HITACHI S-820/80计算机计算到1,073,741,799位；1999年Takahashi & Kanada 用 HITACHI SR8000计算机计算到206,158,430,000位；pi 的最新计算纪录由两位日本人Daisuke Takahashi 和Yasumasa Kanada 所创造。

他们在日本东京大学的IT 中心，以Gauss-Legendre 算法编写程序，利用一台每秒可执行一万亿次浮点运算的超级计算机，从日本时间1999年9月18日19:00:52起，计算了37小时21分04秒，得到了pi 的206,158,430,208(3*236)位十进制精度，之后和他们于1999年6月27日以Borwein 四次迭代式计算了46小时得到的结果相比，发现最后45位小数有差异，因此他们取小数点后206,158,430,000位的p 值为本次计算结果。

这一结果打破了他们于1999年4月创造的68,719,470,000位的世界纪录。

随着数学的发展，数学家们在进行数学研究时有意无意地发现了许多计算pi 的公式。

Machin 公式 115239164arg arctg tg π=- 357211()(1)35721n n x x x x a rctg x x n --=-+-++--这个公式由英国天文学教授John Machin 于1706年发现。

他利用这个公式计算到了100位的pi 值。

Machin 公式每计算一项可以得到1.4位的十进制精度。

因为它的计算过程中被乘数和被除数都不大于长整数，所以可以很容易地在计算机上编程实现。

Bailey-Borwein-Plouffe 算法 014211()1681848586n n n n n n π∞==---++++∑这个公式简称BBP 公式，由David Bailey, Peter Borwein 和Simon Plouffe 于1995年共同发表。

它打破了传统的p 的算法，可以计算p 的任意第n 位，而不用计算前面的n-1位。

这为p 的分布式计算提供了可行性。

因此使用这个公式计算，与计算常数e （2.718282……）的公式非常类似，编程很容易实现。

二、要求：编写程序计算常数π到小数点后任意位，并将数据写入文件中。

参考程序：#include <iostream.h>#include <fstream.h>#include <iomanip.h>#include <time.h>int B=10000; /* Working base */int LB=4; /* Log10(base) */void Initial( int *pi, int *pi1,int *pi2, int *pi3, long n) {for (long i=0;i<n;i++){pi[i]=0; pi1[i]=0; pi2[i]=0; pi3[i]=0; }pi1[0]=4;pi2[0]=2;pi3[0]=1;}// Is the big real x equal to zero ?int IsZero (long n, int * x,int s) {long i;int carry=0, xi;int *y=new int[n];for ( i=0; i<n; i++) { y[i]=x[i];xi= y[i]+carry*B;y[i] = xi/(8*s+1); carry =xi%(8*s+1); }for (i=0; i<n; i++){ if (y[i] ) {delete []y; return 0;}} delete []y;return 1;}// Addition of big reals : x += yvoid Add (long n, int*x, int yy[],int s) {int carry=0, xi;int *y=new int[n];for (long i=0; i<n; i++) { y[i]=yy[i];xi= y[i]+carry*B;y[i] = xi/(8*s+1); carry =xi%(8*s+1); }carry=0;for ( i=n-1; i>=0; i--) {x[i] += y[i]+carry;if (x[i]<B) carry = 0;else { carry = 1; x[i] -= B; }}delete []y;}void Minus1(long n, int*x, int yy[],int s) {int carry=0, xi;int*y=new int[n];for (long i=0; i<n; i++) { y[i]=yy[i];xi= y[i]+carry*B; y[i] = xi/(8*s+4); carry = xi%((8*s+4)); }carry=0;for (i=n-1; i>=0; i--) {x[i]-= y[i]+carry;if (x[i]<0) {carry = 1;x[i]+= B; }else { carry = 0; }}delete []y;}void Minus2(long n, int*x, int yy[],int s) {int carry=0, xi;int *y=new int[n];for (long i=0; i<n; i++) { y[i]=yy[i];xi= y[i]+carry*B;y[i] = xi/(8*s+5); carry = xi%(8*s+5); }carry=0;for ( i=n-1; i>=0; i--) {x[i]-= y[i]+carry;if (x[i]<0) {carry = 1;x[i]+= B; }else { carry = 0; }}delete []y;}void Minus3(long n, int*x, int yy[],int s) {int carry=0,xi;int*y=new int[n];for (long i=0; i<n; i++) { y[i]=yy[i]; xi= y[i]+carry*B;y[i] = xi/(8*s+6); carry = xi%(8*s+6); }carry=0;for (i=n-1; i>=0; i--) {x[i]-=y[i]+carry;if (x[i]<0) {carry = 1;x[i]+= B; }else { carry = 0; }}delete []y;}// Division of the big real x by the integer dvoid Div (long n, int *x) {int carry=0, xi;for (int i=0; i<n; i++) { xi= x[i]+carry*B; x[i] = xi/16; carry = xi%16; }}// Print the big real xvoid Print (long n, int*x) {long i;cout<<x[0]<<".";for (i=1; i<n; i++) {cout<<setfill('0')<<setw(4)<<x[i];if (i%25==0) cout<<"\t"<<i*4<<endl;}cout<<endl; }void Savepi(long n,int *x) //保存PI到文件中{ofstream ofile("d:\\pi.txt");if(!ofile){cout<<"打开文件失败！！！"<<endl;return ;}long i;ofile<<x[0]<<".";for (i=1; i<n; i++) {ofile<<setfill('0')<<setw(4)<<x[i];if (i%25==0) ofile<<" "<<i*4<<endl;}ofile<<endl;ofile.close( );}// Computation of the constant pi example;void main ( ) {long NbDigits;cout<<"Input the number of to calculate:"; cin>>NbDigits;long int size=3+NbDigits/LB;int*pi = new int[size];int*pi1 = new int[size];int*pi2 = new int[size];int*pi3 = new int[size];long k=0;time_t T1,T2; T1=time(NULL);Initial(pi,pi1,pi2,pi3,size);while ( !IsZero (size, pi1,k) ){Add (size, pi, pi1,k);Minus1(size, pi, pi2,k);Minus2(size, pi, pi3,k);Minus3(size, pi, pi3,k);k++;Div (size, pi1); Div (size, pi2); Div (size, pi3);}Print (size-2, pi); /* Print out of e */]Savepi(size-2,pi);cout<<"k="<<k<<endl;T2=time(NULL);cout<<endl<<"Calculate time: "<<difftime(T2,T1)<<" s"<<endl; delete []pi;delete []pi1;delete []pi2;delete []pi3;}运行结果：Input the number of to calculate: 10003.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348 253421170679 100 821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489 5493038196 200 442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726 024******* 300 724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951 9415116094 400 330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381 8301194912 500 983367336244065664308602139494639522473719070217986094370277053921717629317675238467481846 7669405132 600 000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279 6892589235 700 420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328 1609631859 800 502445945534690830264252230825334468503526193118817101000313783875288658753320838142061717 7669147303 900 598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909 2164201989 1000同时在D盘里生成了一个文件，文件名字为“ pi.txt ”。