Fortranで実数の等値比較を手軽かつ確実に行う方法

# 概要
Fortranで実数（浮動小数点数）の等値比較をする場合は，`ieee_next_after`もしくは`ieee_next_up`および`ieee_next_down`を用いて範囲の上下限値を設定し，その中に値が入っているかを確認すると手軽かつ確実です．

# 環境
- Windows 11 64 bit
- gfortran 11.2
- Intel Fortran classic 2021.4.0

# 注意
この記事では，浮動小数点数を実数と書きます．

# 実数の等値比較
Fortranでは，実数の等値比較は`==`で実行できます．
```Fortran
use, intrinsic :: iso_fortran_env
implicit none
real(real32) :: f, g

f = 2.8311e-1
g = 2.8311e-1
print *, f == g ! T（True）が表示される
```

しかし，プログラムを学習していく過程で，実数の比較は`==`で行ってはいけないと学習します．`==`は，上記例のように二つの実数がビットレベルで等しければ真を返しますが，計算機上では実数は計算過程で丸め誤差が入るため，演算した値の等値性を`==`で比較すると，多くの場合等しくないという結果が得られます．

そこで，`==`による等値比較ではなく，二つの実数の差がある範囲に収まっていたら等しいとする技法を用いたりします．

```Fortran
use, intrinsic :: iso_fortran_env
implicit none
real(real32) :: f, g, tol

f = 0.283109874
g = 0.283109993
tol = epsilon(f) ! 1.19209290E-07程度
print *, (abs(g - f) <= tol) ! T（True）が表示される
```

このとき，二つの実数の差をどう決めるかという問題が生じます．二つの実数の差として，マシンイプシロンが用いられることが多いように思います．マシンイプシロンは，
```math
1.0+\varepsilon > 1.0
```
となる最小の$\varepsilon$のことを指します．しかし，マシンイプシロンは，定義によれば数値のオーダーが1付近でのみ用いるべきであり，`f`や`g`が大きい/小さい場合には，想定した通りの比較結果にならない場合があります．

```Fortran
f = 3.40282266e+38
g = 3.40282347e+38
tol = epsilon(f) ! 1.19209290E-07程度
print *, (abs(g - f) <= tol) ! Fが表示される
```

また，極端に大きい/小さい`f`や`g`を想定すると，`g - f`がオーバー/アンダーフローすることがあり，手軽に比較を行えません．

そこで，ULPの概念を導入することにより，このような煩わしさと無関係に実数の等値比較を行います．

# ULP
ULP（Unit in the Last PlaceあるいはUnit of Least Precision）は，実数の精度などを議論する際に用いられる概念であり，ある実数と，その実数のビット表現を1だけ変化（増加）させたときに得られる数との差がULPとなります．

例で用いた数値で試してみます．
```Fortran
real(real32) :: f, g
integer(int32) :: mold_int32
    
f = 1.0
g = 1.0 + epsilon(f)
print '(3E16.8e2)', f, g, h
print '(B32.32)', transfer(f, mold_int32) ! 00111111100000000000000000000000
print '(B32.32)', transfer(g, mold_int32) ! 00111111100000000000000000000001
print '(E16.8e2)', g - f ! 0.11920929E-06
```
確かに1.0にマシンイプシロンを足すと，ビット表現では1だけ増加していることが分かります．定義通り，マシンイプシロンとは，1.0における1 ULPであることが実質的にも確認できます．

```Fortran
real(real32) :: f, g
integer(int32) :: mold_int32
    
f = real(B"00111110100100001111001111000111", kind=real32)
g = real(B"00111110100100001111001111001000", kind=real32)
print '(3E16.8e2)', f, g, g - f ! 0.28310987E+00  0.28310990E+00  0.29802322E-07
```
一方で，0.28310987あたりで1 ULP離れた値を計算し，その差を求めると，0.29802322E-07となり，マシンイプシロンよりかなり小さいことが分かります．

