DRM法によるラマヌジャン型公式π計算FORTRANプログラム
2003/01/10 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
DRM法でのラマヌジャン型π計算ルーチンである。桁数nに対して、演算量はn･log(n)2〜3
に比例する。πを求める公式の2.2 ラマヌジャン型公式 | を参照。
本計算には多数桁計算部品として
乗除算 | 、
加減算 | 及び
π用I/O | の
ルーチンが必要である。
時間測定用ルーチンと
XCLOCK | して使用している。
このルーチンの仕様はHitachi
SR8000 FORTRANに合わせてある。この部分はシステムに合わせて変更が必要。
2. プログラム
C=================================================================C
C DRM公式によるπ計算プログラム ( Ramanujan 版 , FORTRAN )
C 2003/01/18 Yasunori Ushiro ( 後 保範 )
C=================================================================C
C Pi Value Computation Program by Ramanujan with DRM Method C
C Pi = K0*K1*SQRT(K2)/(K3*(Y+K4)) C
C Pi = 0 C
C do i=1,ML C
C Pi = (C + A*Y)/B with A,B,C by DRM Method C
C end do 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), C(0:ND), X(0:ND), PI(0:ND)
DIMENSION W(0:4*ND+1), WK(6*ND)
DIMENSION KD(0:5,3), XX(2,3)
C Pi Ramanujan Constant
C K0 K1 K2 K3 K4 K5
DATA KD/128,10005,10005, 3, 13591409, 5, ! ID=1
1 1, 9801, 2, 4, 1103, 3, ! ID=2
2 1, 3528, 1, 1, 1123, 3/ ! ID=3
DATA XX/3.0,53360.0, 4.0,99.0, 2.0,882.0/
COMMON /VBASE/BASE,NORID
C Set Input Parameter
CALL PIIN(3,ND,ID,N,NC,NN,L10)
C Set MLP
ALG = DLOG(10.0D0)/(XX(1,ID)*DLOG(XX(2,ID)))
MLP = 8.0*L10*(NC+4)*ALG/N + 1
WRITE(6,320) MLP
WRITE(1,320) MLP
CALL XCLOCK(T1,5) ! Start Time
ERR = 0.0
C Main Loop
VV = N/8
DO k=1,MLP
VS = (MLP-k)*VV + 1.0 ! Start Value
CALL XCLOCK(T2,5)
C Computing A,B,C by DRM Method
CALL DRMRMA(A,B,C,ID,VS,N,ER1,W,WK)
C Pi = (C+A*Pi)/B or Pi = C/B (k=1)
IF(k.EQ.1) THEN
CALL FVSET(0,W,N)
W(0) = -N ! W=Floating(0)
ELSE
CALL FMTCFM(A,PI,W,N,ER2,WK) ! W=A*Pi
END IF
CALL FADD(C,W,X,N,1) ! X=C+W
CALL FDIV(X,B,PI,N,ER3,W,WK) ! Pi=X/B
ERR = MAX(ERR,ER1,ER2,ER3)
CALL XCLOCK(T3,5)
WRITE(6,310) k,T3-T2
END DO
C Pi=K0*K1*SQRT(K2)/(K3*(K4+K5*Pi))
K2 = KD(2,ID)
IF(K2.EQ.1) THEN
CALL FVSET(1,X,N) ! X=Floating(1)
ELSE
CALL FCSQRT(K2,X,N,1,ER1,W,WK) ! X=SQRT(K2)
END IF
CALL FCMT(X,KD(0,ID),A,N) ! A=K0*X
CALL FCMT(A,KD(1,ID),B,N) ! B=K1*A
CALL FVSET(KD(4,ID),A,N) ! A=K4
CALL FCMT(PI,KD(5,ID),C,N) ! C=K5*Pi
CALL FADD(A,C,X,N,1) ! X=A+C
CALL FCMT(X,KD(3,ID),C,N) ! C=K3*X
CALL FDIV(B,C,PI,N,ER2,W,WK) ! Pi=B/C
ERR = MAX(ERR,ER1,ER2)
C Output
CALL PIOUT(PI,NC,L10,T1,ERR)
C
STOP
310 FORMAT(1H ,' Ramanujan Step=',I3,' Time(s)=',F8.2)
320 FORMAT(1H ,' DRM Steps=',I3)
END
C==================================================================C
SUBROUTINE DRMRMA(A,B,C,ID,VS,NN,ERR,W,WK)
C==================================================================C
C Outer Routine of DRM Method for Ramanujan C
C Set Initial Value to A,B,C C
C Compute A,B,C by DRM and 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 C(0:NN) R*8, Work, Data X C
C ID I*4, In, Type of Ramanujan (1,2,3) 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), C(0:NN)
DIMENSION W(0:4*NN), WK(6*NN)
DIMENSION IKA(6,3), IKB(12,3), CK(2,3)
C A-Data
DATA IKA/2,-1,6,-5,6,-1, 2,-1,4,-3,4,-1, 2,-1,4,-3,4,-1/
C B-Data
DATA IKB/3,0, 3,0, 1,0, 0,106720, 0,106720, 0,106720,
1 2,0, 1,0, 1,0, 0,39204, 0,198, 0,198,
2 2,0, 1,0, 1,0, 0,3528, 0,3528, 0,1/
C C-Data
DATA CK/545140134.D0,13591409.D0, 26390.D0, 1103.D0,
1 21460.D0, 1123.D0/
C Set Initial Value to A,B,C
N = 8 ! Number of Elements for A,B,C
kk = 1 ! Start Address of A,B,C
SG = 1 ! Sign of C
DO k=1,NN/N
IF(ID.NE.2) SG = -SG ! Change Sign of C
C Set A-Data
VS1 = VS + 1.0
CALL VRMSET(A(kk),VS1,IKA(1,ID),N,3) ! Set A
C Set B-Data
CALL VRMSET(B(kk),VS,IKB(1,ID),N,6) ! Set B
C Set C-Data
CALL SETV(0,C(kk),N)
C(kk+N-1) = VS
CALL NORM(C(kk),N) ! Normalize C
CALL CMLT(C(kk),CK(1,ID),N) ! C=C*IKC(1)
C(kk+N-1) = C(kk+N-1) + CK(2,ID) ! C=C+IKC(2)
CALL NORM(C(kk),N) ! Normalize C
CALL CMLT(C(kk),SG,N) ! C=C*IS
CALL NORM(C(kk),N) ! Normalize C
C Count Up
VS = VS + 1.0
kk = kk + N
END DO
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 DRMRM1(A(1),B(1),C(1),N,M,ER1,W,WK) ! DRM Step
ERR = MAX(ERR,ER1)
N = N*2
M = M/2
END DO
C Change A,B,C from Fixed to Floating-Type
CALL FLOAT(A,NN)
CALL FLOAT(B,NN)
CALL FLOAT(C,NN)
C
RETURN
END
C==================================================================C
SUBROUTINE VRMSET(A,VS,IB,N,NO)
C==================================================================C
C Set Initial Value for A-Data and B-Data C
C ARGUMENTS : C
C A(N) R*8, I/O, Input and Output Data C
C VS R*8, I/O, Start Value of Computation Terms C
C IB(2,LP) I*4, In, Constant Table C
C N I*4, In, Number of Elements in A C
C NO I*4, In, Number of Constant Table 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), IB(2,NO)
C
CALL SETV(0,A,N) ! Clear A
A(N) = 1.0 ! A=1 by Fixed Type
DO k=1,NO
X = IB(1,k)*VS + IB(2,k)
CALL CMLT(A,X,N) ! A=A*X by Fixed Type
END DO
C
RETURN
END
C==================================================================C
SUBROUTINE DRMRM1(A,B,C,N,M,ERR,W,WK)
C==================================================================C
C Inner Routine of DRM Method for Ramanujan C
C A(1-2,k) = A(1,k)*A(2,k) C
C B(1-2,k) = B(1,k)*B(2,k) C
C C(1-2,k) = A(1,k)*C(2,k) + B(2,k)*C(1,k) , for k=1,M 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 C(N,2,M) R*8, I/O, Data C 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), C(N,2,M)
DIMENSION W(N,4), WK(6*N)
C Main Loop for M
ERR = 0.0
DO k=1,M
C C(1-2,k) = A(1,k)*C(2,k) + B(2,k)*C(1,k)
CALL MTCFM(A(1,1,k),C(1,2,k),W(1,1),N,ER1,WK) ! W1=A1*C2
CALL MTCFM(B(1,2,k),C(1,1,k),W(1,3),N,ER2,WK) ! W3=B2*C1
CALL ADD(W(1,1),W(1,3),C(1,1,k),2*N,1) ! C=W1+W3
C A(1-2,k) = A(1,k)*A(2,k)
CALL MTCFM(A(1,1,k),A(1,2,k),W,N,ER3,WK) ! W=A1*A2
CALL COPY(W,A(1,1,k),2*N) ! A=W
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=B1*B2
CALL COPY(W,B(1,1,k),2*N) ! B=W
ERR = MAX(ERR,ER1,ER2,ER3,ER4)
END DO
C
RETURN
END