DRM法によるATAN公式π計算FORTRANプログラム
2003/01/10 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
DRM法によるATAN公式π計算ルーチンである。桁数nに対して、演算量はn･log(n)2〜3
に比例する。πを求める公式の2.1 Arctan公式 | を参照。
本計算には多数桁計算部品として
乗除算 | 、
加減算 | 及び
π用I/O | の
ルーチンが必要である。
時間測定用ルーチンと
XCLOCK | して使用している。
このルーチンの仕様はHitachi
SR8000 FORTRANに合わせてある。この部分はシステムに合わせて変更が必要。
2. プログラム
C=================================================================C
C DRM公式によるπ計算プログラム ( ATAN 版 , FORTRAN )
C 2003/01/18 Yasunori Ushiro ( 後 保範 )
C=================================================================C
C Pi Value Computation Program by Atan with DRM Method C
C Pi = M1*Atan(1/X1) + M2*Atan(1/X2) + M3*Atan(1/X3) C
C Atan(1/X) = 1/X - 1/(3*X**3) + 1/(5*X**5) - 1/(7*X**7) +.. C
C Using FMT for High-Precision Multiplication C
C-----------------------------------------------------------------C
C Written by Yasunori Ushiro, 2003/01/18 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), X(0:ND), Y(0:ND), PI(0:ND)
DIMENSION E(0:ND), W(0:4*ND+1), WK(6*ND)
DIMENSION INO(5), IDATA(2,4,5)
C Pi ATAN Constant
DATA INO/2,3,3,4,4/
DATA IDATA/16, 5, -4,239, 0,0, 0,0,
1 48,18, 32, 57, -20,239, 0,0,
2 32,10, -4,239, -16,515, 0,0,
3 176,57, 28,239, -48,682, 96,12943,
4 48,49, 128,57, -20,239, 48,110443/
COMMON /VBASE/BASE,NORID
C Set Parameter (N;FFT Size, NC;Element Number; NN;Decimal Pi)
CALL PIIN(2,ND,ID,N,NC,NN,L10) ! 2;Atan with DRM
C Set Computation Parameter
CALL XCLOCK(T1,5) ! Start Time
ERR = 0.0
CALL SETV(0,PI,N+1) ! PI=float(0)
DO k=1,INO(ID)
CALL XCLOCK(T2,5)
IX = IDATA(2,k,ID) ! Get IX ; ATAN(1/IX)
IB = IDATA(1,k,ID) ! Get IB ; IB*ATAN(1/IX)
XX = IX
ALG = DLOG(10.0D0)/DLOG(XX) ! Terms by a Digits
MLP = 2.0*L10*(NC+4)*ALG/N + 1 ! Number of DRM steps
C Y=ATAN(1/IX)
CALL DRMATN(A,B,X,E,Y,IX,N,MLP,ER1,W,WK)
ERR = MAX(ERR,ER1)
C Pi=Pi+IB*Y
CALL FCMT(Y,IB,X,N) ! X=IB*Y
CALL FADD(PI,X,Y,N,1) ! Y=Pi+X
CALL COPY(Y,PI,N+1) ! Pi=Pi+IB*Y
CALL XCLOCK(T3,5)
WRITE(6,310) IX,MLP,T3-T2
END DO
C Output Pi
CALL PIOUT(PI,NC,L10,T1,ERR)
C
STOP
310 FORMAT(1H ,' ATAN(1/',I6,'), DRM Step=',I3,', Time(s)=',F8.2)
END
C==================================================================C
SUBROUTINE DRMATN(A,B,X,E,Y,IX,NN,MLP,ERR,W,WK)
C==================================================================C
C Compute Y=ATAN(1/IX) by DRM Method C
C Set Y,E,VS and Loop for MLP C
C Compute Y=Y+A/(B*E) C
C------------------------------------------------------------------C
C ARGUMENTS : C
C A(0:NN) R*8, Work, Data A C
C B(0:NN) R*8, Work, Data B C
C X(0:NN) R*8, Work, Data X C
C E(0:NN) R*8, Work, Data E C
C Y(0:NN) R*8, Out, Data Y=ATAN(1/IX) C
C IX I*4, In, Value IX for ATAN(1/IX) C
C NN I*4, In, Number of Data Gloup C
C MLP I*4, In, Number of Loops C
C ERR R*8, Out, Maximum Rounding Error to Convert Int. C
C W(0:4*NN) R*8, Work, Work for DRM C
C WK(6*NN) R*8, Work, Work for Multiplication C
C------------------------------------------------------------------C
C Written by Yasunori Ushiro, 2003/01/18 C
C (Hitachi Ltd. and Waseda University) C
C==================================================================C
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION A(0:NN), B(0:NN), X(0:NN), Y(0:NN)
DIMENSION W(0:4*NN), WK(6*NN)
C Set Y,E
CALL SETV(0,Y,NN+1) ! Y=0
CALL FVSET(1,E,NN) ! E=Float(1)
CALL FCDIV(E,IX,E,NN) ! E=Float(IX)
VS = 1.0
C Main Loop
ERR = 0.0
DO k=1,MLP
C Compute Y=Y+A/(B*E), Compute A,B,X and E=E*X
CALL DRMAT2(A,B,X,E,Y,IX,VS,NN,ER1,W,WK)
ERR = MAX(ERR,ER1)
END DO
C
RETURN
END
C==================================================================C
SUBROUTINE DRMAT2(A,B,X,E,Y,IX,VS,NN,ERR,W,WK)
C==================================================================C
C Outer Routine of DRM Method for ATAN C
C Set Initial Value to A,B,X C
C Compute A,B,X by DRM and E=E*X C
C Y=Y+A/(B*E) by Floating C
C------------------------------------------------------------------C
C ARGUMENTS : C
C A(0:NN) R*8, Work, Data A C
C B(0:NN) R*8, Work, Data B C
C X(0:NN) R*8, Work, Data X C
C E(0:NN) R*8, I/O, Data E C
C Y(0:NN) R*8, I/O, Data Y C
C IX R*8, In, Value IX for ATAN(1/IX) C
C VS R*8, I/O, Start Value of Computation Terms C
C NN I*4, In, Number of Data Gloup C
C ERR R*8, Out, Maximum Rounding Error to Convert Int. C
C W(0:4*NN) R*8, Work, Work for DRM C
C WK(6*NN) R*8, Work, Work for Multiplication C
C------------------------------------------------------------------C
C Written by Yasunori Ushiro, 2003/01/18 C
C (Hitachi Ltd. and Waseda University) C
C==================================================================C
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION A(0:NN), B(0:NN), X(0:NN), Y(0:NN), E(0:NN)
DIMENSION W(0:4*NN), WK(6*NN)
C Set Initial Value to A,B,X
N = 4
CALL SETV(0,A(1),NN) ! Clear A
CALL SETV(0,B(1),NN) ! Clear B
AV = 1.0
DO k=1,NN/N
B(k*N) = VS
CALL NORM(B(k*(N-1)+1),N) ! B(k)=VS+2*(k-1)
A(k*N) = AV ! A(k)=(-1)**(k-1)
VS = VS + 2.0
AV = -AV
END DO
CALL SETV(0,X(1),N) ! Clear X
VX = IX
X(N) = VX*VX
CALL NORM(X(1),N) ! X=IX**2
C Compute A,B,X by DRM
M = NN/(2*N)
CALL LOG2F(2*M,MEP) ! MEP=log2(M)
ERR = 0.0
DO k=1,MEP
CALL DRMAT1(A(1),B(1),X(1),N,M,ER1,W,WK)
ERR = MAX(ERR,ER1)
N = N*2
M = M/2
END DO
C Change Floating A,B,X from Fixed
CALL FLOAT(A,NN)
CALL FLOAT(B,NN)
CALL FLOAT(X,NN)
C E=E*X
CALL FMTCFM(X,E,W,NN,ER1,WK) ! W=E*X
CALL COPY(W,E,NN+1) ! E=E*X
C Y=Y+A/(B*E)
CALL FMTCFM(B,E,W,NN,ER2,WK) ! W=B*E
CALL COPY(W,X,NN+1) ! X=B*E
CALL FDIV(A,X,B,NN,ER3,W,WK) ! B=A/X
CALL FADD(Y,B,W,NN,1) ! W=Y+A/X
CALL COPY(W,Y,NN+1) ! Y=Y+A(B*E)
ERR = MAX(ERR,ER1,ER2,ER3)
C
RETURN
100 FORMAT(5F10.1)
END
C==================================================================C
SUBROUTINE DRMAT1(A,B,X,N,M,ERR,W,WK)
C==================================================================C
C Inner Routine of DRM Method for ATAN C
C W = X(1)*B(2,k) C
C A(1-2,k) = A(1,k)*W + A(2,k)*B(1,k) C
C B(1-2,k) = B(1,k)*B(2,k) for k=1,M C
C X(1-2) = X(1)*X(1) C
C------------------------------------------------------------------C
C ARGUMENTS : C
C A(N,2,M) R*8, I/O, Data A C
C B(N,2,M) R*8, I/O, Data B C
C X(N,2) R*8, I/O, Data X C
C N I*4, In, Number of Computation Elements C
C M I*4, In, Number of Data Gloup C
C ERR R*8, Out, Maximum Rounding Error to Convert Int. C
C W(N,4) R*8, Work, Work for DRM C
C WK(6*N) R*8, Work, Work for Multiplication C
C------------------------------------------------------------------C
C Written by Yasunori Ushiro, 2003/01/18 C
C (Hitachi Ltd. and Waseda University) C
C==================================================================C
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION A(N,2,M), B(N,2,M), X(N,2)
DIMENSION W(N,4), WK(6*N)
C Main Loop for M
ERR = 0.0
DO k=1,M
C W=X*B(2,k)
CALL MTCFM(X,B(1,2,k),W,N,ER1,WK)
IF(W(N,1).NE.0.0D0) THEN
WRITE(6,*) '** Error Stop **'
STOP
END IF
C A(1-2,k)=A(1,k)*W+A(2,k)*B(1,k)
CALL MTCFM(A(1,1,k),W(1,2),W(1,3),N,ER2,WK) ! W(3)=A(1)*W(1)
CALL MTCFM(A(1,2,k),B(1,1,k),W,N,ER3,WK) ! W(1)=A(2)*B(1)
CALL ADD(W,W(1,3),A(1,1,k),2*N,1) ! A=W(1)+W(3)
C B(1-2,k)=B(1,k)*B(2,k)
CALL MTCFM(B(1,1,k),B(1,2,k),W,N,ER4,WK) ! W=B(1)*B(2)
CALL COPY(W,B(1,1,k),2*N) ! B=W
ERR = MAX(ERR,ER1,ER2,ER3,ER4)
END DO
C X(1-2)=X*X
CALL COPY(X,W,N) ! W=X
CALL MTCFM(W,W,X,N,ER1,WK) ! X=W*W
ERR = MAX(ERR,ER1)
C
RETURN
END