首页 > 代码库 > 使用Machin公式计算

使用Machin公式计算

使用Machin公式计算,并使用百亿进制+末项位数控制,这里可算出数万位(比最简PI快80倍),源代码约40行,在本网页中。计算公式 PI=16arctg(1/5)-4arctg(1/239),其中arctg(x)=x-x^3/3+x^5/5-x^7/7+x^9/9...令X=x^2并提取公因式得:arctg(x)=x(1-X(1/3-X(1/5-X(1/7-X(…,只需迭代b=1/(2n+1)-b*X,n=N,…,3,2,1,0,最后算b*x即得arctg(x)要想快速计算几十万位或几百万位,应使用C++或汇编语言,要取得几千万位请使用Ramanujan公式,要算几亿位或几十亿位请用几何算术平均值算法//PI计算javascript程序,Machin+百亿进制优化//2006.12 许剑伟 莆田十中function add(a,b,n){ //多精度a对多精度b的相加算法(小学加法)for(var i=n-1,f=0;i>=0;i--){a[i]+=b[i]+f;if(a[i]>=10000000000) a[i]-=10000000000,f=1; else f=0;}}function sub0(a,b,r,n){ //多精度a对多精度b的相减算法(小学减法)for(var i=n-1,f=0;i>=0;i--){r[i]=a[i]-b[i]-f;if(r[i]<0) r[i]+=10000000000,f=1; else f=0;}}function div(a,b,n){ //多精度a与单精度b相除算法(小学除法)for(var i=0,f=0,c;i<n;i++){c=a[i]+f*10000000000;a[i]=Math.floor((c+0.1)/b);f=c%b;}}function dao(a,f,b,n){ //倒数(f/b)a[0]=Math.floor(f/b); f=f%b;for(var i=1,c;i<n;i++){c=f*10000000000;a[i]=Math.floor((c+0.1)/b);f=c%b;}}function set(a,v,n){ for(var i=0;i<n;i++) a[i]=0; a[0]=v; a.length=n;} //给数组置0并给首位置初值v//以下计算圆周率,计算公式:Machin PI=16arctg(1/5)-4arctg(1/239)var a=new Array(),b=new Array(),c=new Array(); //三个工作数组,a存PI,b存arctg,c是临时数组function arctg(k,v,zf,N){//求v*arctg(k),zf表示结果累加到a时的正负号for(var i=Math.round(N*23.1/Math.log(k*k)),n=i,n2;i>=0;i--){n2=Math.round((n-i)*N/n)+1; //末项计算位数控制if(n2>N) n2=N;dao(c,v,2*i+1,n2);div(b,k*k,n2);sub0(c,b,b,n2);}div(b,k,N);if(zf>0) add(a,b,N);else sub0(a,b,a,N);}function pi(N){ //N为计算的位数,本程序所得最后5位可能有错set(a,0,N); set(b,0,N); //PI结果数组及arctg数组,初值为0arctg(5,16,1,N);arctg(239,4,-1,N);for(var i=1;i<N;i++) a[i]=String(10000000000+a[i]).substr(1,10); //补足10位return a.join("");}function js(){ ca.innerHTML="最后5位可能有错:PI="+pi(Nw.value-0+1); } //在网页上输出

链接:http://www.fjptsz.com/xxjs/xjw/rj/112/pi23.htm

使用Machin公式计算