浮動小数点数の等価比較を行っても問題ありません。

浮動小数点数の等価比較は問題ありません 2026 年 4 月 14 日

注: この記事のタイトルはあえてクリックベイト（クリック誘導）です。著者もその主張には同意しますが、より正直な表現としては、「浮動小数点数を Epsilon（微小量）を用いて比較するのは OK ではありませ」といったものになるでしょう。

多くの人が耳にしている mantra（決まり文句）に「浮動小数点数の厳密な等価比較は絶対に行ってはいけない。必ず Epsilon 比較などの何らかの許容誤差を用いる必要がある」というものがあると思います。例えば、

bool approxEqual(float x, float y) {
    return abs(x - y) < 1e-4f;
}

という具合です。

私はコードを書くことについて 15 年以上経験を積み、幾何学、グラフィックス、物理シミュレーションなど浮動小数点数を日常的に扱うような業務に従事してきた者として言えることは、「Epsilon 比較」が実際に有効な解決策となるケースは極めて稀（せいぜい 1〜2 つ）であるということです。ほとんど常に、より良き解決策が存在します。それはコードを書き直すことか、単に

x == y

のように浮動小数点数を比較することかのどちらかです。そして、ほとんど常に Epsilon 比較という選択肢は、最悪の選択肢の一つになっているのです。

ここでは、「Epsilon を加えるのが第一反応になるような」多くの例を示しつつ、実際にははるかに優れており、しばしばはるかにシンプルである別の解決策が存在することを示します。まずは、浮動小数点数そのものについて話しましょう。

コンテンツ

浮動小数点数はブラックボックスではない Epsilon 比較という考え方は、「浮動小数点数が計算の神様によって不正確な結果を生じさせることがある何らかのランダムで不確かなブラックボックスだ」という一般的な認識から生じているように思えます。実際には、それは非常に確定性が高く（コンパイラオプションや CPU フラグなどを除く）、厳格に標準化されたシステムです。

浮動小数点数は本質的に不正確であり、すべての実数を表現できないからです。実際、限られたメモリ量では不可能で、数学の性質上、実数（あるいは有理数だけでも十分なほど）があまりにも多いからです。そのような数のそれぞれに対して固定した（有限ではなく、むしろ「有界な」）ビット数を割り当てる必要がある以上、表現可能な数の集合は有限であること（具体的に (2^{\text{bits}}) 個まで）を受け入れなければならず、それ以外の数については近似を扱うしかありません。

ただし、ここで浮動小数点数の表現法についての講義に転じるつもりはありません；ウィキペディアでも十分説明されていますので。より詳細な情報源としては、David Goldberg による古典的な記事や、同様の考えに基づく最新の資料があります。重要なのは、この「不正確さ」が動作の不確実性やランダム性を意味するわけではない、ということです！

例えば、二つの浮動小数点数に対する単一の算術演算（加算、乗算など）では、入力を厳密な値とみなした場合の真の答えに最も近い浮動小数点数を生成しなければならないと規定されています（同点の場合には、いくつかの丸め規則が適用されます）。結果が近似であっても、それはあくまで可能な限り真の答えに近いことが保証されており、非常に決定論的です。

しかし、これにより通常の数学的公式が常に浮動小数点数で成立するわけではないということになります。例えば、加算と乗算は結合法則（(a + b = b + a)）に従うことは保証されていますが、結合性については必ずしも成り立たないのです：((a+b)+c \neq a+(b+c)) となる場合があります。このような例を見つけるのは比較的容易であり、32 ビット浮動小数点数で動作する例は以下の通りです：

// 出力: 0.89999998
std::cout << std::setprecision(8) << ((0.2f + 0.3f) + 0.4f) << '\n';

// 出力: 0.90000004
std::cout << std::setprecision(8) << (0.2f + (0.3f + 0.4f)) << '\n';

演算結果が等しくなくても、非常に近いため（差は約 (6 \times 10^{-8})）、浮動小数点規格にはこの程度の誤差の大きさを予測するための保証があります。

では、問題は何か？ ある計算の結果は近似であることを受け入れ、ある期待値とも近似して比較するのでは、それほど悪くもありません。

それとも違うのですか？

Epsilon 比較の問題点 Epsilon 比較には私が直面している主な問題点は 3 つあります：

1. ハッカ的な一時しのぎの解決策であること
2. 極めてデバッグしにくい問題に派生しやすいこと
3. 当初の問題を本質的に解決していないことが多いこと

