Modern Fortran Explained （赤本） Exercise 9-1 の注 その2の巻

Modern Fortran Explained -2018- p.495

昨日に引き続き、MFE（赤）の二次方程式の問題です。
M.Metcalf の FORTRAN77 例と MFE 中の Fortran90 の解答例を実行してみます。

Modern Fortran Explained: Incorporating Fortran 2018 (Numerical Mathematics and Scientific Computation)

FORTRAN77

注意すべき点は、装置番号を * にした標準出力で carriage control を使っていることで、Intel Fortran の場合はコンパイル時のオプション /ccdefault:fortran で有効にすることが出来ます。

なお行頭１文字目を出力装置制御に使う carriage control は Fortran2008 で廃止になりました。

Fortran Primer

Encyclopedia of Physical Science and Technology, Volume 5, Third Edition

ソース・プログラム

      PROGRAM QROOTS
*        Solution of quadratic equation
*         (Based on code in FORTRAN66)
      CHARACTER*6 ANAME
      COMPLEX COMP
*     
*     Read title and number of solutions from terminal
    1 READ(*, *, END = 999) ANAME, N
      WRITE(*, '('' ROOTS OF QUADRATIC EQUATIONS FROM '',A)') ANAME
      DO 21 I = 1, N
         READ(*, *) A, B, C
         WRITE(*, '(A, I2/3(A,F8.2:TR12))')
     +            '0set No. ', I, ' A = ', A, 'B = ', B, 'C = ', C
         IF (A.EQ.0.) THEN
            IF (B.NE.0.) THEN
               WRITE(*, '('' LINEAR'',TR25,A,F10.3)') 'X = ', -C/B
            ELSE
               WRITE(*, '('' NO ROOTS'')')
            END IF
         ELSE
            D = B**2 - 4.*A*C
            IF (D.LT.0.) THEN
               COMP = CMPLX(-B/(2.*A), SQRT(-D)/(2.*A))
               WRITE(*, '('' COMPLEX'',TR21,''R(X1)= '',F10.3,TR11,
     +                ''I(X1)= '',F10.3/T30,''R(X2)= '',F10.3,TR11,
     +                ''I(X2)= '',F10.3)') COMP, CONJG(COMP)
            ELSE
               SQRTD = SQRT(D)
               REAL1 = (-B + SQRTD)/(2.0*A)
               REAL2 = (-B - SQRTD)/(2.0*A)
               WRITE(*, '('' REAL '',TR25,2(A,F10.3:TR13))')       
     +               'X1 = ', REAL1, 'X2 = ', REAL2
            END IF         
         END IF                
   21 CONTINUE
      WRITE(*, '(''0END OF '',A)') ANAME
      GO TO 1
  999 END

出力結果

 J.DOE 4
ROOTS OF QUADRATIC EQUATIONS FROM J.DOE
        2.        4.        1.

set No.  1
A =     2.00            B =     4.00            C =     1.00
REAL                          X1 =     -0.293             X2 =     -1.707
        2.        4.        3.

set No.  2
A =     2.00            B =     4.00            C =     3.00
COMPLEX                     R(X1)=     -1.000           I(X1)=      0.707
                            R(X2)=     -1.000           I(X2)=     -0.707
        0.      5.67    -11.83

set No.  3
A =     0.00            B =     5.67            C =   -11.83
LINEAR                         X =      2.086
    536.28   -275.61      2.11

set No.  4
A =   536.28            B =  -275.61            C =     2.11
REAL                          X1 =      0.506             X2 =      0.008

END OF J.DOE
^Z
続行するには何かキーを押してください . . .

Fortran90

ソース・プログラム

実数解の時に、 のところで桁落ちが小さくなるような符号の組み合わせを選び、解と係数の関係からもう一方の解を求めているのが近代的進化点です。

                x1 = -(b + sign(sqrt(d), b)) / (2.0*a)
                x2 = c / (x1*a)

!  MFE2018  Solution 9.1 p.494 
!
    program qroots
         
        implicit none
        real :: a, b, c, d, x1, x2
        
        read (*, *) a, b, c
        write(*, *) ' a= ', a, ' b= ', b, ' c= ', c
        if (a == 0.0) then 
            if (b /= 0.0) then 
                write(*, *) ' Linear: x = ', -c/b
            else 
                write(*, *) ' No roots!'
            end if
        else
            d = b**2 - 4.0*a*c
            if (d < 0.0) then
                write(*, *) ' Complex: ', -b/(2.0*a), '+-', sqrt(-d)/(2.0*a)
            else
                x1 = -(b + sign(sqrt(d), b)) / (2.0*a)
                x2 = c / (x1*a)
                write(*, *) ' Real roots:', x1, x2
            end if
        end if
    
    end program qroots

出力結果

単発の計算を４回やっています。

        2.        4.        1.
  a=    2.000000      b=    4.000000      c=    1.000000
  Real roots:  -1.707107     -0.2928932

        2.        4.        3.
  a=    2.000000      b=    4.000000      c=    3.000000
  Complex:   -1.000000     +-  0.7071068

        0.      5.67    -11.83
  a=   0.0000000E+00  b=    5.670000      c=   -11.83000
  Linear: x =    2.086420

    536.28   -275.61      2.11
  a=    536.2800      b=   -275.6100      c=    2.110000
  Real roots:  0.5061559      7.7733188E-03