FORTRAN66 のプログラムを書いたことはあまりないんですが、制限がきつくて思いがけないところでしかられます。
とりあえず、台形公式とシンプソンの公式でを求めるという典型的例題を解いてみました。
IBM の16進浮動小数のためか単精度では意外に精度が出ないので倍精度にしてみました。ラインプリンタ出力にはプリンタページングのプリンタ制御文字が入っているようで黒い豆腐■が、プログラム出力の前後に出ています。
■台形公式
IBM FORTRAN IV H コンパイラ用 JOB カード
//FORTHCLG JOB (001),'FORTRAN H INTEGRAL', // CLASS=A,MSGCLASS=A,MSGLEVEL=(1,1) //INTEGRAL EXEC FORTHCLG,REGION.FORT=384K //FORT.SYSLIN DD UNIT=SYSDA //FORT.SYSABEND DD SYSOUT=A //FORT.SYSIN DD * C C TRAPEZOIDAL INTEGERAL PI C DOUBLE PRECISION X0, X1, X, DX DOUBLE PRECISION SUM, SUM1, PI DOUBLE PRECISION F NSTEP = 50 X0 = 0.0D0 X1 = 1.0D0 DX = (X1 - X0) / NSTEP C SUM0 = F(X0) + F(X1) C SUM1 = 0.0D0 N = NSTEP - 1 DO 10 I = 1, N X = X0 + I * DX SUM1 = SUM1 + F(X) 10 CONTINUE C SUM = DX * ( 0.5D0 * SUM0 + SUM1 ) C PI = 4.0D0 * SUM WRITE(6, 100) PI 100 FORMAT(6H PI = , F20.15) STOP END C********************************** DOUBLE PRECISION FUNCTION F(X) DOUBLE PRECISION X F = 1.0D0 / (1.0D0 + X * X) RETURN END //
実行結果
PI = 3.141525986923249
■シンプソンの公式
IBM FORTRAN IV H コンパイラ用 JOB カード
//FORTHCLG JOB (001),'FORTRAN H INTEGRAL', // CLASS=A,MSGCLASS=A,MSGLEVEL=(1,1) //INTEGRAL EXEC FORTHCLG,REGION.FORT=384K //FORT.SYSLIN DD UNIT=SYSDA //FORT.SYSABEND DD SYSOUT=A //FORT.SYSIN DD * C C SIMPSON INTEGERAL PI C DOUBLE PRECISION X, X0, X1, DX DOUBLE PRECISION SUM, SUM0, SUM1, SUM2, PI DOUBLE PRECISION F NSTEP = 50 X0 = 0.0D0 X1 = 1.0D0 DX = (X1 - X0) / NSTEP C SUM0 = F(X0) + F(X1) C SUM1 = 0.0D0 N = NSTEP - 2 DO 10 I = 2, N, 2 X = X0 + I * DX SUM1 = SUM1 + F(X) 10 CONTINUE C SUM2 = 0.0D0 N = NSTEP - 1 DO 20 I = 1, N, 2 X = X0 + I * DX SUM2 = SUM2 + F(X) 20 CONTINUE C SUM = DX * ( SUM0 + 2.0D0 * SUM1 + 4.0D0 * SUM2 ) / 3.0D0 C PI = 4.0D0 * SUM WRITE(6, 100) PI 100 FORMAT(6H PI = , F20.15) STOP END C********************************** DOUBLE PRECISION FUNCTION F(X) DOUBLE PRECISION X F = 1.0D0 / (1.0D0 + X * X) RETURN END //
実行結果
PI = 3.141592653587252