2 番目の点がおそらく最も難しいです。プログラムのある部分がマンハッタン距離で (1 \times 10^{-4}) の範囲内にあれば 2D の点を等しいと扱う一方で、別の部分は L-∞ 距離（座標の最大値）で (1 \times 10^{-6}) の範囲内で点を等しいと扱い、入力された点自体は別のアルゴリズムによって生成されており、その結果インバリアントが崩壊し、ラインレンダリングが壊れます。しかしそれは特定のシナリオ、特定のデータ、夜中に満月の時だけ発生する極めてニッチな不具合です。デバッグは幸運を祈るしかないでしょう。

ラインレンダリングの一時的な誤りは大きな問題にならないこともありますが、プログラムのクラッシュにいたるまで様々に現れることがあり、他の幾何学的コードと組み合わさると本当に悪あがきになります。私は多くの此类のケースに出会っており、不一致な Epsilon 比較が非常に一般的な根本原因でした。

問題は、これらの Epsilon がほぼ常に恣意的に推測される点にあります。多数の Epsilon から正しいものを選ぶ方法はありません。複数のそのような比較がある場合、Epsilon の組み合わせで正しく動作することを保証することはできません。

Epsilon 比较のもう一つの問題は、それが推移的ではないことです。これは技術的なクスクス（細工）のように聞こえますが、実際にはほとんどのアルゴリズムは「比較」のような操作が推移的であると前提しており、非推移的な比較を使うと簡単に壊れます（ナンセンスを生んだりクラッシュしたりします）。

では、どうすればよいのでしょうか？ 問題を考える必要があります！驚きですね。なぜ最初から浮動小数点数を比較しようとするのか？ 私たちのアルゴリズムを信頼していないのかもしれませんか？ データを信頼していないのかもしれませんか？ CPU を信頼していないのかもしれませんか？ 一つの正解はありませんので、いくつかの例を見てみましょう。

研究事例：グリッドベースの移動 仮に、ターン制ゲームでユニットがグリッド上を移動する状況を考えます。ユニットは一定の移動ポイントを持ち、ターンに複数回移動できますが、ユーザー体験（UX）のためには前回の移動が完了するまで次の移動を実行できないように制限します。

通常、ユニットの位置をそのままターゲットのグリッドセルに瞬移させるのではなく、何らかの方法で補間しているため、プレイヤーが次の移動ターゲットを選択することを許可する前に、移動が完了した正確なタイミングを確認する必要があります。

単に一定時間の待ち時間を課すことも（多くのケースでは完全に有効な解決策です！）、しかし異なるユニットには異なるアニメーションとそのタイミングがあり、アクセシビリティ設定でアニメーションを減らす機能などもあるため、アニメーション時間への依存は良いアイデアではないと判断されます。

そこで、移動はユニットの位置がターゲットセルの中心と一致する正確に完了した瞬間であると気づきます。以下のようなコードを書きます：

void update() {
    if (selectedUnit.position != targetCell.center)
        return;

    // フレーム更新処理
}

移動が実行された後でも、

selectedUnit.position != targetCell.center

が決して真になることがなく、ゲームが実際にはフリーズします。緩和関数付きの典型的なリニア補間を用いたコード：

vec2f getPosition(float time) {
    return lerp(start, end, easing(time));
}

の場合、時間

time == 1.f

のときに結果が

end

に等しくなることはありません。浮動小数点数の丸め誤差により十分大きいためです。実は、一部の補間手法ではその条件が決して成り立たないことも珍しくありません！

「Damn, ストイプな浮動小数点数！」と溜息を吐きつつ、あなたは浮動小数点数の聖杯である Epsilon 比較を追加します：

void update() {
    if (distance(selectedUnit.position, targetCell.center) > 1e-4)
        return;

    // フレーム更新処理
}

これで問題は解決するようですが、なぜこれが悪手なのか？理由は様々です：

• 特定の Epsilon（1e-4）は、補間手法を変更した場合に崩れてしまう可能性があります。
• いつものように Epsilon を使うと、コードを読む人が混乱し疑念を抱いてしまいます。
• プレイヤーを待たせることはすべき悪手の一つです。

ではどうすれば解決できるのでしょうか？一つの方法は「受容半径」を使用することです。ユニットがターゲットセル中心から例えば 0.25 の範囲に入った時点でアニメーションを停止させ、プレイヤーに次のコマンドを入力させることができます。その後、大量のテストによって最適な値を見出します。

これは Epsilon と何が違うのか？違いはむしろある程度バックアップされた値になっている点で、ランダムな値ではないことです。ここでの本質的な問題は、データモデルとその表現方法を混同している点です。ゲームの状態機械の内部動作が 3D モデルの位置を気にする必要はありません。それは一見難しそうに思えます（それが UX が成立する理由でもあります）ですが、十分に可能です。

