LoginSignup
6
3

More than 5 years have passed since last update.

@implicit_none

FortranでNaNを定数として宣言する何通りかの方法

Last updated at Posted at 2019-02-11

概要

FortranでNaN (Not a Number)を定数として宣言する方法を調べました．

使用環境

コンパイラ バージョン
intel Parallel Studio XE Composer Edition for Fortran 19.0.0.117
PGI Visual Fortran for Windows 18.7
gfortran 7.3.0 (for Ubuntu on WSL)

更新履歴

2019年2月12日　NaNの説明が間違っていたとの指摘を受けたので修正．誤字修正

背景

プログラム中で，2次元空間中のある2点の勾配を求める必要が生じました．勾配を求める関数を作っている際に，勾配を求めずに関数を抜ける場合があるので，その際に「求めていない」ことを表す値を返そうと思い立ちました．
しかし，勾配は±無限大の範囲を取り得るので，それ以外の値としてはNaNくらいしか思いつきませんでした．

勾配を求めずに関数を抜ける場合にNaNを代入し，かつ呼び出しもとでNaNかどうかを判断するために，Nanを定数として定義しました．その時の試行錯誤をまとめています．

方法

Wikipedia1のNaNの項目によると，IEEE 754で定義される浮動小数点において，指数部のビットが全て1で，かつ仮数部が0以外の任意のビット列のときにNaNになります．つまり，仮数部の少なくとも一つのビットが1のときにNaNになります．IEEE 754で定義される浮動小数点において，符号ビットを除く上位7ビットが1のときにNaNになります．

また，NaNを返す可能性のある演算として，以下の処理があげられていました2．

  • 0/0
  • 負数の平方根
  • 負数の対数
  • |1|より大きい値のアークサイン，アークコサイン

これらをふまえ，以下の方法を試しました．

  1. -1の変換
  2. 指数部全てと仮数部の最上位ビットが1の整数の変換
  3. 0/0
  4. $\sqrt{-1}$
  5. $\log{-1}$
  6. $\sin^{-1}2$
  7. $\cos^{-1}2$
  8. Fortranのieee_value()とieee_quiet_nanによる変換

1,2については，整数をtransfer関数でbinary64に変換します．NaNは符号ビットを除く上位7ビットが1であればよく，その他のビットは何でもかまいません．そして，整数の-1は全ビットが1なので，これを実数として扱えばNaNになる要件を満たしています．
6,7については，引数の絶対値が1より大きければよいので，単純に2を引数にしました．
8については，Fortran 2003から導入された関数とクラスで，次のようにして使います．

関数 機能
ieee_value(値, クラス) IEEEで定義されるクラス型を，第1引数の型とkindで返す

クラス型（type(IEEE_CLASS_TYPE)）として，ieee_arithmeticモジュールでは以下の定数が定義されています．

  • IEEE_SIGNALING_NAN
  • IEEE_QUIET_NAN
  • IEEE_NEGATIVE_INF
  • IEEE_NEGATIVE_NORMAL
  • IEEE_NEGATIVE_DENORMAL
  • IEEE_NEGATIVE_ZERO
  • IEEE_POSITIVE_ZERO
  • IEEE_POSITIVE_DENORMAL
  • IEEE_POSITIVE_NORMAL
  • IEEE_POSITIVE_INF
  • IEEE_OTHER_VALUE

ここでは，IEEE_QUIET_NANを使います．

数値実験

NaNが生成されるかの確認

まずはNaNが生成されるかを確かめるために，print文で結果を表示してみました．

program main
    use,intrinsic :: iso_fortran_env
    use,intrinsic :: ieee_arithmetic
    implicit none

    print *,transfer(-1_int64,0.0_real64)
    print *,transfer(Z'7FF8000000000000',0.0_real64) !0111111111111000000000000000000000000000000000000000000000000000
    print *,0.0_real64/0.0_real64
    print *,sqrt(-1.0_real64)
    print *,log(-1.0_real64)
    print *,asin(2.0_real64)
    print *,acos(2.0_real64)
    print *,ieee_value(0.0_real64, ieee_quiet_nan)
end program main

IntelコンパイラとPGIコンパイラは無事にコンパイルに成功し，実行するとNaNが8個出力されます．計算結果では絶対に見たくない光景ですね．

                     NaN
                     NaN
                     NaN
                     NaN
                     NaN
                     NaN
                     NaN
                     NaN

gfortranでは，コンパイルエラーがでます．

