幾何学的決定論において浮動小数点数が失敗する理由:整数演算の必要性

計算幾何学の世界では、わずか1ビットの違いが、操作の成功と壊滅的な失敗を分けることがあります。多くの開発者にとって、同じ入力に対して同じコードを実行すれば、実行されるマシンに関係なく同じ結果が得られるという想定があります。しかし、exact-poly の作成者が発見したように、浮動小数点数が関与する場合、これは必ずしも真ではありません。

多角形の重なりテストをデバッグしている際、著者は、頂点が凸か凹かを判定する関数が、ローカル環境では動作するものの、サーバー上では失敗することを発見しました。原因はロジックのバグではなく、異なるアーキテクチャが浮動小数点演算をどのように処理するかという点にありました。x86では、コンパイラは丸めステップを減らすために Fused Multiply-Add (FMA) を使用していましたが、WASMでは使用していませんでした。ゼロ付近に位置する頂点のこのわずかな丸めの差異が、符号の反転を引き起こし、ドミノ倒しのような影響を与えて多角形の分解を完全に変えてしまったのです。

IEEE 754 再現性の錯覚

多くの開発者は、IEEE 754 の互換性が一貫した動作を保証すると信じています。しかし実際には、この標準は主に浮動小数点数がどのように格納されるかを定義しており、異なる実行環境間でどのように動作するかを定義しているわけではありません。いくつかのメカニズムが再現性を損なう可能性があります:

  • 中間レジスタ: x87 FPU はしばしば値を 80 ビット精度で保持し、メモリに退避する際にのみ丸めを行いますが、ARM や WASM ではそうでない場合があります。
  • Fused Multiply-Add (FMA): fma(a, b, c) は、乗算と加算を単一の丸めステップで行い、個別の演算よりも正確ですが異なる結果を生成します。
  • 再結合 (Reassociation): -ffast-math のようなフラグを使用するコンパイラは、(a + b) + ca + (b + c) に書き換えることがあり、これにより丸めの順序が変わります。
  • Denormals: 「flush-to-zero」フラグはプロセスごとに密かに切り替わることがあり、非正規化数の扱いを変えてしまいます。

あるコメント投稿者が指摘したように、JVM (via strictfpStrictMath) のような一部の環境は厳密な再現性を保証しようと試みますが、浮動小数点演算の根本的な性質は近似に留まります。ロックステップ形式のゲームシミュレーションや ZK witness 生成のように、絶対的な決定論が必要なアプリケーションにとって、これらの近似は不十分です。

凸分解における「1ビット」の問題

幾何学において、cross_sign(A, B, C) 関数は意思決定のための主要なシグナルです。これは、回転が左、右、または共線(一直線上)であるかを決定します。3つの点がほぼ共線であるとき、あるマシンでは +1e-12 を返し、別のマシンでは -2e-13 を返すことがあります。

この符号の反転は、凸分解が離散的なプロセスであるため、極めて重要です。頂点は、凹であるか、そうでないかのどちらかです。もし符号が反転すれば、別の頂点が凹として識別され、異なる開始カットが行われ、全く異なる分解グラフが生成されます。エプシロン(例:「比較する前にマイクロメートル単位で丸める」)を使用しても、この問題は解決しません。なぜなら、どのようなエプシロン領域であっても、最終的には実際の多角形を捉えてしまい、動作の分岐を引き起こすからです。

解決策:正確な整数演算

x86、ARM、および WASM においてビット単位で同一の結果を得るために、exact-poly は浮動小数点数を完全に放棄し、整数の使用に切り替えています。座標を固定のスケール(例:測地学における 1 ユニット = 1 マイクロメートル)にマッピングすることで、ライブラリは外積の符号が近似ではなくビット・ロジックの問題であることを保証します。

整数予算の管理

整数演算には、スケールの慎重な選択が必要です。i64::MAX (約 $9.2 \times 10^{18}$) を使用する場合、開発者は測定単位と最大座標サイズのバランスを取る必要があります。測地学の場合、$10^6$ のスケールを使用すれば、地球の周囲の長さ全体をカバーしつつ、マイクロメートル精度の精度を維持するのに十分な余裕を持たせることができます。

i128 によるオーバーフローの防止

外積の公式 — $(bx - ax)(cy - ay) - (by - ay)(cx - ax)$ — は、容易に i64 の限界を超えてしまいます。これを防ぐために、exact-poly は厳格なルールを採用しています:座標は i64 として格納されますが、あらゆる乗算は i128 で実行されます。

重要なのは、ライブラリが減算を行う 前に 値を拡張することです:(bx as i128) - (ax as i128)。最初に i64 で減算を行うと、たとえ i128 コンテナに格納されるとしても、オーバーフローが発生してゴミデータが生じる可能性があります。

堅牢な分解カスケードの実装

単一の凸分解アルゴリズムがすべての入力に対して完璧であるわけではないため、exact-poly はアルゴリズムの「カスケード」を使用しています。最初のものが失敗した場合、次のものが引き継ぎます:

  1. ExactPartition: 凹頂点を優先する貪欲な分割器。
  2. Bayazit: Steiner points (中点頂点) を挿入できる再帰的なアプローチ。
  3. EarClip + Hertel-Mehlhorn: 三角形分割に続いて、パーツの数を減らすためのマージステップ。

3つすべてが失敗した場合、ライブラリはリングを回転させ(開始頂点を変更し)、再試行します。このヒューリスティックは、、病的なケースは特定の頂点走査順序に結びついていることが多いため、これによって解決されます。

離散数学の境界条件の扱い

整数への移行は、連続的な浮動小数点数数学とは異なる、新しい課題をもたらします:

  • 中点の精度: 単純な整数除算 (a + b) / 2 は、下位ビットを切り捨ててしまいます。ライブラリは、点のドリフトを防ぐために round_div2 を使用しています。
  • Steiner Point Snapping: 偶然の共線性を防ぐために、中点が既存の頂点に近すぎる場合は、既存の頂点にスナップさせる必要があります。
  • Bookkeeping: 最終的な出力の頂点数が入力と一致するように、合成的な頂点をタグ付けして管理する必要があります。これにより、ダウンストリームの検証失敗を防ぎます。

結論:精度よりも合意

浮動小数点数は精度(precision)のために設計されていますが、整数演算は合意(agreement)のために設計されています。2つの異なるプロセスが全く同じ答えに到達しなければならない場合、丸め方向が異なれば、浮動小数点数の精度は無意味です。中間計算に i128 を活用し、浮動小数点数を完全に回避することで、exact-poly は、64ビット整数をサポートするあらゆるアーキテクチャにおいて、幾何学的演算が同一であることを保証する決定論的な基盤を提供します。

Sources