AGM法のπ計算FORTRANプログラム
2003/01/10 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
AGM法を利用したπ計算ルーチンである。桁数nに対して、演算量ははn･log(n)2に
比例する。πを求める公式の3.1 ガウス・ルジャンドル公式 | を参照。
本計算には多数桁計算部品として
乗除算 | 、
加減算 | 及び
π用I/O | の
ルーチンが必要である。
時間測定用ルーチンと
XCLOCK | して使用している。
このルーチンの仕様はHitachi
SR8000 FORTRANに合わせてある。この部分はシステムに合わせて変更が必要。
2. プログラム
C=================================================================C
C AGM公式によるπ計算プログラム(O(N*Log2(N)**2)版 , FORTRAN
C 2003/01/06 Yasunori Ushiro ( 後 保範 )
C=================================================================C
C Pi Value Computation Program by AGM ( Gauss-Legendre Method ) C
C A=1, B=1/SQRT(2), T=1/4, S=1 C
C do k=1,log2(N) C
C W=A, A=(W+B)/2, B=SQRT(W*B) C
C T=T-S*(A-W)**2, S=2*S C
C end do C
C Pi = (A+B)**2/(4*T) C
C Using FMT for High-Precision Multiplication C
C-----------------------------------------------------------------C
C Written by Yasunori Ushiro, 2002/09/08 C
C (Hitachi Ltd. and Waseda University) C
C=================================================================C
IMPLICIT REAL*8(A-H,O-Z)
PARAMETER(ND=270000)
DIMENSION A(0:ND), B(0:ND), S(0:ND), PI(0:ND), W(0:ND)
DIMENSION C(0:ND), FW(0:4*ND+1), WK(6*ND)
COMMON /VBASE/BASE,NORID
C Set Input Parameter
CALL PIIN(1,ND,ID,N,NC,NN,L10)
C Initial Value Set (A=1, B=1/SQRT(2), T=1/4, S=1)
CALL XCLOCK(T1,5) ! Start Time
CALL FVSET(1,A,N) ! A=1
CALL FCSQRT(2,B,N,-1,ERR,FW,WK) ! B=1/SQRT(2)
CALL FCDIV(A,4,PI,N) ! T=1/4
CALL FVSET(1,S,N) ! S=1
C Iteration until Conversion
CALL LOG2F(N*L10,LOOP)
DO k=1,LOOP
C W=A, A=(W+B)/2, B=SQRT(W*B)
C WRITE(1,*) 'k=',k
CALL XCLOCK(T2,5)
CALL COPY(A,W,N+1) ! W=A
CALL FADD(W,B,FW,N,1) ! FW=W+B
CALL FCDIV(FW,2,A,N) ! A=(W+B)/2
CALL FMTCFM(W,B,FW,N,ERR1,WK) ! FW=W*B
CALL COPY(FW,C,N+1) ! C=W*B
CALL FSQRT(C,B,N,1,ERR2,FW,WK) ! B=SQRT(W*B)
C T=T-S*(A-W)**2, S=2*S
CALL FADD(A,W,C,N,-1) ! C=A-W
CALL FMTCFM(C,C,FW,N,ERR3,WK) ! FW=C**2
CALL COPY(FW,C,N+1) ! C=(A-W)**2
CALL FMTCFM(S,C,FW,N,ERR4,WK) ! FW=S*(A-W)**2
CALL COPY(PI,C,N+1)
CALL FADD(C,FW,PI,N,-1) ! T=T-S*(A-W)**2
CALL FCMT(S,2,C,N) ! C=2*S
CALL COPY(C,S,N+1) ! S=2*S
ERR = MAX(ERR,ERR1,ERR2,ERR3,ERR4)
CALL XCLOCK(T3,5)
WRITE(6,310) k,T3-T2
END DO
C Pi = (A+B)**2/(4*T)
CALL FADD(A,B,W,N,1) ! W=(A+B)
CALL FMTCFM(W,W,FW,N,ERR1,WK) ! FW=W**2
CALL COPY(FW,W,N+1) ! W=(A+B)**2
CALL FCMT(PI,4,C,N) ! C=4*T
CALL FDIV(W,C,PI,N,ERR2,FW,WK) ! PI=(A+B)**2/(4*T)
ERR = MAX(ERR,ERR1,ERR2)
C Output
CALL PIOUT(PI,NC,L10,T1,ERR)
C
STOP
300 FORMAT(1H ,' ** Pi Computation Start by AGM (Gauss-Legendre) ** '/
1 ' Number of Pi Decimal-Digits =',I8)
310 FORMAT(1H ,' AGM Step=',I3,' Time(s)=',F8.2)
320 FORMAT(1H ,' Total Time(s)=',F8.2,' Rounding Error=',E10.2)
500 FORMAT(1H ,I8,10(1X,I5.5,I5.5))
END