浮動小数点数の10進数演算：楽しみと利益のために

多くの人が、米ドルやセントのような10進数の計算を、ほとんどのプログラミング言語の浮動小数点数で行うべきではないと知っています。これは、10進数をその浮動小数点数として正確に表現できないため、丸め誤差に遭遇するからです。有名な例として、IEEE倍精度浮動小数点数では、0.1 + 0.2 は 0.3 ではなく、0.30000000000000004 になります。一方で、私はPythonのREPLを使って、常にレシートの10進数を合計しています。もちろん、私のリスクは低いです。私は1ドルから100ドルの範囲のものを合計しており、微細な誤差を手動で丸めてから合計をどこかにコピーすることを知っています。しかし、実際には、それらの誤差が全く現れないこともかなりよくあります。0.01の倍数を1.00まで全てのペアで合計し、結果が大きすぎる（青）、小さすぎる（赤）、または正しいかどうかを見ると、クールなパターンが得られます：図1：浮動小数点数が0.01の2つの倍数の合計にどのように影響するか、1.00まで。このパターンはどこから来るのでしょうか？浮動小数点数、簡潔に言うと倍精度浮動小数点数 x は、1つの符号ビット、11個のエキスポーネントビット、および52個のフラクションビット（この順序で）で構成されており、合計64ビットです。符号ビットは、+ または - の符号を提供します。エキスポーネントビットは、-1022から1023までの整数 E を表します。フラクションビットは、2^52 未満の非負整数 F を表し、これはシグニフィカンド 1 + F/2^52 ∈ [1, 2) に対応します。この常に存在する 1 は「隠れビット」と呼ばれます。浮動小数点数の値は x = ± 2^E(1 + F/2^52) です。したがって、最後のフラクションビットの実効値は 2^(E-52) となり、これは x の ulp、「unit in the last place」と呼ばれる量です。x に最も近い2つの浮動小数点数は、1つの ulp だけ x から異なります。ただし、F = 0 で x がちょうど2のべき乗であるという例外的なケースを除きます。例：「隠れビット」52ビット float(0.1) = 0.00011001100110011001100110011001100110011001100110011010₂ ulp(float(0.1)) = 0.00000000000000000000000000000000000000000000000000000001₂ この説明は目的に十分ですが、0、サブノーマル数、無限大、NaNなど、他の多くのケースを無視しています。これらは私が説明しなかったエキスポーネントビットの2つの値を使用します。簡単のため、他の精度の浮動小数点数も考慮しません。加算の解剖学例として（qntmなどに倣って）、Python REPLで 0.1 + 0.2 と入力したときに何が起こるかをステップごとに見てみましょう。まず、Pythonの式 0.1 はある浮動小数点数に評価されます。具体的には、IEEE標準で要求されるように、0.1に最も近い浮動小数点数です。その正確な数値を float(0.1) と書くことにします。同様に、Pythonの式 0.2 は別の浮動小数点数、float(0.2) に評価されます。次に、+ が2つのステップで評価されると想像できます。まず、Pythonは正確な合計 float(0.1) + float(0.2) を計算します。第二に、それを最も近い浮動小数点数に丸めます。結果の値は float(float(0.1) + float(0.2)) です。（これは文字通りには動作しません。最初のステップの正確な合計はどこにも実現されませんが、数学的には正確です。）実際には、私がこれを非常に詳細に調べるまで気づかなかった微妙な点があります：float(0.1) + float(0.2) は、その2つの最も近い浮動小数点数に等しく近いのです！この場合、Pythonは偶数のシグニフィカンドを持つ浮動小数点数に丸めます。この場合、切り上げられます。最後に、この値を表示するために、Pythonはそれを10進数に変換する必要があります。この変換は驚くほど微妙であり、IEEE標準で正確に指定されているわけではありません！ float(float(0.1) + float(0.2)) は正確に 0.3 ではありませんが、非常に近いので、「0.3」と表示しても不合理ではありません。別の選択肢は、それを正確に「0.3000000000000000444089209850062616169452667236328125」と表示することです。さまざまな丸め量で表示することも想像できます：0.300000000000000044、または0.30000000000000004441、などです。float(float(0.1) + float(0.2)) = float(0.30000000000000005) であるという事実を考慮して、例えば 0.30000000000000005 を考えることさえできます。つまり、Pythonの式 0.1 + 0.2 == 0.30000000000000005 は真です。この変換の微妙さは、Python 3.1 の特定のパッチまで、REPLに 1.1 と入力すると、Pythonは入力を 1.1000000000000001 と表示していたという事実によって証明されています。この10進数変換がどのように機能すべきかについての標準的な説明は、1990年のSteeleとWhiteによって形式化されました。彼らは3つの基準を定めています：1. 10進数表現はラウンドトリップ可能であること：再度入力すると、同じ浮動小数点数が得られるはずです。これは「0.3」という出力を除外します。2. 基準1を条件として、10進数表現は可能な限り短くすること。これは「0.30000000000000004441」のような出力を除外します。3. 基準1と2を条件として、10進数表現は浮動小数点数に可能な限り近くなること。これは「0.30000000000000005」のような出力を除外します。このアルゴリズムに従うことで、0.1 + 0.2 == 0.30000000000000004 となる理由を理解できます。一般化しましょう一般的なケースを扱います。正確な正の10進数量 a と b を加算しようとしていると仮定します。その正確な合計は c です。まあ、完全には一般的ではありません。これらの値は「妥当な金額」であり、正で70兆ドル未満であると仮定します。これはレシートのファイルに十分だと思います。それより少し上（246 = 70,368,744,177,664）になると、浮動小数点数はセントの倍数よりも疎になり、これは良くありません。問題は、float(float(a) + float(b)) が float(c) とどのように比較されるかです。それらの差を Δ := float(float(a) + float(b)) - float(c) とします。誤差関数 error(x) := float(x) - x を導入することで、これを推論できます。すると、Δ は次のように書き換えられます： Δ = error(a) + error(b) + error(float(a) + float(b)) - error(c)。(*) 各誤差項をバウンドできます。浮動小数点数 x の ulp（unit in the last place）は、最後のビットの値であることを思い出してください。記法を少し乱用して、x が浮動小数点数として正確に表現できない場合でも ulp(x) と書くことを許可し、これは ulp(float(x)) を意味すると理解します。したがって、float(x) ± ulp(x) も浮動小数点数であり、x から float(x) 自体よりも正確に近くないはずです（そうでなければ、float(x) はより近い値に評価されていたでしょう）。これは、すべての（妥当な）x について、|error(x)| ≤ ulp(float(x))/2 が成り立つことを意味します。さらに、x が2つの最も近い浮動小数点数のちょうど中間に位置する場合にのみ等号が成り立ちますが、x が妥当な金額である場合にはこれは成り立ちません。このバウンドを (*) の各項に適用すると、|Δ| < 2ulp(c) と結論付けられます。さらに、Δ は c の近くの2つの浮動小数点数の差であるため、ulp(c) の倍数です。これにより、Δ ∈ {-ulp(c), 0, +ulp(c)}、つまり結果は答えから最大1 ulpずれる可能性があると結論付けられます。しかし、ここで Δ の中間バウンドをよりタイトにする導出を示します。簡単のため、a ≤ b と仮定します。すると、float(c) + ulp(c) - float(b) は、結果の ulp が b と c の両方の ulp 以下であるため、表現可能な浮動小数点数です。したがって、少なくともこれは a の利用可能な近似です。そしてそれは過大評価です： a = c - b ≤ float(c) + ulp(c)/2 - float(b) + ulp(c)/2 = float(c) - float(b) + ulp(c)。したがって、