この場合、私の意見では最良の解決策は、ユニットの移動が完了するのを待たずにユーザーがコマンドを入力できるようにすることです。ユーザーがグリッド上のセルをクリックすると、レンダリング側に移動が発生したことが通知され、一方ゲームの状態モデル内部ではユニットがすでに次のセルにいると考えるようになります。

これを 구현する方法はいくつかありえます：要求されたアニメーションをキューに入れ順次再生するとか、ターゲット値の急激な変化に頑健な補間/アニメーションスキームを使用するなどです。重要なのは、これが実現可能で比較的シンプルであり、Epsilon を必要としないことです。

研究事例：スフェリカル・ラインアー インターポレーション (Slerp) スフェリカル・ラインアー インターポレーション（球面上の点、すなわち単位ベクトル）を補間する手法で、補間されたベクトルは固定された回転速度を持つ軌道に従います。単純な

normalize(lerp(a, b, t))

では不十分であり、結果としてのベクトルは両端では遅く、中央では速くなります。Slerp の視覚化がこちらです。この関数はしばしばクォータニオン（4 次元複素数）専用実装としてしか見ませんが、任意の次元のベクトルにも有用です（クォータニオンの場合は (q) と (-q) が同一の回転を表すため少し異なります）。

入力が正規化されたベクトルであれば、以下のシンプルで直感的な式に帰結します：

vec3 slerp(vec3 a, vec3 b, float t) {
    float angle = acos(dot(a, b));
    return (sin((1 - t) * angle) * a + sin(t * angle) * b) / sin(angle);
}

このコードは時々 NaN を生成します。理由は次の通り：

• ベクトルが正規化されていても、ドット積の結果は ([-1, 1]) の範囲を外れることがあり、その場合

  acos

  は NaN を返します。
• 入力ベクトルが非常に近い場合、

  angle

  がゼロになり、ゼロ除以ゼロが発生して再び NaN が生成されます。

最初の問題は扱いやすい：

acos

の引数を

clamp

で囲むだけです。

vec3 slerp(vec3 a, vec3 b, float t) {
    float angle = acos(clamp(dot(a, b), -1.f, 1.f));
    return (sin((1 - t) * angle) * a + sin(t * angle) * b) / sin(angle);
}

二つ目の問題はより微妙です。幸运的是、角度が十分に小さい場合、スフェリカル・ラインアー インターポレーションは通常の直線補間と同じ振る舞いになるため、このケースを検出して通常の

lerp(a, b, t)

に切り替えることができます。ただし、「どの程度小さくてよいのか？」という問題があります。

ここで適当な Epsilon を使うのが誘惑的です。

vec3 slerp(vec3 a, vec3 b, float t) {
    float angle = acos(dot(a, b));
    if (angle < 1e-4) {
        return lerp(a, b, t);
    }
    return (sin((1 - t) * angle) * a + sin(t * angle) * b) / sin(angle);
}

しかし、これによりコードは本来の精度よりも低下してしまいます（結果としてのベクトルが期待通り近いたとしても）。そしてまた、ランダムな Epsilon を使うと「なぜその正確な値を選んだのか？」という問いが生じます。もっとよくできるはずです！

GLM や Eigen では妥当なチェックを行います：

vec3 slerp(vec3 a, vec3 b, float t) {
    float d = dot(a, b);
    if (d > 1.f - FLT_EPSILON) {
        return lerp(a, b, t);
    }
    float angle = acos(dot(a, b));
    return (sin((1 - t) * angle) * a + sin(t * angle) * b) / sin(angle);
}

ここで

FLT_EPSILON

は正確に (2^{-23}) で、浮動小数点加算において (1 + x \neq x) となる最小の単精度浮動小数点数です。正直なところ、乱数 Epsilon と比べてさほど良くなっているとは思えません——

FLT_EPSILON

を使うだけでは、これが実際の場合に有効な閾値であることを証明する限り、私たちが直面する問題は解決されません。

コードを分析してみましょう。主な問題は角度がゼロになることです。二次的な問題は角度が小さすぎると精度誤差が生じることです。

まず精度について考えます。

angle

は単なる任意の数字ではありません——それは

acos

を呼んだ結果であり、特に

acos

の引数が 1 に近い場合に心配しています。実際には角度そのものよりも、

\sin(angle)

に対してはあまり気にしません。

さて、(\sin(\operatorname{acos}(x))) は単に (\sqrt{1-x^2}) です。ここでは (x = \text{dot}(a, b)) です。(x = 1 - \varepsilon) の場合、

[ \sin(\operatorname{acos}(x)) = \sqrt{1-x^2} = \sqrt{1-(1-\varepsilon)^2} = \sqrt{2\varepsilon-\varepsilon^2} \approx \sqrt{2\varepsilon} ]