```Fortran
real(real32) :: f, g
integer(int32) :: mold_int32
    
f = real(B"01111111011111111111111111111110", kind=real32)
g = real(B"01111111011111111111111111111111", kind=real32)
print '(3E16.8e2)', f, g, g - f ! 0.34028233E+39  0.34028235E+39  0.20282410E+32
```
同様に，最大値付近では，1 ULPは0.20282410E+32となります．

色々と例を見ましたが，ここで言いたいことは，二つの実数の差を決めるときに，このULPの概念を使えば良いということです．例えば，二つの実数の差が1 ULP以内であれば等しいと見なす，等です．これはパラメータあるいはプリプロセッサマクロとして定義しておくと便利でしょう．

では肝心の1 ULP（あるいはN ULP）離れた値をどう求めるかということですが，それには`ieee_next_after`あるいは`ieee_next_up`, `ieee_next_down`関数を用います．

# `ieee_next_*`
`ieee_next_after`関数は，`ieee_arithmetic`モジュールで定義されており，二つの引数`x, y`を取り，`x`から`y`方向に1 ULP離れた値を返します．`x`は実数（real32, real64, real128）です．`y`には，基本的に`huge()`もしくは`-huge()`を指定することになるでしょう．

```Fortran
use, intrinsic :: iso_fortran_env
use, intrinsic :: ieee_arithmetic
real(real32) :: f, g, h

f = 0.28310987
g = ieee_next_after(f, huge(f))
h = ieee_next_after(f, -huge(f))
print '(3E16.8e2)', f, g, h ! 0.28310987E+00  0.28310990E+00  0.28310984E+00
print '(B32.32)', transfer(f, mold_int32) ! 00111110100100001111001111000111
print '(B32.32)', transfer(g, mold_int32) ! 00111110100100001111001111001000
print '(B32.32)', transfer(h, mold_int32) ! 00111110100100001111001111000110
```

上の例では，`f`から正の最大値あるいは負の最大値方向に 1 ULP離れた値を求めています．ほとんどの場合，`ieee_next_after(x, huge(x))`は`ieee_next_up(x)`で，`ieee_next_after(x, -huge(x))`は`ieee_next_down(x)`で代用できます．

`ieee_next_after`のうれしいところは，実行の過程でオーバー/アンダーフローを発生させないようにできるところです．下の例では，`g`を32 bit実数の正の最大値とし，`f`は負の最大値方向に1 ULP離れた値，`h`は正の最大値方向に1 ULP離れた値を求めています．

```Fortran
use, intrinsic :: iso_fortran_env
use, intrinsic :: ieee_arithmetic
real(real32) :: f, g, h

g = huge(g)
f = ieee_next_after(g, -huge(g))
h = ieee_next_after(g, huge(g))
print '(3E16.8e2)', f, g, h ! 0.34028233E+39  0.34028235E+39  0.34028235E+39
print '(B32.32)', transfer(f, mold_int32) ! 01111111011111111111111111111110
print '(B32.32)', transfer(g, mold_int32) ! 01111111011111111111111111111111
print '(B32.32)', transfer(h, mold_int32) ! 01111111011111111111111111111111
```

`h`のビット表現を見ると`g`と同じであり，オーバーフローして負の最大値になっていません．これは，二つの引数`x, y`が等しいため，`x`が正の最大値の場合は値がそのまま返されるためです．一方，`ieee_next_up`を使うと，`h`はオーバーフローのシグナルなしで`Infinity`を返します．



# 実数の比較
`ieee_next_after`を用いて等値比較をするオペレータを定義してみます．

まず，二つの実数が等しいと判定する差のULPを`WITHIN_ULP`というパラメータとして定義します．
次に，基準となる値から正負の最大値方向に`WITHIN_ULP`だけ離れた実数を計算するための関数を，`get_lower_bound`および`get_upper_bound`として定義します．それぞれの関数の中では，`do`文を用いて1 ULPずつ離れた値を取得し，戻り値の値を更新していきます．もし，`ieee_next_after`によって得られた値と基準となる値が等しければ，負の最大値もしくは正の最大値であると判断し，離れた値を求める必要はないとして，関数から抜けるようにしています．ただし，これは必須の処理ではなく，`ieee_next_after`の負荷と`if`文の負荷を考慮して，軽い方を選べばよいと思います．

