Orz日記 by Akio Morita

_ [Fortran]累乗演算子

GCC-4.8の累乗演算子は、非正規化数の取扱いが正しく無いようです 以下のコードの実行結果は、 GNU Fortran (FreeBSD Ports Collection) 4.8.0 20120715 (experimental) 上で

     integer :: i1022 = -1022
     integer :: i1023 = -1023
     integer :: i1024 = -1024
     integer :: i1025 = -1025
     write(*,*)2.d0**(-1022), 2.d0**i1022, 2.d0**( 0)*2.d0**i1022
     write(*,*)2.d0**(-1023), 2.d0**i1023, 2.d0**(-1)*2.d0**i1022
     write(*,*)2.d0**(-1024), 2.d0**i1024, 2.d0**(-2)*2.d0**i1022
     write(*,*)2.d0**(-1025), 2.d0**i1025, 2.d0**(-3)*2.d0**i1022
     end

 2.22507385850720138E-308  2.22507385850720138E-308  2.22507385850720138E-308
 1.11253692925360069E-308  1.11253692925360069E-308  1.11253692925360069E-308
 5.56268464626800346E-309   0.0000000000000000       5.56268464626800346E-309
 2.78134232313400173E-309   0.0000000000000000       2.78134232313400173E-309

コンパイル時の定数畳み込みや正規化数同士の乗算で発生する非正規化数は問 題ないが、累乗演算子で直接発生する非正規化数がアンダーフロー扱い になっている模様

libmのpow関数は、非正規化数を正しく返す模様なので、gfortran 側の問題？

追記

     real(8) :: i1022 = -1022.d0
     real(8) :: i1023 = -1023.d0
     real(8) :: i1024 = -1024.d0
     real(8) :: i1025 = -1025.d0
     write(*,*)2.d0**(-1022), 2.d0**i1022, 2.d0**( 0)*2.d0**i1022
     write(*,*)2.d0**(-1023), 2.d0**i1023, 2.d0**(-1)*2.d0**i1022
     write(*,*)2.d0**(-1024), 2.d0**i1024, 2.d0**(-2)*2.d0**i1022
     write(*,*)2.d0**(-1025), 2.d0**i1025, 2.d0**(-3)*2.d0**i1022
     end

は問題なく、

 2.22507385850720138E-308  2.22507385850720138E-308  2.22507385850720138E-308
 1.11253692925360069E-308  1.11253692925360069E-308  1.11253692925360069E-308
 5.56268464626800346E-309  5.56268464626800346E-309  5.56268464626800346E-309
 2.78134232313400173E-309  2.78134232313400173E-309  2.78134232313400173E-309

を返すので、integerを右オペランドに受ける累乗演算子の実装の問題と確定