となります。(\varepsilon^2) 項は無視でき、(\sqrt{2}) は約 1 に近い定数です。ここでは主に (\sqrt{\varepsilon}) が重要です：たとえ (\varepsilon) が

FLT_EPSILON

程度に小さくても、

\sin(angle)

の項はそれの平方根程度の値になります。小さな数に対して平方根を取るとその値が大幅に増大します（例えば、

FLT_EPSILON

は約 (1.2 \times 10^{-7}) ですが、その平方根は約 (3.5 \times 10^{-4})、つまり約 2000 倍大きくなります）。

私が言いたのは、ここでは

acos

の引数が 1 に近いケースにおいて、角度自体がゼロに近いことよりも、

acos

の引数の表現可能な値同士の差分の方が精度問題の重大な原因となるということです。

実際には、浮動小数点数で表現可能で 1 よりも小さな最小値は

1 - FLT_EPSILON

ではなく、

1 - FLT_EPSILON / 2.f

です。これは私の主張（

FLT_EPSILON

定数がここである程度恣意的に選ばれた）を部分的に支持します。

特別なことを望まない限り、精度の観点から角度が小さいこと自体は無視でき、NaN を回避することに焦点を当てることができます。このコードにおいて引数がある値である場合、NaN が発生するのはゼロ除以ゼロによるもののみで、これは角度がゼロの場合のみ起こり得ます。これは

dot(a, b) >= 1

の場合にのみ発生します。前述した通り、

dot(a, b)

の次の小さい表現可能な値は

1 - FLT_EPSILON / 2.f

で、対応する角度は約 (\sqrt{\text{FLT_EPSILON}})（約 (3.5 \times 10^{-4})）であり、これは全範囲 [(10^{-38} .. 10^{38})] に比較してかなり大きな浮動小数点数です。したがって、中途半端な実装であっても

acos

はゼロを返さないことが確実です。

つまり、

dot(a, b) < 1

の場合、実際には問題ありません！確認すべきは

dot(a, b) == 1

であり、先に追加した

clamp

と組み合わせて言えば

dot(a, b) >= 1

です：

vec3 slerp(vec3 a, vec3 b, float t) {
    float d = dot(a, b);
    if (d >= 1.f) {
        return lerp(a, b, t);
    }
    float angle = acos(dot(a, b));
    return (sin((1 - t) * angle) * a + sin(t * angle) * b) / sin(angle);
}

（注：ドット積が (-1) に近いケースは意図的に除外しています。その場合、補間問題は一意の解を持たず、別途扱うべきです。）

あるいは、

sin(angle)

を直接チェックすることもできます。ただし、これは角隅ケースで高価な

acos

呼び出しを節約できなくなります：

vec3 slerp(vec3 a, vec3 b, float t) {
    float angle = acos(clamp(dot(a, b), -1.f, 1.f));
    float s = sin(angle);
    if (s == 0.f) {
        return lerp(a, b, t);
    }
    return (sin((1 - t) * angle) * a + sin(t * angle) * b) / s;
}

研究事例：ベクトル長を計算する 私は「ベクトルの（ユークリッド）長さを計算することは非常に一般的な有用な演算だ」と言ったとしても、あなたが驚くことはありません。実装も比較的単純です：

float length(vec3 v) {
    return sqrt(dot(v, v));
}

またはドット積をインライン展開した場合：

float length(vec3 v) {
    return sqrt(v.x * v.x + v.y * v.y + v.z * v.z);
}

これを目にしたことは十回以上あるでしょう。何が悪いのでしょうか？特に何もありません。ただし、非常に小さなベクトル（長さが約 (10^{-19}) のもの）の場合にはゼロを返し、大きなベクトルの場合には有意な精度損失が生じます。結果が有効な浮動小数点数値であるにもかかわらず、その平方根の内部にある式（二乗和）が浮動小数点で適切に表現できるほど小さくなりすぎる可能性があります。

ではどうしましょう？そのようなベクトルをゼロと区別できなくても実際の使い道はほぼありません。しかし、「長さゼロなのは零ベクトルだけ」という不変式を持つことは極めて便利であり、特に浮動小数点を慎重に扱う場合に多くのコードで記述を簡素化できます。

実はこれを実現するための非常にシンプルな方法があります：ベクトル座標の絶対値のうち最大値 (M) を計算し、

M * length(v / M)

を返しますというものです。浮動小数点保証により、少なくとも一つの座標は 1 より大きいので、平方根内の式は最小で 1、最大で (D)（vec3 の場合は次元数 3）です。

したがって、コードは以下のようになります：

