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

概要

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を定数として定義しました．その時の試行錯誤をまとめています．

方法

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

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

• 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があります．これもWikipedia⁶によると，一般に仮数部の最上位ビットで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#エンコード ↩