各種基底によるπ計算I/O用FORTRANプログラム
2003/02/15 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
各種基底でのπ計算用パラメータ入力と、πの値を出力するFORTRANのプログラム。
従来法方式によるAtan公式を利用する場合以外のπ計算ルーチンから使用される。
2. プログラム
C==================================================================C
SUBROUTINE XPIIN(MID,ND,ID,N,NC,NN,L10,LBASE,ISB)
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 Digits of Pi C
C L10 I*4, Out, Number of Digits in a Element C
C LBASE I*4, Out, Base Value (3=< LBASE =<62) C
C ISB I*4, Out, 0: Char. 0 Start, 10; Char. a Start C
C------------------------------------------------------------------C
C Written by Y. Ushiro ( Waseda Uni. and Hitachi ) 2003.02.15 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 Check MID
IF(MID.LT.0 .or. MID.GE.4) MID = 0
C Base Read
WRITE(6,*) ' Type In Base Value ( from 3 to 62 )'
READ(5,*) LBASE
C ISB Read
WRITE(6,*) ' Type In 0 or 1 (0; 0 start, 1; a start)'
READ(5,*) ISB
IF(ISB .NE. 0) ISB = 1
ISB = ISB*10
IF(LBASE.LT.3 .or. LBASE.GT.62-ISB) LBASE = 26
C Constant Set
B2L = DLOG10(LBASE*1.0D0)
L10 = 5.7/B2L
BASE = LBASE**L10
NORID = -1 ! Nearest Rounding
NDMY = 4
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
WRITE(6,310) LBASE,L10
WRITE(1,310) LBASE,L10
C
RETURN
300 FORMAT(1H ,' ** Pi Computation by ',A18,' ** '/,
1 ' ID=',I3,', Number of Pi=',I8,', FMT size=',I7)
310 FORMAT(1H ,' Base=',I2,', Digit-Number in a Element=',I2)
END
C==================================================================C
SUBROUTINE XPIOUT(PI,NC,L10,LBASE,ISB,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 LBASE I*4, In, Base Value (3=< LBASE =<62) C
C ISB I*4, In, 0: Char. 0 Start, 10; Char. a Start 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.02.15 C
C==================================================================C
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION PI(0:NC+3)
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.'
LP = 60/L10
DO k=1,NC,LP ! Output Pi Value
KEND = MIN(LP,NC-k+1) ! 60 Characters in a Line
KETA = (k-1)*L10 + 1
WRITE(1,500) KETA
DO i=1,KEND
NPI = PI(i+k)
CALL CHNGBS(NPI,L10,LBASE,ISB)
END DO
END DO
C
RETURN
300 FORMAT(1H ,' Total Time(s)=',F8.2,' Rounding Error=',E10.2)
500 FORMAT(/,I8,' ',\)
END
C==========================================================C
SUBROUTINE CHNGBS(IN,NCR,NBASE,IS)
C==========================================================C
IMPLICIT REAL*8(A-H,O-Z)
CHARACTER*1 IOUT(32)
CHARACTER*1 CTB1(0:61)
DATA CTB1/'0','1','2','3','4','5','6','7','8','9'
1 ,'a','b','c','d','e','f','g','h','i','j'
2 ,'k','l','m','n','o','p','q','r','s','t'
3 ,'u','v','w','x','y','z','A','B','C','D'
4 ,'E','F','G','H','I','J','K','K','M','N'
5 ,'O','P','Q','R','S','T','U','V','W','X'
6 ,'Y','Z'/
C Change to Charcter
IV = IN
DO k=NCR,1,-1
IDV = IV/NBASE
IC = IV - IDV*NBASE
IOUT(k) = CTB1(IC+IS)
IV = IDV
END DO
C Write
WRITE(1,100) (IOUT(i),i=1,NCR)
C
RETURN
100 FORMAT(A1,\)
END