float length(vec3 v) {
    float M = max(abs(v.x), abs(v.y), abs(v.z));
    vec3 u = v / M;
    return M * sqrt(u.x * u.x + u.y * u.y + u.z * u.z);
}

アセンブリ命令数は約 2 倍増えますが、これはホットループ内ではない限り、追加の精度と安全性に値するでしょう。

唯一の問題は、ベクトルが文字通りゼロ（単なる小ささではなく）の場合にはこの関数が NaN を返してしまうことです。しかし、少なくとも一つの座標がゼロでない場合（サブノーマル値であっても）、(M) は厳密に正になり、アルゴリズムは正しく動作します（ naïve なバージョンよりも低い精度で決して動作しません）。したがって、正しい直し方は文字通り (M) をゼロと比較することです：

float length(vec3 v) {
    float M = max(abs(v.x), abs(v.y), abs(v.z));
    if (M == 0.f) return 0.f;
    vec3 u = v / M;
    return M * sqrt(u.x * u.x + u.y * u.y + u.z * u.z);
}

もし

if (M < 1e-4)

と書いたら、この関数が存在する意義をほぼ完全に破壊してしまいます。

研究事例：線形方程式系を解く 「線形方程式系を解くことは非常に一般的で有用な演算だ」と私が言ったとしても、あなたが驚くことはありません。物理学・工学問題の半分近くは、ある種の線形方程式系の求解に還元されます（もう半分は固有値方程式）。

一般的な非疎行列系に対して、基本的に標準的な方法は Gauss-Jordan 消去法のみです（実際には QR 分解も可使用できます——知見によるとより遅いですが数値的に安定します）。一般アルゴリズムは少し長くなりますが、それほど複雑ではなく、単に for ループの数と係数に対する算術処理だけです。

重要なのは、アルゴリズムはある点で非常に不安定になる可能性があることです：ある時点で残りの行列部分の左上要素を取り出して、その要素で一行全体を割ります。これが小さい場合、打ち消し合いらしい誤差が拡大され、ゼロの場合には NaN が発生します。

ただし、もし「partial pivoting（部分選定）」と呼ばれることを代わりに行う場合、実用上は安定になります（不安定性は極めて稀であり理論的なものです）。現在の列で絶対値が最大のコエフィシエントを持つ行を見つけ、それを元の上部の行の代わりに使用します。アルゴリズムを少し複雑にはなりますが、実用上非常に使いやすくなります。

問題なのは、今でも割り算を行うことであり、非常に小さい数やゼロによって割られる可能性があることです。この問題を巧妙に回避する手法はありません：行列が特異（singular）な場合、数学的に解が存在しないことが分かっています。

したがって、まず私たちの関数は

optional<vector>

を返し、そしてその割り算をチェックする必要があります：

optional<vector> solve(matrix const& m, vector const& v) {
    // Gauss 消去法のコード...

    for (int i = 0; i < m.columns(); ++i) {
        
        // 残りの行の中から最大の m[j][i] を見つける
        float M = 0.f;
        for (int j = i; j < m.rows(); ++j)
            M = max(M, abs(m[j][i]));

        if (M < 1e-4)
            return std::nullopt;

        // アルゴリズムを継続
    }
}

さて、当然ながらこの

1e-4

は何処から来るのか——単に問題解決せず黙らせるための典型的な Epsilon です。特異に近いた行列系に対しては、私たちが非常に恣意的な閾値に基づいて失敗と報告するでしょう。より悪辣なのは、この計算した (M) 値は行列の本質的な性質ではなく、アルゴリズムで得た中間値に過ぎないことです。

この閾値は絶対的にユーザー自身が提供すべきものです（デフォルトパラメータもあるかもしれません）が、それでも疑問があります。行列が特異であるか、それに近いのかを適切にチェックする方法は条件数計算や特異値の検査であり、恣意的な中間値に対して恣意的な閾値を設定することではありません！

私の主張は：行列の特異度がどの程度かを知ることは私たちの仕事ではなく、ユーザーの仕事です。Gauss-Jordan 消去法を実装する者としての私たちの仕事は、可能な限り最良の答えを提供するか、そうでない場合は失敗と報告することです。我々のコードが本当に失敗できるのはゼロ除以ゼロによる場合のみであり、それ以外のあらゆるものは合理的な答えを生じさせます（行列によっては非常に不正確でも）。

したがって、私の意見では正しいコードは次のようになります：

// ...

if (M == 0.f)
    return std::nullopt;

// ...

正直に言うと、これですら大きな数で小さな数を割って無限大を得ることを防いでいない——しかし Epsilon はいずれの場合も保護できません。