```Fortran
        lower_bound = f
        do i = 1, ulp
            lower_bound = ieee_next_after(lower_bound, -huge(f))
            if (lower_bound == f) return
            ! lower_boundとfが等しければfは負の最大値なので，
            ! 離れた値をこれ以上探す必要はない．
        end do
```

次に，二つの値の等値比較を行う関数を定義します．基準値を`f`，近似値を`g`とします．`g`が基準値`f`から±`WITHIN_ULP`だけ離れた数の範囲に入っていれば，`f`と`g`が等しいとしています．

```Fortran
approx_r32 = (get_lower_bound(f, WITHIN_ULP) <= g .and. g <= get_upper_bound(f, WITHIN_ULP))
```

この関数は，`interface`文によって`.approx.`という演算子を用いて呼び出せるようにしていますので，
```Fortran
if (g .approx. f) then ...
```
などと使用できます．

32 bit実数に対する実装の全景を示すと，下記のようになります．それぞれの関数をコピペして，引数の型を`real(real64)`や`real(real128)`に変更すれば，すべての実数の型に対応できます．

~~なお，今回は実数ということで`==`演算子を上書きできませんでしたが，ユーザ定義派生型に対しては，`.approx.`ではなく`==`とするべきでしょう．~~

```Fortran
module floating_point_number_comparison
    use, intrinsic :: iso_fortran_env
    use, intrinsic :: ieee_arithmetic
    implicit none
    private
    public :: operator(.approx.)

    integer(int32), private, parameter :: WITHIN_ULP = 2
        !! 二つの実数が等しいと判定するULP（±）

    interface operator(.approx.)
        procedure :: approx_r32
    end interface

    interface get_lower_bound
        procedure :: get_lower_bound_r4
    end interface

    interface get_upper_bound
        procedure :: get_upper_bound_r4
    end interface
contains
    !>`g`が`f`の± WITHIN_ULP以内であれば真，それ以外は偽を返す．
    logical function approx_r32(g, f)
        implicit none
        real(real32), intent(in) :: f
            !! 基準値
        real(real32), intent(in) :: g
            !! 比較する近似値

        approx_r32 = (get_lower_bound(f, WITHIN_ULP) <= g .and. g <= get_upper_bound(f, WITHIN_ULP))
    end function approx_r32

    !>`f`から負の最大値方向に`ulp`ULP離れた数を返す．
    !>`f`が負の最大値の場合は負の最大値を返す．
    pure function get_lower_bound_r4(f, ulp) result(lower_bound)
        implicit none
        real(real32), intent(in) :: f
        integer(int32), intent(in) :: ulp

        real(real32) :: lower_bound
        integer(int32) :: i

        lower_bound = f
        do i = 1, ulp
            lower_bound = ieee_next_after(lower_bound, -huge(f))
            if (lower_bound == f) return
            ! lower_boundとfが等しければfは負の最大値なので，
            ! 離れた値をこれ以上探す必要はない．
        end do
    end function get_lower_bound_r4

    !>`f`から正の最大値方向に`ulp`ULP離れた数を返す．
    !>`f`が正の最大値の場合は負の最大値を返す．
    pure function get_upper_bound_r4(f, ulp) result(upper_bound)
        implicit none
        real(real32), intent(in) :: f
        integer(int32), intent(in) :: ulp

        real(real32) :: upper_bound
        integer(int32) :: i

        upper_bound = f
        do i = 1, ulp
            upper_bound = ieee_next_after(upper_bound, huge(f))
            if (upper_bound == f) return
            ! uppper_boundとfが等しければfは正の最大値なので，
            ! 離れた値をこれ以上探す必要はない．
        end do
    end function get_upper_bound_r4
end module floating_point_number_comparison
```

# まとめ
実数の比較という，プログラミング学習初期からずっとつきまとってくる問題に対して，ULPという概念を導入し，`ieee_next_after`を使った，手軽かつ確実な実装を示しました．

実数と1 ULPだけ親しくなれた気がしませんか？