概要
Fortranで実数（浮動小数点数）の等値比較をする場合は，ieee_next_afterもしくはieee_next_upおよびieee_next_downを用いて範囲の上下限値を設定し，その中に値が入っているかを確認すると手軽かつ確実です．
環境
- Windows 11 64 bit
- gfortran 11.2
- Intel Fortran classic 2021.4.0
注意
この記事では，浮動小数点数を実数と書きます．
実数の等値比較
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）が表示される
しかし，プログラムを学習していく過程で，実数の比較は==で行ってはいけないと学習します．==は，上記例のように二つの実数がビットレベルで等しければ真を返しますが，計算機上では実数は計算過程で丸め誤差が入るため，演算した値の等値性を==で比較すると，多くの場合等しくないという結果が得られます．
そこで，==による等値比較ではなく，二つの実数の差がある範囲に収まっていたら等しいとする技法を用いたりします．
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）が表示される
このとき，二つの実数の差をどう決めるかという問題が生じます．二つの実数の差として，マシンイプシロンが用いられることが多いように思います．マシンイプシロンは，
1.0+\varepsilon > 1.0
となる最小の$\varepsilon$のことを指します．しかし，マシンイプシロンは，定義によれば数値のオーダーが1付近でのみ用いるべきであり，fやgが大きい/小さい場合には，想定した通りの比較結果にならない場合があります．
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となります．
例で用いた数値で試してみます．
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であることが実質的にも確認できます．
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となり，マシンイプシロンよりかなり小さいことが分かります．
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()を指定することになるでしょう．
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離れた値を求めています．
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文の負荷を考慮して，軽い方を選べばよいと思います．
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が等しいとしています．
approx_r32 = (get_lower_bound(f, WITHIN_ULP) <= g .and. g <= get_upper_bound(f, WITHIN_ULP))
この関数は，interface文によって.approx.という演算子を用いて呼び出せるようにしていますので，
if (g .approx. f) then ...
などと使用できます．
32 bit実数に対する実装の全景を示すと，下記のようになります．それぞれの関数をコピペして，引数の型をreal(real64)やreal(real128)に変更すれば，すべての実数の型に対応できます．
なお，今回は実数ということで==演算子を上書きできませんでしたが，ユーザ定義派生型に対しては，.approx.ではなく==とするべきでしょう．
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だけ親しくなれた気がしませんか？