研究事例：レイとボックスの交差 (Ray-Box Intersection) レイトレーサー（ボクセルやそうでなくても）を作成するか、マウスによるオブジェクトピッキングを実装するか、あるいは何百万もの他のタスクを行う場合、ray-box 交差ルーチンが必要になります。

アルゴリズム自体は比較的単純です：レイがボックスに入る時間（パラメータ (t)）と出る時間を計算します。前者が後者より小さい場合が交差点です。入り時間は各座標軸における入り時間の最大値として計算され、同様に通り時間は最小値として計算されます。

この Scratchapixel の記事はこれをよく説明しています。コードはおよそ以下のようになります：

void sort(float & x, float & y) {
    if (x > y) swap(x, y);
}

pair<float, float> intersect(ray r, box b) {
    float txmin = (b.min.x - r.origin.x) / r.direction.x;
    float txmax = (b.max.x - r.origin.x) / r.direction.x;
    float tymin = (b.min.y - r.origin.y) / r.direction.y;
    float tymax = (b.max.y - r.origin.y) / r.direction.y;
    float tzmin = (b.min.z - r.origin.z) / r.direction.z;
    float tzmax = (b.max.z - r.origin.z) / r.direction.z;

    sort(txmin, txmax);
    sort(tymin, tymax);
    sort(tzmin, tzmax);

    float tmax = min(txmax, min(tymax, tzmax));
    float tmin = max(txmin, min(tymin, tzmin));

    return {tmin, tmax};
}

注:

sort

を少し書き換えることで、minss/maxss SSE 命令を使用を強制でき、したがって

intersect

関数をブランチレスにできます。

非常に短く可読な関数で、99.99% のケースでは完璧に動作します。それでも時折、

direction.x

がゼロになるレイに遭遇し、除算が無限大または NaN を与えてしまいます。Oops!

そこで、いくつか Epsilon を追加しましょう？いいえ。まず、

abs(direction.x) < 1e-4

の場合にどうするか考えましょう？それでも

direction.x

で割らずに適切に交差を見つけられるコードを追加する必要がありますし、Epsilon なしで機能するかもしれません。

コードを分析してみましょう！

direction.x == 0

の場合何が起きますか？レイは YZ 平面に平行なので、真の交差点はその平面上にあります。しかし X 座標を無視することはできません：もしレイの始点

[b.min.x, b.max.x]

区間の外側にあるなら交差はなく、そうでなければ Y と Z 座標を検査する必要があります。

実際に

direction.x == 0

の場合何が起こるか？A を

direction.x

で割った結果は、A < 0 なら (-\infty)、A > 0 なら (\infty) を返します。私たちのケースでは A は

b.min.x - r.origin.x

または

b.max.x - r.origin.x

です。もしレイの始点

r.origin.x

が区間

[b.min.x, b.max.x]

の内部にあるなら、

txmin == -INF

で

txmax == INF

になります。始点がこの区間外にある場合、両方の

txmin

と

txmax

はどちらも (-\infty) または

\(+\infty

) です。

交差がありかつ始点が区内にある場合、

tmin

の計算は実質的に

txmin == -INF

の値を無視します。なぜなら常に

max(-INF, x) == x

であるからです。同様に、

tmax

も実質的に

txmax == INF

を無視します。この場合、コードは YZ 平面で正しい交差点を計算します。

しかし、始点が区内にない場合、

txmin == txmax == INF

で

tmin == INF

となりつつ

tmax

が有限値になる（→ 交差なし）、または

txmin == txmax == -INF

で

tmax == -INF

となりつつ

tmin

が有限値になる（→ 再び交差なし）となります。

つまり、コードは実際には

direction.x == 0

の場合正しく動作しています！Epsilon は不要です。無限大を適切に処理しています。それともそうでしょうか？

忘れ去ったことがあります：IEEE754 ではゼロには正の +0 と負の -0 の二種類があります。負のゼロで割ると無限大になりますが、符号も反転します！したがって、

txmin == INF

で

txmax == -INF

となり交差がないと報告されるケースも発生します。

幸いなことに、私たちの

sort(txmin, txmax)

呼び出しがこの問題を解決し、これを

txmin == -INF

で

txmax == INF

のケースに置換します！再び、コードはこの場合でも正しく動作しています。

上記の Scratchapixel リンクでは、この問題に対する解法のうち、多くの if ステートメントを含むものが紹介されていますが、これはブランチレスにはならず、早期リターンも可能となっています。

実際にはまだ修正していないケースがあります：

r.direction.x == 0

および

r.origin.x == b.min.x

の両方が真の場合、ゼロ除以ゼロになり、NaN が発生します。残念ながら、このコードは自動的にこれを解決しません：NaN を何と比較しても常に偽を返すため、

sort

,

