π計算I/O用FORTRANプログラム
2003/01/10 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
π計算用パラメータ入力と、πの値を出力するFORTRANのプログラム。
従来法方式によるAtan公式を利用する場合以外のπ計算ルーチンから使用される。
2. プログラム
C==================================================================C
SUBROUTINE PIIN(MID,ND,ID,N,NC,NN,L10)
C------------------------------------------------------------------C
C Read and Set for Pi Computation Parameter C
C MID I*4, In, Computing Method C
C MID = 0 ; Original Atan(1/X) C
C 1 ; AGM Method C
C 1 ; Atan(1/X) with the DRM C
C 2 ; Ramanujan with the DRM C
C ND I*4, In, Size of Work Dimension (ND >= N) C
C ID I*4, Out, ID in the Computing Method C
C N I*4, Out, Number of Computing Elements for Pi C
C NC I*4, Out, Number of Output Elements for Pi C
C NN I*4, Out, Computating Decimal-Digits of Pi C
C L10 I*4, Out, Number of Decimal-Digits in a Element C
C------------------------------------------------------------------C
C Written by Y. Ushiro ( Waseda Uni. and Hitachi ) 2003.01.25 C
C==================================================================C
IMPLICIT REAL*8(A-H,O-Z)
COMMON /VBASE/BASE,NORID
DIMENSION NSIZE(0:3)
DATA NSIZE/5,1,5,3/ ! Number of Case in MID
CHARACTER*28 MSG1(0:3)
CHARACTER*52 MSG2(0:3)
DATA MSG1/' Type In ID, (ID=1,2,3,4,5)',
1 ' Type In ID, (ID=1) ',
2 ' Type In ID, (ID=1,2,3,4,5)',
3 ' Type In ID, (ID=1,2,3) '/
DATA MSG2/' 1;Machin, 2;Gauss, 3;Klingen., 4;Stoemer, 5;Takano',
1 ' 1;Gauss-Legendre (Now Only) ',
2 ' 1;Machin, 2;Gauss, 3;Klingen., 4;Stoemer, 5;Takano',
3 ' 1;Chudnovsky, 2;Fast-Ramanujan 3;Free-Sqrt '/
CHARACTER*12 FILE1,FILET(0:3)
DATA FILET/'AtanORG?.txt', 'AGM-2nd?.txt',
1 'AtanDRM?.txt', 'RamnDRM?.txt'/
CHARACTER*18 NAME(0:3)
DATA NAME/'Atan Original ', 'AGM Gauss-Legendre',
1 'DRM with Atan ', 'DRM with Ramanujan'/
C Constant Set
L10 = 5 ! L10=5
BASE = 10.0D0**L10 ! BASE=10**5
NORID = -1 ! Nearest Rounding
NDMY = 4
C Check MID
IF(MID.LT.0 .or. MID.GE.4) MID = 0
C ID Read
WRITE(6,*) MSG1(MID)
WRITE(6,*) MSG2(MID)
READ(5,*) ID ! Read ID (Pi Formulation)
IF(ID.LE.0 .or. ID.GT.NSIZE(MID)) ID = 1
C Pi Decimal-Digits Read
WRITE(6,*) ' Type In N ( N ; Number of Pi )'
READ(5,*) NN ! Input Pi Decimal-Digits
N = 2*((NN+2*L10-1)/(2*L10)) + NDMY
NC = N - NDMY ! Computation Decimal-Digits
NN = NC*L10 ! Reset Input NN
CALL LOG2F(N,NB) ! NB=log2(N)
N = MAX(2**NB,16)
IF(N-NDMY.LT.NC) N = 2*N ! N is Number of Elements
IF(N.GT.ND) THEN
WRITE(6,*) ' ** Error Stop ** N > ND'
STOP
END IF
C Output File Open
FILE1 = FILET(MID)
FILE1(8:8) = CHAR(ID+ICHAR('0'))
OPEN(UNIT=1,FILE=FILE1)
WRITE(6,300) NAME(MID),ID,NN,N
WRITE(1,300) NAME(MID),ID,NN,N
C
RETURN
300 FORMAT(1H ,' ** Pi Computation by ',A18,' ** '/,
1 ' ID=',I3,', Pi-Decimal=',I8,', FMT size=',I7)
END
C==================================================================C
SUBROUTINE PIOUT(PI,NC,L10,T1,ERR)
C------------------------------------------------------------------C
C Output Pi Value on File-1 C
C PI(0:NC+3) R*8, In, Pi Decimal-Value with Floating C
C NC I*4, In, Number of Output Elements for PI C
C L10 I*4, In, Number of Decimals in a Element C
C T1 R*8, In, Start Time C
C ERR R*8, In, Rounding Error C
C------------------------------------------------------------------C
C Written by Y. Ushiro ( Waseda Uni. and Hitachi ) 2003.01.25 C
C==================================================================C
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION PI(0:NC+3)
DIMENSION NPI(16)
COMMON /VBASE/BASE,NORID
C
NORID = 1 ! Positive Rounding
CALL NORM(PI(1),NC+3)
CALL XCLOCK(T4,5) ! End Time
WRITE(6,300) T4-T1,ERR ! Time and Rounding Error
WRITE(1,300) T4-T1,ERR
WRITE(1,*) ' Pi = 3.'
DO k=1,NC,16 ! Output Pi Value
KEND = MIN(16,NC-k+1) ! 80 Decimal in a Line
KETA = (k-1)*L10 + 1
DO 10 i=1,KEND
10 NPI(i) = PI(i+k)
WRITE(1,500) KETA,(NPI(i),i=1,KEND)
END DO
C
RETURN
300 FORMAT(1H ,' Total Time(s)=',F8.2,' Rounding Error=',E10.2)
500 FORMAT(1H ,I8,10(1X,I5.5,I5.5))
END