http://www.pluto.ai.kyutech.ac.jp/plt/matumoto/pi_small/node9.html
http://www.ecs.cst.nihon-u.ac.jp/~hiroshi/JavaScript/pi.html
{32bit delphi版です。
コンソールアプリとして Dcc32 -CC PiHex.pas で作成します。
フリーパスカルなら PPC386 PiHex.pas です。
}
program PiHex;
uses SysUtils,math;
{ 割り算 ( (16^n)*b/x) mod (2^32) }
function DivM(n,b0,x:Cardinal; f:integer; var fraction:Extended):Cardinal;
var n0:Cardinal;
var a:Cardinal;
var b:Cardinal;
var l:Cardinal;
begin
n0:=n;
a:=0;
b:=b0;
l:=0;
n:=n0*4;
while( (n>=1) and ((x and 1)=0) ) do
begin
dec(n);
x := x shr 1;
end;
while(n>=1) do
begin
a := a shl 1;
b := b shl 1;
while(b>=x) do begin a:=a+1 ; b:=b-x;end;
dec(n);
inc(l);
if(b=b0)then
begin
if(l>=32)then
begin
n:=n mod l;
end;
end;
end;
while(b>=x) do begin a:=a+1 ; b:=b-x;end;
DivM:= a;
if f <0 then fraction := fraction - int(b)/x { 少数点以下 }
else fraction := fraction + int(b)/x; { 少数点以下 }
end;
{π= Σ (4/(8*k+1)-1/(4*k+2)-1/(8*k+5)-1/(8*k+6))/(16^k) }
function Pis(n:integer) : Cardinal;
var s:Cardinal;
var k:integer;
var fr:Extended;
begin
s :=0;
fr:=0;
for k:=0 to n do
begin
s:=s +DivM( n-k,4,8*k+1, 1,fr);
s:=s -DivM( n-k,1,4*k+2,-1,fr);
s:=s -DivM( n-k,1,8*k+5,-1,fr);
s:=s -DivM( n-k,1,8*k+6,-1,fr);
while fr >=1 do {小数点以下の累積結果を小数点以上に反映}
begin
fr := fr - 1;
s := s + 1;
end;
while fr <0 do
begin
fr := fr + 1;
s := s - 1;
end;
end;
Pis:=s;
end;
procedure PisPrint ;
var i,j:integer;
var s1,s2:Cardinal;
begin
writeln('3.');
for i:=1 to 10000000 do
begin
j:=i*4;
s1:=Pis(j+4) shr 16; {４桁ずらして求めて、同じかどうか比較 }
s2:=Pis(j ) and $ffff;
if s1 = s2 then write(IntToHex(s1,4))
else break;
write(' ');
end;
end;
begin
PisPrint;
end.
#include <stdio.h>
/* πを１６進で求める */
/* 割り算 ( (16^n)*b/x) % (2^32) */
unsigned DivM(unsigned n0,unsigned b0,unsigned x,int f,long double *fraction )
{
unsigned a =0;
unsigned b =b0;
unsigned l =0;
unsigned n =n0*4;
while( (n>=1) && ((x & 1)==0) ) { --n; x >>= 1; };
while(n>=1)
{
a <<= 1u;
b <<= 1u;
l++;
n--;
if( b >= x ){ do { ++(a) ; b -= x;} while( b >= x );
if(b==b0)
{
if(l>=32) n = n % l;
}
}
}
if(b>=x) { a=a+b/x ; b=b % x;}
if(f <0 ) (*fraction) -=( (long double) b / (long double)x);
else (*fraction) +=( (long double) b / (long double)x);
return a;
}
/*π= Σ (4/(8*k+1)-1/(4*k+2)-1/(8*k+5)-1/(8*k+6))/(16^k) */
unsigned Pis(int n)
{
int k;
unsigned s;
long double fr;
s =0;
fr=0;
for(k=0 ;k<= n;k++)
{
s=s +DivM( n-k,4,8*k+1, 1,&fr);
s=s -DivM( n-k,1,4*k+2,-1,&fr);
s=s -DivM( n-k,1,8*k+5,-1,&fr);
s=s -DivM( n-k,1,8*k+6,-1,&fr);
while( fr >=1 ) /* 小数点以下の累積結果を小数点以上に反映 */
{
fr = fr - 1;
s = s + 1;
}
while( fr <0 )
{
fr = fr + 1;
s = s - 1;
}
}
return s;
}
void PisPrint()
{
int i,j;
unsigned s1,s2;
printf("%s\n","3.");
for (i=1 ; i<10000000 ;i++)
{
j=i*4;
s1=Pis(j+4) >> 16u; /*４桁ずらして求めて、同じかどうか比較 */
s2=Pis(j ) & 0xffff;
if( s1 == s2 ) printf("%04X ",s1) ;
else break;
}
}
main()
{
PisPrint();
return 0;
}
program pai;
uses SysUtils;
///////////////////////////////
//2000 桁くらのπの計算
//ppc386 -So wpi.pas
//dcc32 -CC wpi.pas
const _N100 = 100; //1000でもOK
type wlong = array[0..1003] of smallint;
/////////////////////////////////////
// 結果の印刷
procedure wprint(var d:wlong ; gd:integer);
var i:integer;
var s:string;
begin
for i:=low(d) to high(d)-gd do
begin
if _N100 = 100 then s:=format('%2d',[d[i]])
else if _N100 = 1000 then s:=format('%3d',[d[i]])
else if _N100 = 10000 then s:=format('%4d',[d[i]]);
if i<>low(d) then
if s[1]=' 'then
begin s[1]:='0'; if s[2]=' 'then
begin s[2]:='0';
if s[3]=' 'then s[3]:='0';
end;
end;
write (s);
if i=0 then write( '.');
end;
writeln;
end;
/////////////////////////////////////
// ゼロクリア
procedure w0(var a:wlong);
var i:integer;
begin
for i:=high(a) downto low(a) do a[i]:=0;
end;
/////////////////////////////////////
// 正規化
procedure wc(var a:wlong);
var r,w:integer;
var i:integer;
begin
r:=0;
for i:=high(a) downto low(a) +1 do
begin
w :=a[i]+r; r:=0;
while(w<0 ) do begin w:=w+_N100; r:=r-1;end;
while(w>_N100) do begin w:=w-_N100; r:=r+1;end;
a[i]:=w;
end;
end;
/////////////////////////////////////
// 割算 a = a / b
function wdiv(var a:wlong; b:integer):boolean;
var r,w:integer;
var i:integer;
begin
r:=0;
result:=true;
for i:=low(a) to high(a) do
begin
w:=a[i]+r;
a[i]:= w div b;
result:=result and ( a[i]=0 ) ;//全てがゼロになったか
r:= w mod b;
r:= r * _N100;
end;
wc(a);
end;
/////////////////////////////////////
// 割算の和 sum =sum + a / b
procedure wdivsum(var sum:wlong;var a:wlong; b:integer);
var r,w:integer;
var i:integer;
begin
r:=0;
for i:=low(a) to high(a) do
begin
w:=a[i]+r;
sum[i]:= sum[i]+(w div b);
r:= w mod b;
r:= r * _N100;
end;
wc(sum);
end;
/////////////////////////////////////
// マチンの公式でπを求める
var a5 :wlong;
var a239:wlong;
var xwpi:wlong;
procedure wwpi;
var k:integer;
begin
w0(xwpi);
w0(a5);
w0(a239);
a5[0]:=16; // 16* arctan(1/5)
k:=1;
repeat
wdiv(a5,5);
wdivsum( xwpi , a5,k);
wdiv(a5,5);
wdiv(a5,5);
k:=k+2;
wdivsum(xwpi, a5 ,-k);
k:=k+2;
until wdiv(a5,5);
k:=1;
a239[0]:=4; // 4* -arctan(1/239)
repeat
wdiv(a239,239); //*x
wdivsum( xwpi , a239,-k);
wdiv(a239,239); //*x
wdiv(a239,239); //*x
k:=k+2;
wdivsum(xwpi, a239 ,k);
k:=k+2;
until wdiv(a239,239); //*x
end;
/////////////////////////////////////
// 計算時間を測定
const kernel32 = 'kernel32.dll';
function GetTickCount:Cardinal; external kernel32 name 'GetTickCount';
var t0:Cardinal;
var t1:Cardinal;
var t2:Cardinal;
begin
t0:=GetTickCount;
wwpi;
t1:=GetTickCount;
wprint(xwpi,3);
t2:=GetTickCount;
writeln('exec time=', t1-t0,'ms');
writeln('print time=', t2-t1,'ms');
writeln('total time=', t2-t0,'ms');
end.
java によるπ
http://www.matsusaka-u.ac.jp/~okumura/java/keisan-big.html
c によるπhttp://www.nara-edu.ac.jp/~asait/visual_c/sample0/pi.htm
c によるπ(プログラム入門編)http://www.first.tsukuba.ac.jp/~yuya/mathematics/pi/index-j.html
C インラインアセンブラhttp://www.osk.3web.ne.jp/~akasaki/