min

,

max

の実装方法によっては、

tmin

または

tmax

が NaN となり、結果として交差がないと扱われます。これはあなたのユースケースによっては問題ないかもしれませんが、これを真の交差として報告したい場合、単純な修正方法が分かりません。

r.direction.x == 0 && r.origin.x == b.min.x

をチェックし、この場合は

tmin = -INF

に設定すればよいですが、コードはブランチレスではなくなります。ただし実用上これは極めて稀に発生するため、これで済むかもしれません :)

研究事例：凸包を計算する ほとんどの 2D 凸包アルゴリズムは最終的には、三つの点 (A, B, C) について最後の点 (C) が線分 (AB) の右側にあるかをチェックすることになります。等価的には、(A) から (B) へ移動しその後 (B) から (C) へ向かうとき、左回転するか右回転するかの判断を行います。これがこの述語がしばしば

left_turn

と呼ばれる理由です。

これは二本のベクトル (AB) と (AC) の相対方向をチェックすることに帰結し、単純な行列式で計算できます：

float det(vec2 a, vec2 b) {
    return a.x * b.y - a.y * b.x;
}

bool left_turn(vec2 a, vec2 b, vec2 c) {
    return det(b - a, c - a) > 0.f;
}

すべての浮動小数点処理をこの述語にまとめることができて非常にcool です——凸包アルゴリズム自体は座標に関心せず、この

left_turn

述語を何らかのブラックボックスとして扱える抽象アルゴリズムとして動作できます。

ただし、計算幾何学のほとんどすべてのアルゴリズムはエッジケースに対して極めて敏感であり、このようなエッジケースがアルゴリズムの不変性を崩すと完全に壊れてしまいます。

left_turn

述語の一つの不変性は、点を巡回（cyclic）したときに変化しないことです：

left_turn(A,B,C) == left_turn(B,C,A) == left_turn(C,A,B)

三つのほぼ一直線上にある点を見つけやすいです。これが非常に頻繁に発生するため、少なくとも 1000 点程度を処理する必要がある凸包アルゴリズムはこれを考慮する必要があります。

解決策は知っていますよね？— ただここに Epsilon を付け加えるだけです：

bool left_turn(vec2 a, vec2 b, vec2 c) {
    return det(b - a, c - a) > 1e-4f;
}

もちろんこれは機能しません。浮動小数点誤差は簡単にこの述語が上記の巡回不変性を破るのを招きます。

この問題を解決する多くの方法があります：

• 入力値を固定グリッドに丸めて、整数で作業する
• 任意精度の有理数算術を使って真の結果を得る
• CPU の丸めフラグを利用した区間算術を使って結果を計算し、失敗したら任意精度に頼る
• 正確な IEEE754 丸則を知り、誤差の上界を計算してナィーブな結果と比較し、誤差が結果より小さい場合は戻す；さもなければより高精度の方法を使う
• 2Sum アルゴリズムの某种variationを使って結果を計算し、失敗したらより高精度の方法を使う

これらのうち、最初の選択肢（固定グリッドに丸める）はおそらく最もシンプルで実用的です。ただし、それができない場合や未乱された入力データを扱う必要がある場合、より複雑な方法が必要になります。いずれにせよ、Epsilon だけはあなたを救うことはできません—— sooner or later あなたのアルゴリズムはナンセンスを返したり、クラッシュしたりします。

いずれのケースでも、そのような述語は三つの可能な値（左、右、共線）を返すものとして考えられる方が良いでしょう（-1, 0, 1 のように'90s の感覚なら）。これは三路比較と非常に似ています。

研究事例：ユーザ入力をサニタイズする 最後の二つのケースでは、私は Epsilon が実際に良い解決策だと考えます！

例えば、幾何可視化ライブラリや、何か種の地測量エンジンを書いているとします。ユーザはあなたが制御できないポリラインを供給します。あなたの任務はそれに計算を適用し、レンダリングすることです。

ポリラインのレンダリングは比較的複雑なトピックです；特に、我々は通常、いくつかの種類のコネクション（接合部）やキャップを使用してポリラインを三角分割します。これらを計算するには、特定の線分への法線や二つの連続する線分間の角度など様々なものが必要です。そして当然のことですが、二つの連続する点の位置が等しいか、単にあまりにも近ければ全て壊れます——法線と角度の計算が大きく逸れ、我々の愛すべきラインに醜いアーティファクトを生じます。

ここでは「連続して等しい点は許されない」と言うことも（あるライブラリやアルゴリズムには合理的な要件ですが）、ユーザ向けのエンジン/フレームワークに対してはそれほどではないかもしれません。畢竟、等しい点をフィルタリングするのは容易です——文字通り

