Fortranの累乗計算で虚数が結果として得られるような場合に注意すること

結論

虚数が出てくるような累乗の計算の場合、計算前の値がたとえ実数であるとしても
計算前の値を虚数(e.g. $ a+0i $)のように書き換えてから累乗の計算を行わないと、
結果として$ NaN+0i $という通常欲しくないような値が得られる

例

• $ \sqrt{-100} $ の計算
  → このとき $ \sqrt{-100} = 0+10i $ という結果が(数学的には)得られるはず
  しかしFortranでは-100を例えばreal(8)の変数$x$として格納している状態で

  print *, x ** 0.5
  

  を計算すると$ NaN+0i $ が結果として得られてしまう

  つまりキャストしてから演算しないと計算前の変数の型の下で計算してしまう

  解決策としては型変換関数を用いるか別の型の変数に代入して計算をするかのどちらかをすればよい

検証

• コンパイラ
  • ifort (IFORT) 19.0.5.281 20190815
  • GNU Fortran (GCC) 4.8.5 20150623 (Red Hat 4.8.5-36)
• プログラム

main.f90

program main
    implicit none
    real(8):: minus
    complex*16:: minus_cmp, ans
    ans = 0; minus = -100; minus_cmp = -100  ! Initialization
    ans = minus_cmp**5.0d-1  ! √-100
    print *, "complex", ans  ! Complex*16を使って計算
    ans = minus**5.0d-1      ! √-100
    print *, "real   ", ans  ! Real*8を使って計算
end program main

• 結果(ifort)

ifort main.f90 && ./a.out
 complex (0.000000000000000E+000,10.0000000000000)
 real    (NaN,0.000000000000000E+000)

• 結果(gfortran)

gfortran main.f90 && ./a.out 
 complex (  6.1232339957367673E-016,  10.000000000000002     )
 real    (                      NaN,  0.0000000000000000     )

一度complex*16にしてから $ \sqrt{-100} $を計算するのとreal*8で計算するのでは結果が異なり、
数学的にはcomplex*16で計算したほうが正しい結果を与えることが確認できた!

感想

虚数が出てくる計算には気を付けよう!!