使用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数组,初值为0
arctg(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