std::unique

の単一呼び出しで済みます。しかし、あまりにも近いが完全には等しくない点についてはどうすればよいでしょうか？

数学的には、点が互いに非常に近くても問題ではありません—これらは出力に非常に薄い三角形を生み出すだけで、メッシュのトポロジーを崩すことはありません。ただし、実際には浮動小数点の精度不足によりこれらの薄い三角形を正確に計算できず（これが醜いアーティファクトの原因です）、一方では、今回は最も高い精度で何かを計算するのではなくデータを可視化することが目標です。実際には誰もそのような薄い三角形を見ないでしょう——単にそれらがあまりにも薄すぎるためです。

つまり、レンダリングしないこと自体が問題ないということです。生成さえしなくても問題ないということです。入力の点をフィルタリングしても問題ないということです。しかし「どれくらい近いと問題なのか」？

再び乱数の Epsilon をつけるのは通常悪い解決策ですが、この特定のケースではおそらく悪くなく、ただまあまあです。この特定のユースケースに適した良い Epsilon を見出す方法があります：

• 最大スケーリング/ズームと画素サイズを知り、半画素未満にマップされる距離を計算し、それをフィルタリングの Epsilon として使用
• アルゴリズムが何を doing か、不正確さがどこから来るかを理解し、受け入れ可能な浮動小数点誤差をもたらす点間の最小分離距離を見積もり、それを Epsilon として使う
• バイナリ検索を使用して、入力データセットに対してテストすることで良い Epsilon を見出す
• 数十年の浮動小数点経験を活かし、良い Epsilon を推測する（ただし文書化を忘れずに）

最後の選択肢を選ばない限り、私は正しい値が

1e-4

または

1e-6

ではないことを保証します。

研究事例：テストケースを書くこと これはおそらく最も明白なケースです。もし私が数学ライブラリを書いているなら、広範なテストカバレッジを望みます。テストでは関数が期待する値を返すかどうかをチェックする必要があり、つまり結果を何らかのリファレンス値と比較する必要があります。しかし浮動小数点誤差のため、結果はほぼ常に異なります。

このことは二つの方法で行えます。オプション 1 は厳密な IEEE754 の準拠を仮定し、等価比較を用いたテストを書くことです。たとえ丸めが関与しても、それはすべての準拠ハードウェアで同じである必要があります（CPU フラグが適切に設定されている場合）ので、安全に等価テストを実行できます。多くのケースでは、丸めが全く関与しないように入力データを工夫することも可能です！

多くの関数（例えばベクトルの加算、減算、乗算など）では、整数座標（多くは浮動小数点数として正確に表現できます）や 2 のべき乗で割った整数を使用できます。いくつかのケースではこれが十分ではありません：例えばベクトル長を計算する場合でも、座標が整数であっても結果が無理数になる可能性があります。

この場合でも少しチートしてピタゴラスの三つ組を利用して、浮動小数点数で正確に表現できる長さを持つテストケースを作成できます：

length(vec2(3.f, 4.f)) == 5.f
length(vec2(5.f, 12.f)) == 13.f
length(vec2(0.4375f, 1.5f)) == 1.5625f

これは非常に手間がかかります。すべてのコードに対して特別に複雑になるテストを工夫する必要があります。オプション 2 は何でしょうか？

Epsilon を使う！比較的小さい Epsilon（例えば

FLT_EPSILON

またはそのようなもの）を使うことだが、実はこのケースでは良い解決策です。大部分において、数学ライブラリを破ることは非常に難しく、Epsilon 比較がそれを見つけられない可能性は低いです。

ただし、いくつかの関数が追加の不変保証を提供する場合（例えば長さゼロはゼロベクトルに限られる場合や、

left_turn

が引数の巡回下で不変の場合など）には、これらの保証のために専用の特別なテストを書くのが最良です。

Epsilon を使うか使わないか

結局、Epsilon は良いのか悪いのか？通常は悪いが、時には OK である。今この記事を読んだことで、あなたの全 Epsilon に依存するコードが突然壊れると思うでしょうか？おそらくいいえでしょう。

実際には多くの現実のユースケースでは、これらに関心を持つ必要はありません。もしあなたが汎用高品質な数学ライブラリを書いているなら、このことに気を配る義務があります。もしあなたがクイッシュプラットフォームゲームのための玩具物理エンジンを作っているなら、Epsilon はあなたに適しているでしょう。そして壊れたとしても、

1e-4

を

2e-4

に置き換えておけばよいのです。

エンジニアリング問題に対する真の答えは、カルトやランネット上の人物に従うことではなく、自分の頭で考えることです。驚きですね、私も。しかし、何か学到了かと思いました！