使用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