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

[^nan]: https://ja.wikipedia.org/wiki/NaN

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

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

[^nanop]: 他にも$0\times \infty$などもありましたが，$\infty$を定義しないといけないので候補から外しました．

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

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`文で結果を表示してみました．

```Fortran
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個出力されます．計算結果では絶対に見たくない光景ですね．

```console
                     NaN
                     NaN
                     NaN
                     NaN
                     NaN
                     NaN
                     NaN
                     NaN
```

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

```console
$ 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個の演算をコメントアウトして実行しました．

```console
                       NaN
                       NaN
                       NaN
```

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

```Fortran
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`|○[^rem]|○[^warn]|×|
|`sqrt(-1.0_real64)`|×|○[^oprnd]|×|
|`log(-1.0_real64)`|×|○[^oprnd]|×|
|`asin(2.0_real64)`|×|○[^oprnd]|×|
|`acos(2.0_real64)`|×|○[^oprnd]|×|
|`ieee_value(0.0_real64, ieee_quiet_nan)`|×|×|×|

[^rem]: 結果がNaNであるとのRemarkが出る
[^warn]: 結果がNaNであるとのWarningが出る
[^oprnd]: 不正なオペランドであるとのWarningが出る

問題なくコンパイルできた`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.`を返します．

```Fortran
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
```

```console
 T
```

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

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

[^wikienc]: https://ja.wikipedia.org/wiki/NaN#エンコード

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