$ gfortran main.f90
main.f90:8:22:

     print *,0.0_real64/0.0_real64
                      1
Error: Division by zero at (1)
main.f90:9:17:

     print *,sqrt(-1.0_real64)
                 1
Error: Argument of SQRT at (1) has a negative value
main.f90:10:16:

     print *,log(-1.0_real64)
                1
Error: Argument of LOG at (1) cannot be less than or equal to zero
main.f90:11:17:

     print *,asin(2.0_real64)
                 1
Error: Argument of ASIN at (1) must be between -1 and 1
main.f90:12:17:

     print *,acos(2.0_real64)
                 1
Error: Argument of ACOS at (1) must be between -1 and 1

そのため，コンパイルエラーの出た5個の演算をコメントアウトして実行しました．

                       NaN
                       NaN
                       NaN

定数の定義

NaNを定数として定義するには，型宣言の際に右辺にNaNを置かなければなりません．上述の8通りの方法を右辺においてパラメータを宣言してみました．すべてコメントアウトされていますが，一つずつコメントを外し，コンパイルされるか試しました．

program main
    use,intrinsic :: iso_fortran_env
    use,intrinsic :: ieee_arithmetic
    implicit none

    !real(real64),parameter :: NaN = transfer(-1_int64,0.0_real64)
    !real(real64),parameter :: NaN = transfer(Z'7FF8000000000000',0.0_real64)
    !real(real64),parameter :: NaN = 0.0_real64/0.0_real64
    !real(real64),parameter :: NaN = sqrt(-1.0_real64)
    !real(real64),parameter :: NaN = log(-1.0_real64)
    !real(real64),parameter :: NaN = asin(2.0_real64)
    !real(real64),parameter :: NaN = acos(2.0_real64)
    !real(real64),parameter :: NaN = ieee_value(0.0_real64, ieee_quiet_nan)
end program main
処理 Intel PGI gfortran
transfer(-1_int64,0.0_real64) ○ ○ ○
transfer(Z'7FF8000000000000',0.0_real64) ○ ○ ○
0.0_real64/0.0_real64 ○3 ○4 ×
sqrt(-1.0_real64) × ○5 ×
log(-1.0_real64) × ○5 ×
asin(2.0_real64) × ○5 ×
acos(2.0_real64) × ○5 ×
ieee_value(0.0_real64, ieee_quiet_nan) × × ×

問題なくコンパイルできたieee_value(0.0_real64, ieee_quiet_nan)は定数の右辺としては使えませんでした．

選択肢は

  • transfer(-1_int64,0.0_real64)
  • transfer(Z'7FF8000000000000',0.0_real64)

のどちらかです．前者が簡単でよいと思います．

一方，正の無限大はZ'7FF0000000000000'，負の無限大はZ'FFF0000000000000'で定義されるので，NaNに加えて$\pm \infty$を定数として定義する場合には，transfer(Z'7FF8000000000000',0.0_real64)を利用すると区別できてよいかと思います．

補足

ieee_arithmeticモジュールでは，値がNaNか，無限大かなどを判別する関数も提供しています．使い方は簡単で，引数として渡した実数がNaNなら.true.，NaNでなければ.false.を返します．

program main
    use,intrinsic :: iso_fortran_env
    use,intrinsic :: ieee_arithmetic
    implicit none

    real(real64),parameter :: NaN = transfer(-1_int64,0.0_real64)

    !print *, isNan(NaN)
    print *, ieee_is_nan(NaN)

end program main

 T

値がNaNかどうかを調査する関数としてisNaN()関数が有名ですが，これはFortran標準ではなく，PGIコンパイラには用意されていません．

補足2

NaNにはquiet NaN とsignaling NaNがあります．これもWikipedia6によると，一般に仮数部の最上位ビットでsignalingか否かを判定する実装が多かったが，その意味が逆（0のときquiet/1のときsignaling，もしくは1のときquiet/0のときsignaling）の場合があったということです．

結論

NaNはともだち
こわくないよ

  1. https://ja.wikipedia.org/wiki/NaN ↩

  2. 他にも$0\times \infty$などもありましたが，$\infty$を定義しないといけないので候補から外しました． ↩

  3. 結果がNaNであるとのRemarkが出る ↩

  4. 結果がNaNであるとのWarningが出る ↩

  5. 不正なオペランドであるとのWarningが出る ↩

  6. https://ja.wikipedia.org/wiki/NaN#エンコード ↩

6
3
3

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
6
3