IEEE754浮動小数点数と数値誤差 — 0.1+0.2が0.3にならない理由

Pythonを起動して 0.1 + 0.2 と打ち込むと、返ってくるのは 0.30000000000000004 です。0.1 + 0.2 == 0.3 は False になります。小学校で習った足し算がコンピュータの中では成り立たない——初めて見たときは、多くの人がバグを疑います。

しかしこれはバグではありません。CPythonの実装ミスでもなければ、あなたのCPUの欠陥でもありません。世界中のほぼすべてのプロセッサが従っている IEEE 754 という規格に照らして、これは完全に正しい答えです。しかも「たまたま最後の桁がずれた」のではなく、後で見るようにちょうど真ん中でのタイブレークが上に転んだ結果という、規格の細部まで説明できる必然的な現象です。

この現象を「浮動小数点あるある」で済ませてしまうと、痛い目に遭います。実際に起きた/起きうる事故を挙げます。

• 分散が負になる。$V[X] = E[X^2] – (E[X])^2$ という高校で習う公式をそのままコードにすると、平均が大きいデータ（センサーのオフセットが乗った信号、UNIX時刻、緯度経度など）で分散が負の値になります。本記事では実際に $-4.0$ という値を出します。標準偏差を取ろうとして nan が出る典型パターンです。
• 根が見つからない。多項式を展開形で評価してニュートン法や二分法に渡すと、真の根の周りで関数値が完全なノイズになり、根を $\pm 0.07$ 程度でしか特定できなくなります。反復回数を増やしても改善しません。
• 長時間積分が発散する。1万ステップ、100万ステップと足し込む数値積分・モンテカルロ法・オンライン統計量では、1回あたり $10^{-16}$ の誤差が積み上がって効いてきます。

応用先はもっと広く、金融計算（円未満の丸めが総額に効く）、GNSS測位（$10^7$ メートルの座標に対してセンチメートルの精度を要求する）、機械学習の損失計算（log(1+x) や softmax の数値安定化）、CFDや軌道伝播の長時間シミュレーションなど、「大きな数と小さな数を混ぜて扱う」場面のすべてに関わります。逆に言えば、IEEE 754の仕組みを一度きちんと理解しておけば、これらの落とし穴は設計段階で回避できるようになります。

この記事の全体像を1枚にまとめたのが上の図です。上段の「人間の数直線」では数がすきまなく連続に並んでいるので $0.1+0.2$ はぴったり $0.3$ になりますが、下段の「計算機が持てる数」は $0.3$ の近くでも $5.55 \times 10^{-17}$ 間隔のとびとびの点しかありません。しかも $0.1$ と $0.2$ を厳密に足した値は、その点と点のちょうど真ん中に落ちてしまう。どちらに倒すかを決める規則が「同点はきりのよい側へ」であり、その結果 $0.3$ とは反対側の点が選ばれる——これが本記事で最後の1ビットまで追いかける現象です。

本記事の内容

• なぜ $0.1$ が2進法で正確に書けないのかを、実際のビット列で確認する
• IEEE 754 binary64（float64）のビット配置から、機械イプシロン $2^{-52}$ と丸め誤差の上界 $u = 2^{-53}$ を導出する
• 引き算の条件数を導き、桁落ち（catastrophic cancellation）がなぜ致命的かを定量化する
• 二次方程式の解の公式・分散公式 $E[X^2]-(E[X])^2$ で桁落ちを実演し、有理化・Welfordの逐次更新という数値的に安定な代替式を導出する
• 総和の誤差蓄積とカハンの補正総和、np.float32 との比較、== ではなく np.isclose を使う理由

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

• 条件数と数値的安定性の理論 — 本記事の桁落ちの議論は、条件数という共通言語の上に乗っています
• 【Python】数値計算で微分可能な関数の微分値を求める — 刻み幅 $h$ を小さくしすぎると精度が悪化する、という現象の正体が本記事の丸め誤差です
• Numpyのrandomで生成できる乱数のまとめ — 実験用データの生成に使います

なぜ $0.1$ は正確に表せないのか — 2進法の「割り切れなさ」

まず、コンピュータを一旦忘れて、私たちが普段使っている10進法の話から始めます。10進法で $1/3$ を書こうとすると $0.3333\ldots$ と無限に続きます。有限の桁数しか書けない紙の上では、$1/3$ は原理的に正確に書けません。$0.3333$ と書いた瞬間、それは $1/3$ ではない別の数です。

同じことが2進法でも起きます。ただし「割り切れない分数」の顔ぶれが変わります。10進法では分母が $2$ と $5$ の積だけでできている分数（$1/2, 1/4, 1/5, 1/20, \ldots$）が有限桁で書けました。2進法で有限桁になるのは、既約分数にしたときの分母が2のべき乗である分数だけです。

$0.1 = \dfrac{1}{10} = \dfrac{1}{2 \times 5}$ の分母には $5$ が残っています。だから2進法では循環小数になります。実際に計算してみましょう。小数部を2倍して整数部を取り出す、を繰り返します。

$$ \begin{align} 0.1 \times 2 &= 0.2 \quad \to \quad 0 \\ 0.2 \times 2 &= 0.4 \quad \to \quad 0 \\ 0.4 \times 2 &= 0.8 \quad \to \quad 0 \\ 0.8 \times 2 &= 1.6 \quad \to \quad 1 \ (\text{残り } 0.6) \\ 0.6 \times 2 &= 1.2 \quad \to \quad 1 \ (\text{残り } 0.2) \end{align} $$

ここで小数部が $0.2$ に戻りました。2ステップ目とまったく同じ状態です。つまりこれ以降は $0011$ の4桁が永遠に繰り返されます。

$$ 0.1_{(10)} = 0.0001100110011001100\ldots_{(2)} = 0.0\overline{0011}_{(2)} $$

左のグラフは、いま手で追った「2倍して整数部を取り出す」操作を12回ぶん描いたものです。1回目の残りと5回目の残りがどちらも $0.2$ に戻っており（緑の矢印）、そこから先はまったく同じ状態を巡回するだけなので、取り出されるビット $0011$ が永遠に繰り返されることが目で確認できます。右のパネルは各数を2進小数として24桁ぶん並べたもので、$0.5, 0.25, 0.125$ は最初の1桁だけ青（=1）でそこから先はすべて灰（=0）＝有限桁で書き切れているのに対し、$0.1, 0.2, 0.3$ は青と灰の縞模様が端まで途切れず続いています。分母に $5$ が残るかどうかという1点だけで、有限と無限がきれいに分かれています。

これが「$0.1$ は2進法で正確に表せない」ということの中身です。10進法で $1/3$ が書けないのと、まったく同じ理由です。$0.5$ や $0.25$ や $0.125$ は分母が2のべき乗なので正確に表せますが、$0.1, 0.2, 0.3, 0.7$ といった日常でいちばん使う数がことごとく循環小数になってしまう。ここが、10進で暮らす人間と2進で計算する機械のすれ違いの出発点です。

では、この無限に続くビット列を、有限のメモリにどう押し込めるのでしょうか。次はfloat64の中身、64個のビットの割り当てを見ていきます。

IEEE 754 binary64 のビット配置

Pythonの float、C言語の double、NumPyの np.float64 は、すべてIEEE 754の binary64 という同じ形式です。64ビットを次のように3つに割り当てます。

これは科学表記 $1.2345 \times 10^{6}$ の2進版だと思ってください。$1.2345$ にあたるのが仮数、$10^6$ の $6$ にあたるのが指数です。正規化数（$1 \le E \le 2046$）の場合、表す値は次式で決まります。

$$ \begin{equation} x = (-1)^s \times \left(1 + \frac{M}{2^{52}}\right) \times 2^{E – 1023} \end{equation} $$

この式には2つの工夫が仕込まれています。順に説明します。

工夫1: 暗黙の1（ケチ表現）。仮数の先頭に $1 + \cdots$ と書いてあるのは、2進法の科学表記では先頭桁が必ず $1$ になるからです。10進法では $0.0012 = 1.2 \times 10^{-3}$ のように先頭を $1$〜$9$ に正規化しますが、2進法では $0$ でない桁は $1$ しかないので、正規化すると先頭は必ず $1$ です。分かりきっている $1$ をわざわざ記憶する必要はありません。この「ケチ」のおかげで、52ビットの記憶で 53ビット分の精度が得られます。

工夫2: バイアス付き指数。指数は負にもなり得ますが、符号ビットをもう1つ用意する代わりに、実際の指数に $1023$ を足した値 $E$ を格納します。$E$ が符号なし整数として大小比較できるので、浮動小数点数の大小比較を整数比較で代用できる（正の数の場合）という利点があります。実際の指数は $e = E – 1023$ で、正規化数では $-1022 \le e \le 1023$ の範囲を取ります。

$E = 0$ と $E = 2047$ は特別扱いです。$E=0$ は $0$ と非正規化数（subnormal、$2^{-1074}$ まで表せる極小の数）、$E=2047$ は $\pm\infty$ と nan に予約されています。

具体例で確かめます。$1.0$ は $1.0 \times 2^0$ なので $M = 0$、$e = 0$ すなわち $E = 1023 = 01111111111_{(2)}$ です。$2.0$ なら $E = 1024$、$0.5$ なら $E = 1022$。指数部だけが1ずつ動いていることがわかります。

では $0.1$ はどうなるでしょうか。先ほど求めた $0.0\overline{0011}_{(2)}$ を正規化します。最初の $1$ は小数第4位なので、

$$ 0.1 = 1.\overline{1001}_{(2)} \times 2^{-4} $$

となり、$e = -4$、$E = 1019 = 01111111011_{(2)}$ です。仮数部には $\overline{1001}$ の繰り返しを52ビットぶん詰め込みます。$52 = 4 \times 13$ なので、ちょうど $1001$ が13回並んで

$$ M_{\text{切り捨て}} = \underbrace{1001\,1001\,\cdots\,1001}_{13 \text{ 回}} $$

……となるはずですが、実際にPythonで見ると仮数部の末尾は $\ldots 1010$ になっています。

$$ \texttt{0 01111111011 1001100110011001100110011001100110011001100110011010} $$

末尾が $1001$ ではなく $1010$、つまり切り上げられています。53ビット目以降が $1001\ldots$ で、捨てる部分が残り全体の半分より大きいからです。この「捨てる部分をどう処理するか」が丸めであり、次節の主題です。

結果として、Pythonの 0.1 が実際に保持している値は $0.1$ そのものではなく、

$$ 0.1000000000000000055511151231257827021181583404541015625 $$

という、$0.1$ よりわずかに大きい数です。from decimal import Decimal; print(Decimal(0.1)) と打てば誰でも確認できます。同様に 0.2 は

$$ 0.200000000000000011102230246251565404236316680908203125 $$

で、これも $0.2$ より少し大きい。この2つが後で効いてきます。

上段が64ビットの区画割り、下段が実際に格納されているビット列（色つきが1、灰が0）です。$1.0$ と $0.5$ は仮数部が全ビット $0$ で、違うのは指数部の値 $E$ が $1023$ か $1022$ かだけ——「$2$ のべき乗はきれいに表せる」が一目でわかります。対照的に $0.1$ と $0.2$ は仮数部が右端まで縞模様で埋まっており、しかも両者の仮数部は完全に同一で指数部だけが $1019$ と $1020$ で1違う（$0.2 = 2 \times 0.1$ の反映）。そして最下段の $0.3$ と $0.1+0.2$ を比べると、黒枠で囲った末尾の数ビットだけが食い違っています。この1ビットのずれの正体を、次節以降で丸めの規則から説明します。

ここまでで「64ビットに何がどう入っているか」がわかりました。次は「入りきらない部分をどうするか」——丸めのルールと、そこから出てくる誤差の上限を導きます。

丸めと機械イプシロン $2^{-52}$、丸め誤差の上界 $u = 2^{-53}$

実数直線の上に、float64で表せる数を点として打つとどう見えるでしょうか。等間隔ではありません。2のべき乗の区間ごとに、点の間隔が2倍ずつ広がるという、対数目盛の物差しのような並び方になります。

理由は式(1)から明らかです。指数 $e$ を固定すると、仮数 $M$ が $1$ 増えるごとに値は $2^{e-52}$ ずつ増えます。つまり区間 $[2^e, 2^{e+1})$ の中では、表せる数が $2^{e-52}$ 間隔で $2^{52}$ 個並んでいる。この間隔を ULP（Unit in the Last Place、最終桁の重み）と呼びます。

$$ \begin{equation} \mathrm{ulp}(x) = 2^{e-52}, \qquad 2^e \le |x| < 2^{e+1} \end{equation} $$

具体的な数字で見ると迫力があります。

いちばん下の行に注目してください。$10^{16}$ 付近では、隣り合う浮動小数点数の間隔が 2 です。つまり 1e16 + 1 を計算すると答えは 1e16 のまま変わりません（実際にそうなります）。整数が整数として扱えなくなる境界は $2^{53} \approx 9.007 \times 10^{15}$ で、これを超えるとすべての整数を表せなくなります。JavaScriptの Number.MAX_SAFE_INTEGER が $2^{53}-1$ なのはこのためです。

両対数でプロットすると、ULPの実測値（青）は傾き $1$ の直線 $x \cdot 2^{-52}$（黒破線）にぴったり沿って伸びています。細かく見ると青線が階段状にギザギザしているのは、区間 $[2^e, 2^{e+1})$ の中では間隔が一定で、$2$ のべき乗をまたぐ瞬間だけ2倍に跳ぶからです。この直線が「間隔 $=1$」の赤い水平線と交わるのがちょうど $2^{53} \approx 9.0 \times 10^{15}$ で、ここから右では連続する整数を区別できなくなります。「大きな整数IDをfloatで持ってはいけない」という定番の注意が、交点の位置として読み取れます。

丸めのルール。実数 $x$ が2つの浮動小数点数の間に落ちたとき、IEEE 754の既定は round-to-nearest, ties-to-even（最近接丸め、同点なら仮数の最終ビットが偶数の側へ）です。「同点なら偶数へ」というのは統計的な偏りを避けるための工夫で、いつも切り上げにすると誤差が一方向に溜まってしまうからです。

このルールから、丸め誤差の大きさを見積もります。$x$ を実数、$\mathrm{fl}(x)$ をそれに最も近い浮動小数点数とします。最近接丸めなので、$x$ と $\mathrm{fl}(x)$ の距離は間隔の半分以下です。

$$ |\mathrm{fl}(x) – x| \le \frac{1}{2}\,\mathrm{ulp}(x) = \frac{1}{2} \cdot 2^{e-52} = 2^{e-53} $$

一方、$x$ が区間 $[2^e, 2^{e+1})$ にいることから $|x| \ge 2^e$ です。この2つを組み合わせて相対誤差を評価します。分母を $|x| \ge 2^e$ で下から押さえると分数は大きくなるので、

$$ \frac{|\mathrm{fl}(x) – x|}{|x|} \le \frac{2^{e-53}}{2^{e}} = 2^{-53} $$

$e$ がきれいに約分されて消えました。ここが重要です。相対誤差の上界は数の大きさによらず一定であり、その値は

$$ \begin{equation} u = 2^{-53} \approx 1.11 \times 10^{-16} \end{equation} $$

です。この $u$ を unit roundoff（丸めの単位）と呼びます。同じことを別の形に書くと、任意の実数 $x$ について

$$ \begin{equation} \mathrm{fl}(x) = x(1 + \delta), \qquad |\delta| \le u = 2^{-53} \end{equation} $$

と書けます。この $(1+\delta)$ 表記は、以降の誤差解析で何度も使う道具です。

機械イプシロンとの違い。ここでよく混乱が起きます。np.finfo(np.float64).eps を打つと $2.22 \times 10^{-16} = 2^{-52}$ が返ってきて、$u = 2^{-53}$ の2倍です。両者の定義が違うからです。

• 機械イプシロン $\varepsilon = 2^{-52}$ は、$\mathrm{ulp}(1)$、すなわち $1$ の次に大きい浮動小数点数と $1$ の差です。「$1$ の周りでの分解能」を表します。
• unit roundoff $u = \varepsilon/2 = 2^{-53}$ は、丸めによる相対誤差の上界です。

「間隔」と「間隔の半分」の違いだと覚えておけば混乱しません。誤差の上限を見積もるときは $u = 2^{-53}$ を、「どこまで細かい差を区別できるか」を考えるときは $\varepsilon = 2^{-52}$ を使います。ちなみに $2^{53} \approx 9.0 \times 10^{15}$ なので、float64が保持できる10進の有効桁数は $\log_{10} 2^{53} \approx 15.95$、つまり約16桁です。17桁目以降は原理的に存在しません。

丸めのルールが厳密に決まっていることがわかりました。ではこの厳密なルールを、冒頭の 0.1 + 0.2 に適用してみましょう。答えがなぜ $0.30000000000000004$ という中途半端な数になるのか、最後の1ビットまで追跡できます。

$0.1 + 0.2$ が $0.3$ にならない理由を最後の1ビットまで追う

計算は3段階に分かれます。①$0.1$ を丸める、②$0.2$ を丸める、③その2つの和を丸める。誤差が入る余地は3か所です。

①と②は前節で見た通りで、どちらも真の値より少し大きい数になりました。この2つは無限精度で足せます（IEEE 754の加算は「無限精度で足してから丸める」と定義されています）。

$$ \begin{align} \mathrm{fl}(0.1) + \mathrm{fl}(0.2) &= 0.1000000000000000055511151231257827\ldots \\ &\quad + 0.2000000000000000111022302462515654\ldots \\ &= 0.3000000000000000166533453693773481\ldots \end{align} $$

③この和を、$0.3$ 付近の浮動小数点数に丸めます。$0.3$ は区間 $[0.25, 0.5) = [2^{-2}, 2^{-1})$ にいるので、この付近のULPは $2^{-54} \approx 5.5511 \times 10^{-17}$ です。候補となる隣り合う2つの浮動小数点数は

$$ \begin{align} \text{下側} &= \mathrm{fl}(0.3) = 0.2999999999999999888977697537484346\ldots \\ \text{上側} &= 0.3000000000000000444089209850062616\ldots \end{align} $$

先ほどの厳密な和 $0.30000000000000001665\ldots$ は、この2つのどちらに近いでしょうか。差を取ります。

$$ \begin{align} \text{和} – \text{下側} &= 2.77555756156289135\ldots \times 10^{-17} \\ \text{上側} – \text{和} &= 2.77555756156289135\ldots \times 10^{-17} \end{align} $$

完全に同点です。$2.7755\ldots \times 10^{-17}$ はULPのちょうど半分 $2^{-55}$ に一致します。ここで ties-to-even の出番です。下側の仮数の末尾4ビットは $0011$（最終ビットが $1$ で奇数）、上側は $0100$（最終ビットが $0$ で偶数）。規格は偶数側を選ぶので、上側が勝ちます。

$$ \mathrm{fl}(0.1) \oplus \mathrm{fl}(0.2) = 0.30000000000000004\ldots \ne \mathrm{fl}(0.3) $$

これが 0.1 + 0.2 == 0.3 が False になる、最後の1ビットまでの理由です。「なんとなく誤差が出た」のではなく、きっかり真ん中に落ちた和が、偶数丸めのルールで $0.3$ とは反対側に転んだ。1ビット、値にして $5.55 \times 10^{-17}$ のずれです。

上段の数直線で、星印（厳密な和）から左右の候補までの距離がどちらも $2.7756 \times 10^{-17}$ と表示されています。四捨五入で言えば「ちょうど $5$」の状況で、最近接丸めだけでは行き先が決まりません。下段はそこで効く決着ルールで、両候補の仮数部の末尾8ビットを並べたものです。上側の候補は最終ビットが $0$（偶数）、下側＝計算機の $0.3$ は最終ビットが $1$（奇数）。ties-to-even は偶数側を選ぶので上側が採用され、答えは $\mathrm{fl}(0.3)$ から1ビットぶんだけ外れます。

面白いことに、同じことを 0.1 + 0.4 でやると 0.5 と一致しますし、0.1 + 0.2 + 0.3 と 0.3 + 0.2 + 0.1 は違う答えになります。浮動小数点の加算は結合法則を満たしません。$(a+b)+c \ne a+(b+c)$ です。これは並列計算で「足す順番が実行のたびに変わる」と結果が再現しなくなる、という実務上の問題に直結します。

さて、1回の演算での誤差が $u = 2^{-53}$ 以下だとわかりました。$10^{-16}$ の誤差なんて工学的には無視できそうに思えます。ところが、この $10^{-16}$ が一気に $10^{0}$ まで拡大される演算が存在します。それが引き算です。

演算の誤差モデルと、引き算の条件数

IEEE 754は、四則演算と平方根について次を保証しています。「無限精度で計算した正確な結果を、最も近い浮動小数点数に丸めたものを返す」（正しく丸められた演算）。これを式(4)と組み合わせると、実際の演算 $\oplus, \ominus, \otimes, \oslash$ について

$$ \begin{equation} \mathrm{fl}(a \circ b) = (a \circ b)(1 + \delta), \qquad |\delta| \le u \end{equation} $$

というモデルが立ちます。1回の演算で持ち込まれる相対誤差は高々 $u$、という単純明快な保証です。

ここで大事なのは、この保証が入力がすでに正確である場合の話だという点です。実際には $a$ も $b$ も、それ以前の丸めで少し汚れています。$\tilde{a} = a(1+\delta_a)$、$\tilde{b} = b(1+\delta_b)$（$|\delta_a|, |\delta_b| \le u$）とします。この汚れた入力に演算をかけると何が起きるか。掛け算と引き算で比べます。

掛け算の場合。$\tilde{a}\tilde{b} = ab(1+\delta_a)(1+\delta_b) \approx ab(1 + \delta_a + \delta_b)$ なので、相対誤差は高々 $2u$ です。入力の相対誤差がそのまま足し算されるだけで、増幅しません。割り算も同様です。

引き算の場合。ここが問題です。

$$ \tilde{a} – \tilde{b} = a(1+\delta_a) – b(1+\delta_b) = (a – b) + (a\delta_a – b\delta_b) $$

誤差の絶対値は $|a\delta_a – b\delta_b| \le (|a| + |b|)u$ で押さえられます。これを結果 $a-b$ で割って相対誤差にすると、

$$ \begin{equation} \frac{|(\tilde{a}-\tilde{b}) – (a-b)|}{|a-b|} \le \frac{|a| + |b|}{|a-b|} \, u \equiv \kappa \cdot u \end{equation} $$

この

$$ \begin{equation} \kappa = \frac{|a| + |b|}{|a – b|} \end{equation} $$

が引き算の条件数です。分子は入力の大きさ、分母は結果の大きさ。$a$ と $b$ が近いと分母が小さくなり、$\kappa$ が爆発します。

数字を入れます。$a = 1.0000001$、$b = 1.0000000$ なら $\kappa \approx 2/10^{-7} = 2 \times 10^{7}$ です。入力の相対誤差 $10^{-16}$ が、結果では $2 \times 10^{-9}$ にまで拡大されます。$a$ と $b$ が16桁一致するところまで近ければ $\kappa \approx 10^{16}$ となり、$\kappa u \approx 1$、つまり結果に有効数字が1桁も残りません。

この現象を 桁落ち（catastrophic cancellation、破滅的な打ち消し）と呼びます。「近い数どうしの引き算をしてはいけない」という数値計算の第一戒律の根拠がこれです。

注意すべきは、引き算そのものは何も悪いことをしていないという点です。式(5)により、引き算の演算自体が持ち込む誤差は高々 $u$ です。悪いのは、$a$ と $b$ がすでに持っていた誤差が、結果の小ささに対して相対的に巨大になること。桁落ちは「新しい誤差を作る」のではなく「既にあった誤差を白日の下に晒す」操作なのです。だからこそ、対策は「引き算をしない式に書き換える」ことになります。

条件数という物差しが手に入りました。次はこれを実際の計算式に当てはめ、教科書の公式が数値的には使い物にならない場面を見ていきます。

桁落ちの実演①: $\sqrt{x+1} – \sqrt{x}$ と有理化

$f(x) = \sqrt{x+1} – \sqrt{x}$ を考えます。$x$ が大きくなると $\sqrt{x+1}$ と $\sqrt{x}$ は限りなく近づくので、絵に描いたような桁落ちが起きます。

条件数を見積もります。$a = \sqrt{x+1}$、$b = \sqrt{x}$ として式(7)に代入し、$x$ が大きいとき $f(x) \approx 1/(2\sqrt{x})$ であることを使うと、

$$ \kappa = \frac{\sqrt{x+1} + \sqrt{x}}{\sqrt{x+1} – \sqrt{x}} \approx \frac{2\sqrt{x}}{1/(2\sqrt{x})} = 4x $$

$\kappa \approx 4x$ なので、$x = 10^{14}$ では $\kappa \approx 4 \times 10^{14}$、$\kappa u \approx 4.4 \times 10^{-2}$。つまり相対誤差が数%に達する見積もりです。

対策は中学校で習った有理化です。分子と分母に $\sqrt{x+1}+\sqrt{x}$ を掛けます。

$$ \begin{align} \sqrt{x+1} – \sqrt{x} &= \frac{(\sqrt{x+1} – \sqrt{x})(\sqrt{x+1} + \sqrt{x})}{\sqrt{x+1} + \sqrt{x}} \end{align} $$

分子は和と差の積なので $(\sqrt{x+1})^2 – (\sqrt{x})^2 = (x+1) – x = 1$ と、きれいに $1$ になります。したがって

$$ \begin{equation} \sqrt{x+1} – \sqrt{x} = \frac{1}{\sqrt{x+1} + \sqrt{x}} \end{equation} $$

数学的には完全に同じ式ですが、右辺には引き算が1つもありません。足し算と割り算だけなので、式(5)より相対誤差は高々数 $u$ で収まります。実際に測った結果が次の表です（真値は多倍長演算で計算）。

この表を、横軸を条件数 $\kappa$ にしてグラフにしたのが上の図です。素直な式（青丸）は理論の上界 $\kappa u$（黒破線）とほぼ平行に、傾き $1$ で上がっていきます。式(6)の不等式が「上界」としてだけでなく実際の誤差の大きさの予言として機能していることが読み取れます。一方で有理化した式（緑四角）は、条件数が $10^{17}$ まで増えても $10^{-16}$ 台の水平線に張りついたままです。同じ数学的な値を計算しているのに、青と緑の縦の隔たりは右端で16桁ぶんに達しています。

読み取れることが3つあります。第一に、素直な式の誤差は理論の見積もり $\kappa u$ とほぼ同じ勢いで増えていて、$\kappa$ が10倍になれば誤差も10倍になります。式(6)の不等式が実測でそのまま成り立っている、ということです。第二に、有理化した式の誤差は $x$ をどれだけ大きくしても $10^{-16}$ 台に張り付いたままです。第三に、$x = 10^{16}$ では素直な式が $0$ を返します。$x+1$ を計算した時点で $1$ が丸め落ちて $x$ になってしまうからで、相対誤差 $100\%$、有効数字ゼロです。

同じ構造の落とし穴が標準ライブラリにも用意されています。$e^x – 1$ を $x \approx 0$ で計算すると桁落ちしますが、math.expm1(x) を使えば回避できます。$x = 10^{-10}$ で実測すると、math.exp(x)-1 は相対誤差 $8.3 \times 10^{-8}$、math.expm1(x) はほぼ完全です。$\log(1+x)$ も同じで、math.log(1 + 1e-16) は $0.0$ を返しますが math.log1p(1e-16) は $10^{-16}$ を返します。「$\text{expm1}$ / $\text{log1p}$ が標準ライブラリにわざわざある理由」はここにあります。機械学習で交差エントロピー損失を実装するときに log1p や logsumexp を使うのも、まったく同じ動機です。

有理化という書き換えの威力がわかりました。次は、もっと身近な公式——二次方程式の解の公式——で同じことをやってみます。

桁落ちの実演②: 二次方程式の解の公式

中学校で習った $ax^2+bx+c=0$ の解の公式

$$ x = \frac{-b \pm \sqrt{b^2 – 4ac}}{2a} $$

は、そのままコードにすると危険です。$b^2 \gg 4ac$ のとき $\sqrt{b^2-4ac} \approx |b|$ となるため、$-b$ と $\pm\sqrt{\cdot}$ の符号が逆になる側で桁落ちが起きます。

$a=1$、$b=10^8$、$c=1$ で試します。$b^2 = 10^{16}$ に対して $4ac=4$ なので、判別式はほとんど $b^2$ です。$\sqrt{10^{16}-4} \approx 99999999.99999998$ で、これと $-b = -10^8$ を足す（$+\sqrt{\cdot}$ の側）と、8桁が丸ごと打ち消されます。実測すると

• 素直な公式: $x_1 = -7.450580596923828 \times 10^{-9}$
• 真値（多倍長）: $x_1 = -1.0000000000000001 \times 10^{-8}$

相対誤差 25% です。有効数字は1桁も合っていません。

左のグラフは $b$ を $10$ から $10^{11}$ まで振ったときの、小さいほうの根の相対誤差です。教科書どおりの公式（青）は $b$ に比例して誤差が増え、$b = 10^{8}$（星印）で $25\%$、$b \gtrsim 10^{9}$ では相対誤差 $1$ に張りついて有効数字がゼロになります。有理化した式（緑）は $b$ によらず $10^{-16}$ 台のままです。右の棒グラフが桁落ちの現場そのもので、$-b$ と $+\sqrt{b^2-4ac}$ というどちらも $10^8$ の大きさを持つ2数を足した結果が $1.49 \times 10^{-8}$ しかない。$10^8$ から $10^{-8}$ へ、16桁ぶんが打ち消されて消えているので、float64が持つ約16桁の情報が残らないのは当然だとわかります。

対策は先ほどと同じ有理化です。分子・分母に $(-b – \sqrt{b^2-4ac})$ を掛けます。

$$ x_1 = \frac{-b + \sqrt{b^2-4ac}}{2a} \cdot \frac{-b – \sqrt{b^2-4ac}}{-b – \sqrt{b^2-4ac}} $$

分子は和と差の積なので $b^2 – (b^2 – 4ac) = 4ac$ になります。

$$ x_1 = \frac{4ac}{2a\left(-b – \sqrt{b^2-4ac}\right)} = \frac{2c}{-b – \sqrt{b^2-4ac}} $$

$b > 0$ なら分母は $-b – |b|\cdot(\text{ほぼ}1) \approx -2b$ となり、打ち消しが起きません。この式で計算すると $x_1 = -1.0 \times 10^{-8}$、相対誤差 $7.9 \times 10^{-17}$ と、真値と機械精度で一致します。

一般的な実装方針は次の通りです。まず $q = -\dfrac{1}{2}\left(b + \mathrm{sign}(b)\sqrt{b^2-4ac}\right)$ を計算します。$\mathrm{sign}(b)$ を掛けているので、この足し算は必ず同符号どうしの加算になり、桁落ちが起きません。そのうえで

$$ x_1 = \frac{q}{a}, \qquad x_2 = \frac{c}{q} $$

と2つの根を求めます（$x_2$ には解と係数の関係 $x_1 x_2 = c/a$ を使っています）。この形なら、$b$ の符号がどちらでも安全です。数値計算ライブラリの二次方程式ソルバはほぼこの実装になっています。

「教科書の公式をそのまま書いてはいけない」ことが2例で確認できました。次は、統計を扱う人がほぼ全員一度は踏む、もっとも有名な地雷を見ます。

分散の計算式 $E[X^2] – (E[X])^2$ はなぜ壊れるか

分散の定義は

$$ V[X] = \frac{1}{n}\sum_{i=1}^{n}(x_i – \bar{x})^2 $$

ですが、これを展開して得られる

$$ \begin{equation} V[X] = \frac{1}{n}\sum_{i=1}^{n} x_i^2 – \bar{x}^2 = E[X^2] – (E[X])^2 \end{equation} $$

という「1パスで計算できる便利な公式」がよく使われます。データを1回なめるだけで $\sum x_i$ と $\sum x_i^2$ が両方たまるので、ストリーミング処理でも使えて魅力的に見えます。

ところが、この式の最後の演算は引き算です。しかも、平均 $\bar{x}$ が大きく分散が小さいデータでは、$E[X^2]$ と $\bar{x}^2$ がほとんど同じ値になります。式(7)で条件数を評価しましょう。$a = E[X^2]$、$b = \bar{x}^2$、$a – b = V[X]$ なので、$V[X] \ll \bar{x}^2$ のとき

$$ \begin{equation} \kappa = \frac{E[X^2] + \bar{x}^2}{V[X]} \approx \frac{2\bar{x}^2}{V[X]} \end{equation} $$

平均の2乗と分散の比が、そのまま条件数になります。平均が大きいほど、分散が小さいほど危険です。

実験してみます。標準正規乱数（分散 $\approx 1$）に一定のオフセットを加えたデータ $10{,}000$ 点で、3通りの方法で分散を計算しました。真値は多倍長演算で求めています。

オフセット $10^8$ の行で、展開公式が $-4.0$ を返しています。分散は定義上ゼロ以上のはずなのに、負です。標準偏差を取ろうとすれば nan になります。オフセット $10^9$ では $128.0$、真値 $1.03$ の124倍です。しかもこれは異常なデータではありません。UNIX時刻（$1.7 \times 10^9$）、気圧のPa表示、ADCの生カウント値など、$10^8$〜$10^9$ 程度のオフセットは実務でごく普通に現れます。

条件数の見積もりとも合っています。オフセット $10^6$ で $\kappa u \approx 2\times10^{12} \times 1.1\times10^{-16} = 2.2\times10^{-4}$、実測 $2.3 \times 10^{-4}$。ぴたりです。オフセット $10^8$ では $\kappa u \approx 2$ で「相対誤差が $1$ を超える」、つまり答えが完全に無意味になる領域に入っています。実際に負になったのはその通りの帰結です。

オフセットを横軸に振ると、展開公式（青丸）の相対誤差が条件数から見積もった上界（黒破線）にぴったり重なりながら右上がりに増え、オフセット $10^{8}$ あたりで相対誤差 $1$ の赤い線を突き抜けているのが見えます。赤いバツ印が、まさに分散が負の値になった点です。それに対して二パス法（緑）はオフセットをどれだけ大きくしても $10^{-16}$ 台の水平線から動かず、Welford法（紫）はオフセットとともに緩やかに悪化するものの、右端でも $10^{-7}$ 前後にとどまっています。「同じ分散を求める3つの式」の縦の隔たりが、右端では15桁以上に開くことがこの1枚に凝縮されています。

対策1: 二パス法

いちばん確実なのは、定義式をそのまま使うことです。

1. まず全データをなめて $\bar{x} = \frac{1}{n}\sum x_i$ を求める
2. もう一度なめて $\frac{1}{n}\sum (x_i – \bar{x})^2$ を求める

引き算 $x_i – \bar{x}$ はデータ点ごとに行われますが、$x_i$ と $\bar{x}$ が近いのは当たり前で、その差こそが求めたい量です。ここでの桁落ちは「小さい答えを作るために大きい数を引く」のではなく「答えそのものが小さい」だけなので、相対誤差は増幅しません。表を見ると、二パス法の誤差はオフセットによらず $10^{-16}$ 台に収まっています。NumPyの np.var はこの二パス法（正確にはその改良版）を使っているので、np.var を使う限り安全です。

対策2: Welfordの逐次更新

二パス法の弱点は、データを2回読む必要があることです。ストリーミングデータや巨大ファイルでは困ります。そこで使うのが Welford のオンラインアルゴリズムです。データが1点届くたびに、平均と平方和を更新していきます。

平均の更新式を導きます。$m_n$ を $n$ 個目までの平均とすると、

$$ m_n = \frac{1}{n}\sum_{i=1}^{n} x_i = \frac{(n-1)m_{n-1} + x_n}{n} $$

右辺を $m_{n-1}$ を軸に整理します。分子に $-m_{n-1} + m_{n-1}$ を仕込んで $(n-1)m_{n-1} + x_n = n\,m_{n-1} + (x_n – m_{n-1})$ と書き直すと、

$$ \begin{equation} m_n = m_{n-1} + \frac{x_n – m_{n-1}}{n} \end{equation} $$

「今の平均に、新しい点とのずれを $1/n$ だけ効かせる」という形です。ずれ $\delta = x_n – m_{n-1}$ は小さい量なので、$m_{n-1}$ に対する微小な補正として加わります。桁落ちの心配がありません。

平方和の更新式を導きます。$M_{2,n} = \sum_{i=1}^{n}(x_i – m_n)^2$ とおきます。目標は $M_{2,n}$ を $M_{2,n-1}$ から作ることです。差を書き出します。

$$ M_{2,n} – M_{2,n-1} = (x_n – m_n)^2 + \sum_{i=1}^{n-1}\left[(x_i – m_n)^2 – (x_i – m_{n-1})^2\right] $$

右辺の和の中身に、和と差の積 $A^2 – B^2 = (A-B)(A+B)$ を使います。$A = x_i – m_n$、$B = x_i – m_{n-1}$ とすれば $A – B = m_{n-1} – m_n$、$A + B = 2x_i – m_n – m_{n-1}$ です。$A-B$ は $i$ によらない定数なので和の外に出せて、

$$ \sum_{i=1}^{n-1}(A-B)(A+B) = (m_{n-1}-m_n)\sum_{i=1}^{n-1}\left(2x_i – m_n – m_{n-1}\right) $$

$\sum_{i=1}^{n-1} x_i = (n-1)m_{n-1}$ を代入すると、和の部分は $2(n-1)m_{n-1} – (n-1)(m_n + m_{n-1}) = (n-1)(m_{n-1} – m_n)$ になります。したがって

$$ \sum_{i=1}^{n-1}\left[(x_i-m_n)^2 – (x_i-m_{n-1})^2\right] = (n-1)(m_{n-1}-m_n)^2 $$

次に $\delta = x_n – m_{n-1}$ を導入して整理します。式(11)より $m_n – m_{n-1} = \delta/n$ なので上の項は $(n-1)\delta^2/n^2$。また $x_n – m_n = \delta – \delta/n = \delta(n-1)/n$ なので $(x_n-m_n)^2 = \delta^2(n-1)^2/n^2$。2つを足すと

$$ M_{2,n} – M_{2,n-1} = \frac{\delta^2(n-1)^2}{n^2} + \frac{\delta^2(n-1)}{n^2} = \frac{\delta^2(n-1)\left[(n-1)+1\right]}{n^2} = \frac{\delta^2 (n-1)}{n} $$

最後にこれを、実装しやすい形に書き換えます。$x_n – m_n = \delta(n-1)/n$ だったので $\delta \cdot (x_n – m_n) = \delta^2(n-1)/n$ となり、上式の右辺とちょうど一致します。よって

$$ \begin{equation} M_{2,n} = M_{2,n-1} + (x_n – m_{n-1})(x_n – m_n) \end{equation} $$

これがWelfordの更新式です。更新前の平均とのずれと更新後の平均とのずれを掛けて足す、という覚えやすい形になりました。$\delta$ という「小さい量」だけを扱っているので、大きい数どうしの引き算が現れません。分散は $M_{2,n}/n$（不偏分散なら $M_{2,n}/(n-1)$）です。

表に戻ると、Welford法の誤差はオフセット $10^9$ でも $3.9\times10^{-8}$ で、展開公式の相対誤差 $123$ とは比べ物になりません。二パス法よりはやや大きいものの（$n$ に比例して誤差が蓄積するため）、データを1回しか読まずにこの精度というのが価値です。オンライン統計量、ストリーム処理、逐次的な異常検知など、実務での用途は非常に広いです。

分散が壊れる仕組みと、その2種類の直し方がわかりました。次は、引き算とは逆に「大きい数に小さい数を足すと小さいほうが消える」という現象と、その対策を見ます。

情報落ちと総和の誤差 — カハンの補正総和

1e16 + 1 を計算してみてください。答えは 1e16 です。$1$ が消えました。前に見た通り、$10^{16}$ 付近のULPは $2$ なので、$1$ は隣の浮動小数点数に届かず、丸めで消滅します。これを 情報落ち（swamping、絶対値の小さい数が飲み込まれる現象）と呼びます。

左のグラフは $s = 2^{k}$ に $1$ を足したとき「実際に増えた量」を測ったものです。$k \le 52$ までは律儀に $1$ 増えますが、$k = 53$ を境にぴたりと $0$ に落ちます。その付近の刻み幅（黒破線）が $2$ を超えた瞬間で、足したい $1$ が刻みの半分未満になると丸めで消えるという境目がはっきり出ています。右のグラフはより実務的で、初期値の異なる変数に $0.1$ を $1000$ 回足し込んだときの増分の相対誤差です。足し込み先が大きいほど誤差が比例して悪化し、初期値 $10^{15}$ 付近では増分がまるごと消えます。同じ足し算でも「どこに足すか」で精度が決まるわけです。

情報落ちが本当に困るのは、大量の足し算をするときです。$n$ 個の数を素直に順番に足していく状況を考えます。

$$ s_1 = x_1, \quad s_k = \mathrm{fl}(s_{k-1} + x_k) \quad (k = 2,\ldots,n) $$

各ステップで式(5)より $s_k = (s_{k-1} + x_k)(1+\delta_k)$、$|\delta_k| \le u$ です。ここで問題なのは、$k$ が進むと $s_{k-1}$ が足し込んだ全部の合計まで育っている一方、$x_k$ は一定サイズにとどまることです。育った合計に小さな数を足すたびに、その小さな数の下位ビットが少しずつ削られます。

誤差を積み上げると、最終的な総和の誤差は

$$ \begin{equation} \left|\hat{S} – S\right| \lesssim (n-1)\,u \sum_{i=1}^{n}|x_i| \end{equation} $$

で押さえられます。誤差が $n$ に比例する、これが素直な総和（naive summation）の性質です。$n = 10^6$ なら $nu \approx 10^{-10}$ で、有効桁が16桁から6桁ぶん減る計算になります。

実測します。$0.1$ を $n$ 個足す実験です（$0.1$ は前に見た通り真の $0.1$ よりわずかに大きいので、誤差が一方向に溜まる意地悪なケースです）。

素直な総和の誤差は $n$ が10倍になるごとに約100倍に増えています。式(13)の右辺は $\sum|x_i| = 0.1n$ を含むので、誤差の見積もりは $n \cdot u \cdot 0.1n \propto n^2$。すべての項が同符号で誤差が打ち消し合わないため、上界の $n^2$ がそのまま出ているわけです。$n=10^6$ で誤差 $1.3\times10^{-6}$、総和 $100000$ に対する相対誤差 $1.3\times10^{-11}$。倍精度を使っているのに有効数字が11桁しか残っていません。

素直な総和（青丸）は傾き $2$ の参照直線（黒破線）にぴったり重なっており、$n$ が10倍になると誤差が100倍になることが実測で確認できます。np.sum（緑三角）は明らかに緩い傾きで、$n = 10^{7}$ でも青との差は3〜4桁ぶんあります。これが pairwise summation の $O(u \log n)$ という性質です。そしてカハンの補正総和（紫）は $n = 10^{7}$ まで誤差がちょうど $0$ で、対数目盛の下限に張りついたまま動きません。3本の線の傾きの違いが、そのままアルゴリズムの質の違いになっています。

カハンの補正総和の原理

カハンの補正総和（Kahan compensated summation）は、この誤差を劇的に減らします。アイデアはとてもシンプルです——丸めで捨てられた下位ビットを覚えておいて、次の足し算で返す。

補正変数 $c$ に「前回捨てられた分」を保持します。1ステップは次の4行です。

$$ \begin{align} y &\leftarrow x_k – c && \text{(前回捨てた分を今回の項に返す)} \\ t &\leftarrow s + y && \text{(ここで丸めが起き、$y$ の下位が削られる)} \\ c &\leftarrow (t – s) – y && \text{(実際に加算された量 $-$ 加算したかった量 $=$ 捨てられた分)} \\ s &\leftarrow t \end{align} $$

3行目が発明の核心です。$t – s$ は「$y$ のうち実際に $s$ に反映された部分」を表します。浮動小数点の引き算 $t \ominus s$ は、$t$ と $s$ が近いので誤差なしで正確に計算されます（Sterbenzの補題: $s/2 \le t \le 2s$ なら $t-s$ は厳密）。そこから本来足したかった $y$ を引けば、残るのは足しそこねた分の符号を反転したものです。これを $c$ に入れて次のループで $x_{k+1}$ から引けば、捨てた分がきちんと戻ってきます。

c の符号に混乱しやすいので整理すると、$c$ は「$s$ が本来よりどれだけ多いか」を持ちます。だから1行目で $x_k$ から引くことで帳尻が合います。

理論的な誤差限界は

$$ \begin{equation} \left|\hat{S}_{\text{Kahan}} – S\right| \le \left(2u + O(nu^2)\right)\sum_{i=1}^{n}|x_i| \end{equation} $$

です。式(13)の $(n-1)u$ と比べてください。主要項から $n$ が消えています。何個足しても、誤差は「2回ぶんの丸め」程度で済む。表の実測でも、$n = 10^6$ までカハン加算の誤差は 厳密に $0$ でした。

NumPyの np.sum はなぜそこそこ速くて正確なのか

表の右列に注目してください。np.sum の誤差は素直な総和より4桁ほど小さく、しかしカハン加算ほど完璧ではありません。これは NumPy が pairwise summation（対分割総和）を使っているからです。配列を再帰的に半分に分け、小さいブロックまで分割してから足し上げます。加算の木の深さが $\log_2 n$ になるので、誤差は $O(u \log n)$ で済みます。$n$ に比例する素直な総和と比べて圧倒的に良く、しかも補正のための追加演算がないので素直な総和と同じ速さです。「特に何も考えずに np.sum を使えばそこそこ安全」なのは、この実装のおかげです。

さらに厳密さが必要なら、Python標準の math.fsum があります。これは Shewchukのアルゴリズムで、部分和を複数の浮動小数点数の組で厳密に保持し、常に正しく丸められた総和（真値に最も近いfloat64）を返します。速度は落ちますが、答えは一意に定まります。用途に応じて、

• 速度重視・十分な精度: np.sum（pairwise）
• 精度重視・ストリーミング可: カハン加算（またはより強力なNeumaier変種）
• 完全な正しさ: math.fsum

と使い分けるとよいでしょう。

総和の話が済んだところで、もう一つ実務で嫌な形で現れる増幅——反復法における「関数値のノイズ」を見ておきます。

ニュートン法・根の探索で誤差が増幅する仕組み

ニュートン法は $x_{k+1} = x_k – f(x_k)/f'(x_k)$ で根に近づく反復法です。理論的には2次収束、つまり有効桁数が毎回2倍になる強力な手法です。ところが、$f$ の評価に誤差があると、到達できる精度に天井ができます。

根 $x^\ast$ の近くで $f(x) \approx f'(x^\ast)(x – x^\ast)$ と線形近似できます。$f$ の計算に絶対誤差 $\epsilon_f$ が乗っているとすると、計算上の「$f$ がゼロに見える範囲」は

$$ \begin{equation} |x – x^\ast| \lesssim \frac{\epsilon_f}{|f'(x^\ast)|} \end{equation} $$

の幅を持ちます。この幅の中では、$f$ の符号すらランダムに変わります。どんな反復法も、この幅より細かく根を特定できません。$f’$ が小さい（根が重根に近い）ほど、また $f$ の評価誤差が大きいほど、天井は高くなります。

劇的な例を作ります。$p(x) = (x-2)^9$ を展開した多項式

$$ p(x) = x^9 – 18x^8 + 144x^7 – 672x^6 + 2016x^5 – 4032x^4 + 5376x^3 – 4608x^2 + 2304x – 512 $$

を、$x = 2$ の近くで評価します。$x \approx 2$ のとき各項は $2^9 = 512$ 前後の大きさを持ち、係数を掛けると $4032 \times 2^4 \approx 6.5\times10^4$ にもなります。それらが打ち消し合って、答えは $10^{-13}$ 以下の極小の値になるはずです。条件数は $\kappa \sim 10^{5}/10^{-13} = 10^{18}$。$\kappa u \gg 1$ で、答えは完全なノイズです。

実測すると、$x \in [1.9, 2.1]$ を $2001$ 点に区切って展開形 $p(x)$ を評価したところ、符号が515回変化しました。真の $(x-2)^9$ は $x=2$ でただ一度だけ符号を変えるのに、展開形では根が500個以上あるように見えるわけです。$|p(x)| < 10^{-11}$ を満たす $x$ の範囲は $[1.934, 2.068]$、幅 $0.133$。式(15)の意味するところそのままで、根を $\pm 0.067$ 以上の精度で特定することは原理的に不可能です。二分法でも、ニュートン法でも、反復回数をいくら増やしても改善しません。

実際にニュートン法を回すと、$x_0 = 2.3$ から始めて40回反復しても $x_{39} = 2.063$、誤差 $0.063$ にとどまります（重根なので収束が1次に落ちることも重なって、じりじりとしか進みません）。一方 $(x-2)^9$ を因数分解形のまま評価すれば、同じ反復でも桁落ちが起きないぶん素直に近づきます。

左のグラフでは、区間の端（$x = 1.9$ や $x = 2.1$）では展開形（青）と因数分解形（赤）がぴたりと重なっています。関数値が $10^{-9}$ 程度あって誤差に埋もれないからです。ところが黄色く塗った帯の中では、赤い曲線がほぼ $0$ に張りついているのに青が細かく震えており、$\pm 10^{-11}$ の値がすべてノイズであることがわかります。右はその帯の中心付近を拡大したもので、真の値（赤）が水平線にしか見えない一方、展開形（青）は $\pm 2 \times 10^{-11}$ の範囲で無秩序に暴れ、区間全体で符号が 515回も入れ替わっています。根が1つしかない関数が、計算機の上では「根が数百個あるように見える」わけです。

ここから引き出せる実務的な教訓は3つです。

1. 多項式は展開せず、因数分解形やHorner法など、打ち消しの少ない形で評価する。Horner法（$((\ldots(a_nx + a_{n-1})x + \ldots)x + a_0$）は乗算回数が減るだけでなく、丸め誤差の蓄積も抑えられます。
2. 停止判定に $|f(x)| < \varepsilon$ という絶対値の閾値を使わない。$f$ のスケールは問題ごとに違うので、$|f(x_k)| < \varepsilon_{\text{rel}}|f(x_0)|$ のような相対判定か、$|x_{k+1}-x_k| < \varepsilon_{\text{rel}}|x_k| + \varepsilon_{\text{abs}}$ という増分ベースの判定を使います。
3. 到達精度の見積もりを持つ。単根なら $\sqrt{u}$ ではなく $u$ 程度まで、重根度 $m$ の根なら $u^{1/m}$ 程度までしか到達できません。$(x-2)^9$ なら $u^{1/9} \approx 0.017$ で、実測の $0.06$ とオーダーが合います。

理論の道具立てはひと通り揃いました。ここからはPythonで、これまでの現象を自分の手で再現していきます。

Pythonでの実装

実装1: ビット表現を覗く

まず、float64の中身を実際に見るところから始めます。struct モジュールで64ビットの生パターンを取り出し、符号・指数・仮数に切り分けます。

import struct
from decimal import Decimal, getcontext

getcontext().prec = 60  # 多倍長で厳密値を表示するため

def float_bits(x: float):
    """float64 を (符号, 指数部11bit, 仮数部52bit) の文字列に分解する"""
    q = struct.unpack('<Q', struct.pack('<d', x))[0]  # 64bit整数として取り出す
    b = format(q, '064b')
    return b[0], b[1:12], b[12:]

def show(x: float, label: str = ""):
    s, e, m = float_bits(x)
    E = int(e, 2)
    print(f"{label:>10} {x!r:>22}")
    print(f"           符号={s}  指数部={e} (E={E}, e=E-1023={E-1023})")
    print(f"           仮数部={m}")
    print(f"           厳密値={Decimal(x)}")
    print()

for v, lab in [(1.0, "1.0"), (0.5, "0.5"), (0.1, "0.1"),
               (0.2, "0.2"), (0.3, "0.3"), (0.1 + 0.2, "0.1+0.2")]:
    show(v, lab)

出力を見ると、$0.1$ の仮数部が 1001100110011001...1001100110011010 と、$1001$ の繰り返しの末尾だけが $1010$ に切り上がっていることが確認できます。また $0.1$ と $0.2$ の仮数部は完全に同一で、指数部だけが $1019$ と $1020$ で1違う——$0.2 = 2 \times 0.1$ なのだから当然で、ビット列がそれをそのまま反映しています。そして 0.3 と 0.1+0.2 を比べると、仮数部の末尾が ...0011 と ...0100 で、ちょうど1ビットだけ違うことがわかります。前節で追跡した「同点をties-to-evenで上に丸めた」結果が、目に見える形で現れています。

実装2: 丸め誤差の上界と、隣り合う数の間隔

次に、$u = 2^{-53}$ という理論値が実際に守られているかを確認し、ULPの広がり方を可視化します。

import numpy as np
import matplotlib, matplotlib.pyplot as plt

for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
    if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
        plt.rcParams["font.family"] = cand
        break
plt.rcParams["axes.unicode_minus"] = False

u = 2.0 ** -53          # unit roundoff（丸め誤差の相対上界）
eps = 2.0 ** -52        # 機械イプシロン = ulp(1)
print(f"u   = {u:.6e}  (2^-53)")
print(f"eps = {eps:.6e}  (2^-52, np.finfo: {np.finfo(np.float64).eps:.6e})")

# ULP がスケールとともに 2 のべきで広がる様子
xs = np.logspace(-3, 17, 400)
ulps = np.array([np.spacing(x) for x in xs])

plt.figure(figsize=(9, 5))
plt.loglog(xs, ulps, lw=2, label="隣り合う float64 の間隔 (ULP)")
plt.loglog(xs, xs * eps, "--", lw=1.5, label=r"目安 $x \cdot 2^{-52}$")
plt.axhline(1.0, color="crimson", ls=":", lw=1.5, label="間隔=1（整数が飛び始める）")
plt.axvline(2.0 ** 53, color="gray", ls=":", lw=1.5, label=r"$2^{53} \approx 9.0\times10^{15}$")
plt.xlabel("数の大きさ $x$")
plt.ylabel("その付近での表現の粗さ (ULP)")
plt.title("float64 は「相対的に一定」の精度を持つ — ULP は $x$ に比例して広がる")
plt.legend(); plt.grid(alpha=0.3, which="both")
plt.tight_layout(); plt.show()

生成されるのは、前掲した ULP のグラフと同じものです。その傾きが $1$ の直線になっていることが、float64の設計思想そのものです。間隔は数の大きさに比例して広がる——だからこそ、相対誤差はどこでも $u = 2^{-53}$ 以下に保たれます。赤い水平線（間隔 $=1$）と交わる点が $2^{53} \approx 9.0\times10^{15}$ で、ここを超えると連続する整数を区別できなくなります。$10^{16}$ 付近では間隔が $2$、$10^{17}$ では $16$ です。「大きな整数IDをfloatで持ってはいけない」という定番の注意が、この図から直感的に理解できます。

実装3: 桁落ちの相対誤差は条件数で予言できる

理論の主役だった式(6) $\ \text{相対誤差} \le \kappa u\ $ を、実測で確かめます。$\sqrt{x+1}-\sqrt{x}$ を素直な式と有理化した式の両方で計算し、多倍長演算の真値と比べます。

import numpy as np, math
from decimal import Decimal, getcontext
import matplotlib.pyplot as plt

getcontext().prec = 60
u = 2.0 ** -53

def naive(x):    # 桁落ちする素直な式
    return math.sqrt(x + 1.0) - math.sqrt(x)

def stable(x):   # 有理化した式（引き算なし）
    return 1.0 / (math.sqrt(x + 1.0) + math.sqrt(x))

def exact(x):    # 多倍長での真値
    X = Decimal(x)
    return (X + 1).sqrt() - X.sqrt()

ks = np.arange(0, 17, 1.0)
xs = 10.0 ** ks
rel_naive, rel_stable, conds = [], [], []
for x in xs:
    ex = exact(x)
    rn = abs(Decimal(naive(x)) - ex) / ex
    rs = abs(Decimal(stable(x)) - ex) / ex
    kappa = (math.sqrt(x + 1.0) + math.sqrt(x)) / float(ex)   # 式(7)の条件数
    rel_naive.append(max(float(rn), 1e-18))
    rel_stable.append(max(float(rs), 1e-18))
    conds.append(kappa)

plt.figure(figsize=(9, 5.5))
plt.loglog(conds, rel_naive, "o-", lw=2, label=r"素直な式 $\sqrt{x+1}-\sqrt{x}$")
plt.loglog(conds, rel_stable, "s-", lw=2, label=r"有理化 $1/(\sqrt{x+1}+\sqrt{x})$")
plt.loglog(conds, np.array(conds) * u, "k--", lw=1.5, label=r"理論の上界 $\kappa \cdot u$")
plt.axhline(1.0, color="crimson", ls=":", label="相対誤差 1（有効数字ゼロ）")
plt.xlabel(r"引き算の条件数 $\kappa = (|a|+|b|)/|a-b|$")
plt.ylabel("相対誤差")
plt.title("桁落ちの誤差は条件数に比例して増える（有理化すれば増えない）")
plt.legend(); plt.grid(alpha=0.3, which="both")
plt.tight_layout(); plt.show()

生成されるのは前掲した条件数のグラフで、本記事でいちばん伝えたい一枚です。青い丸（素直な式）は、黒い破線（理論の上界 $\kappa u$）とぴったり平行に、傾き $1$ で増えていきます。条件数が10倍になれば誤差も10倍という式(6)の予言が、実測でそのまま成り立っています。$\kappa \approx 10^{16}$ で赤い線（相対誤差 $1$）に達し、答えは意味を失います。一方、緑の四角（有理化）は、条件数がいくら大きくなっても $10^{-16}$ 台の水平線のままです。同じ数学的値を計算しているのに、式の書き方だけで16桁の精度差がつく——これが数値的安定性という言葉の中身です。

実装4: カハン加算 vs 素直な加算の誤差蓄積

最後に、総和の誤差が $n$ とともにどう育つかを比較します。

import numpy as np, math
import matplotlib.pyplot as plt

def naive_sum(a):
    s = 0.0
    for v in a:
        s += v
    return s

def kahan_sum(a):
    s = 0.0
    c = 0.0                 # 前回捨てられた分（の符号反転）
    for v in a:
        y = v - c           # 捨てた分を今回の項に返す
        t = s + y           # ここで丸めが起きる
        c = (t - s) - y     # 実際に足された量 − 足したかった量 = 捨てられた分
        s = t
    return s

Ns = [10**3, 10**4, 10**5, 10**6]
err_naive, err_kahan, err_np = [], [], []
for N in Ns:
    a = np.full(N, 0.1)          # 0.1 は真値よりわずかに大きい → 誤差が一方向に溜まる
    S = math.fsum(a)             # 正しく丸められた真値
    err_naive.append(abs(naive_sum(a) - S))
    err_kahan.append(max(abs(kahan_sum(a) - S), 1e-18))
    err_np.append(max(abs(float(np.sum(a)) - S), 1e-18))

plt.figure(figsize=(9, 5.5))
plt.loglog(Ns, err_naive, "o-", lw=2, label="素直な総和（誤差 $\\propto n^2$）")
plt.loglog(Ns, err_np, "^-", lw=2, label="np.sum（pairwise, $O(u\\log n)$）")
plt.loglog(Ns, err_kahan, "s-", lw=2, label="カハン補正総和（$n$ に依らない）")
plt.xlabel("足し込む項数 $n$")
plt.ylabel("真値からの絶対誤差")
plt.title("総和アルゴリズムによる誤差蓄積の違い（$0.1$ を $n$ 個足す）")
plt.legend(); plt.grid(alpha=0.3, which="both")
plt.tight_layout(); plt.show()

for N, a, b, c in zip(Ns, err_naive, err_kahan, err_np):
    print(f"n={N:>8}  素直={a:.3e}  カハン={b:.3e}  np.sum={c:.3e}")

生成されるグラフは前掲した総和のグラフと同じもので、3本の線の傾きの違いが、そのままアルゴリズムの質の違いです。素直な総和（青丸）は傾きがほぼ $2$ で、$n$ が10倍になると誤差が100倍になります。すべての項が同符号なので誤差が打ち消し合わず、式(13)の上界がそのまま出ています。np.sum（緑の三角）は pairwise summation により、$n=10^6$ でも誤差 $2\times10^{-10}$ と素直な総和の約 $1/6500$ です。カハン加算（紫の四角）は $n$ によらず誤差 $0$ で、対数目盛では下限に張り付いています。「ループで足すな、np.sum を使え」という定番の助言は、速度の話であると同時に精度の話でもあるとわかります。

実装5: 分散の3つの計算法

理論編で見た「分散が負になる」現象を、自分の手で再現します。

import numpy as np
from decimal import Decimal, getcontext
getcontext().prec = 60

def var_expand(x):        # E[X^2] - (E[X])^2（桁落ちする）
    return float(np.mean(x * x) - np.mean(x) ** 2)

def var_twopass(x):       # 定義どおり2回なめる
    m = np.mean(x)
    return float(np.mean((x - m) ** 2))

def var_welford(x):       # 1パスで安定
    n, m, M2 = 0, 0.0, 0.0
    for v in x:
        n += 1
        d = v - m         # 更新前の平均とのずれ
        m += d / n        # 式(11)
        M2 += d * (v - m) # 式(12): (更新前のずれ)×(更新後のずれ)
    return M2 / n

def var_exact(x):         # 多倍長での真値
    n = len(x)
    mm = sum(Decimal(float(v)) for v in x) / n
    return sum((Decimal(float(v)) - mm) ** 2 for v in x) / n

rng = np.random.default_rng(0)
print(f"{'オフセット':>10} {'真値':>14} {'展開公式':>16} {'二パス':>16} {'Welford':>16}")
for shift in [0.0, 1e6, 1e8, 1e9]:
    x = rng.standard_normal(10_000) + shift
    ex = var_exact(x)
    print(f"{shift:>10.0e} {float(ex):>14.10f} {var_expand(x):>16.6f} "
          f"{var_twopass(x):>16.10f} {var_welford(x):>16.10f}")

実行すると、オフセット $10^8$ の行で展開公式が $-4.0$、オフセット $10^9$ で $128.0$ を返すのが確認できます。真値はどちらも $1.0$ 前後なので、符号すら合っていません。一方で二パス法とWelford法は、オフセットが $10^9$ になっても真値と小数点以下7桁以上一致します。「数学的に正しい式」と「計算機で正しい式」は別物だということを、これ以上ないほど明快に示す実験です。ストリーミング処理で分散を追いたい場合は、迷わずWelford法を使ってください。

実装6: float32との比較と、等値比較の代わり

最後に、精度を落とした np.float32 での挙動と、比較のやり方を確認します。

import numpy as np, math, struct

f32 = np.float32
print("float32 の eps =", np.finfo(f32).eps, " (2^-23 =", 2.0**-23, ")")
print("float64 の eps =", np.finfo(np.float64).eps, " (2^-52 =", 2.0**-52, ")")

# float32 では 0.1+0.2 == 0.3 が成り立ってしまう
print("float32: 0.1+0.2 =", float(f32(0.1) + f32(0.2)), " 0.3 =", float(f32(0.3)))
print("float32: 0.1+0.2 == 0.3 ->", bool(f32(0.1) + f32(0.2) == f32(0.3)))
print("float64: 0.1+0.2 == 0.3 ->", 0.1 + 0.2 == 0.3)

# 等値比較の代わりに許容誤差つき比較を使う
print("np.isclose(0.1+0.2, 0.3)   ->", np.isclose(0.1 + 0.2, 0.3))
print("math.isclose(0.1+0.2, 0.3) ->", math.isclose(0.1 + 0.2, 0.3))

# 2つの isclose は「ゼロとの比較」で挙動が違う
print("np.isclose(0.0, 1e-9)      ->", np.isclose(0.0, 1e-9))     # True（atol が効く）
print("math.isclose(0.0, 1e-9)    ->", math.isclose(0.0, 1e-9))   # False（相対のみ）

まず驚くのは、float32では 0.1+0.2 == 0.3 が True になることです。精度が低いほうが「正しい」答えを返すという逆説的な結果ですが、理由は単純で、float32の仮数は23ビットしかなく、丸めの粒が粗いために $0.1+0.2$ の誤差が $0.3$ の粒の中に埋もれてしまうからです。これは「float32のほうが安全」という意味ではまったくなく、等値比較の結果は精度に依存する当てにならない指標だという教訓です。float32の $\varepsilon = 2^{-23} \approx 1.19\times10^{-7}$ は有効数字約7桁で、機械学習の推論などでは十分でも、長時間積分や桁落ちのある計算では危険です。

そして np.isclose と math.isclose の違いに注目してください。ゼロとの比較で結果が分かれます。理由は判定式が違うからです。

$$ \begin{align} \texttt{np.isclose(a,b)} &: \quad |a – b| \le \texttt{atol} + \texttt{rtol}\cdot|b| \\ \texttt{math.isclose(a,b)} &: \quad |a – b| \le \max\left(\texttt{rel\_tol}\cdot\max(|a|,|b|),\ \texttt{abs\_tol}\right) \end{align} $$

np.isclose は既定で rtol=1e-05, atol=1e-08 を持ち、絶対許容誤差 atol があるのでゼロとの比較が機能します。ただし $|b|$ だけを使う非対称な判定なので、np.isclose(a, b) と np.isclose(b, a) が一致しない場合があります。一方 math.isclose は既定で abs_tol=0.0 なので、片方がちょうど $0$ だと相対判定が破綻して常に False を返します。ゼロと比べる可能性があるなら abs_tol を明示的に指定しなければなりません。

使い分けの指針は次の通りです。

• スケールが既知でゼロを含まない量（長さ、質量、確率密度など）→ 相対許容誤差だけで十分。math.isclose(a, b, rel_tol=1e-12)
• ゼロになりうる量（残差、差分、誤差そのもの）→ 絶対許容誤差が必須。math.isclose(a, b, rel_tol=1e-9, abs_tol=1e-12) や np.isclose の既定
• 配列全体を一括判定したい → np.allclose(A, B)、テストなら np.testing.assert_allclose(A, B, rtol=..., atol=...)
• 金額など10進で厳密に扱いたい → そもそも浮動小数点を使わず decimal.Decimal や整数（最小単位で保持）を使う

「浮動小数点を == で比べない」は有名なルールですが、大事なのはその先——許容誤差をどう決めるかです。誤差の見積もりは本記事で見てきた通り、$u$ と条件数と演算回数から立てられます。たとえば $n$ 回の加算を経た量なら $nu$ 程度、条件数 $\kappa$ の引き算を含むなら $\kappa u$ 程度を目安にする。根拠のある許容誤差を置けるかどうかが、数値計算コードの品質を分けます。

実務でのチェックリスト

ここまでの内容を、コードを書くときにその場で使える形にまとめます。

1. 近い数どうしの引き算を探す。式の中に $a – b$ があって $a \approx b$ になりうるなら、有理化・三角関数の加法定理・テイラー展開などで書き換えられないか検討します。$\sqrt{x+1}-\sqrt{x}$、$1-\cos x$、$e^x-1$、$\log(1+x)$ は定番の書き換え対象です。
2. 標準ライブラリの安定版関数を使う。math.expm1、math.log1p、math.hypot（$\sqrt{x^2+y^2}$ をオーバーフローなしで）、scipy.special.logsumexp、np.logaddexp は、すべて桁落ち・オーバーフロー対策済みです。自作しないでください。
3. 統計量は np.var / np.std またはWelford法。$E[X^2]-(E[X])^2$ を手書きしない。共分散行列も同様で、$E[XY]-E[X]E[Y]$ ではなく中心化してから計算します。
4. 総和は np.sum / math.fsum。Pythonの for ループで足し込むのは、遅いだけでなく精度も悪い。
5. データを中心化・スケーリングしてから計算する。オフセットの大きいデータ（時刻、座標、ADCカウント）は、平均や基準点を引いてから処理すると条件数が劇的に下がります。回帰でも最適化でも、前処理の中心化には数値的な意味があります。
6. == を使わず許容誤差つき比較を使い、その許容誤差に根拠を持つ。テストなら np.testing.assert_allclose。
7. 整数として扱うべきものはfloatに入れない。$2^{53}$ を超えるIDや通し番号は、Pythonの int（多倍長）か文字列で持ちます。JSONを経由すると整数がfloatに化けることがあるので要注意です。
8. 金額は Decimal か「最小単位の整数」で。円なら円単位の整数、ドルならセント単位の整数。浮動小数点で会計処理をしてはいけません。
9. 反復法の停止判定は相対 + 絶対の組み合わせ。$|x_{k+1}-x_k| < \varepsilon_{\text{rel}}|x_k| + \varepsilon_{\text{abs}}$ の形にし、無限ループを避けるため最大反復回数も必ず設定します。
10. 疑ったら多倍長で答え合わせ。decimal.Decimal（getcontext().prec で桁数指定）や fractions.Fraction（厳密な有理数）、mpmath で真値を出し、float64の結果と突き合わせます。本記事の実験もすべてこの方法で検証しています。

まとめ

本記事では、IEEE 754浮動小数点数の内部構造から出発して、数値誤差がどこで生まれ、どこで増幅されるのかを追いました。

• $0.1$ が表せない理由は、$1/10$ の分母に $5$ が残るから。2進法で有限桁になるのは分母が2のべき乗の分数だけで、$0.1$ は $0.0\overline{0011}_{(2)}$ という循環小数になる
• float64のビット配置は符号1・指数11・仮数52ビット。暗黙の1により実効53ビットの精度を持ち、値は $(-1)^s(1+M/2^{52})2^{E-1023}$
• 丸め誤差の相対上界は $u = 2^{-53} \approx 1.11\times10^{-16}$。区間ごとの間隔 $2^{e-52}$ の半分を $|x|\ge 2^e$ で割ると $e$ が消えることから導かれる。機械イプシロン $\varepsilon = 2^{-52}$ はその2倍（$\mathrm{ulp}(1)$）で、混同しやすいので注意
• $0.1+0.2$ の厳密和は $\mathrm{fl}(0.3)$ と次の数のちょうど中間に落ち、ties-to-evenで上側に丸められる。1ビット、$5.55\times10^{-17}$ のずれ
• 桁落ちの相対誤差は条件数 $\kappa = (|a|+|b|)/|a-b|$ 倍に増幅される。$\sqrt{x+1}-\sqrt{x}$ も二次方程式の解の公式も、有理化すれば引き算が消えて機械精度まで回復する
• 分散の展開公式 $E[X^2]-(E[X])^2$ は $\kappa \approx 2\bar{x}^2/V[X]$ の条件数を持ち、オフセット $10^8$ のデータで $-4.0$ という負の分散を返す。二パス法かWelfordの逐次更新 $M_{2,n} = M_{2,n-1} + (x_n-m_{n-1})(x_n-m_n)$ を使う
• 総和の誤差は素直に足すと $O(nu)$ で育つ。カハンの補正総和は捨てられた下位ビットを $c = (t-s)-y$ で回収して次に返すことで、誤差を $n$ に依存しない $O(u)$ に抑える。np.sum の pairwise summation は $O(u\log n)$、math.fsum は常に正しく丸められた答えを返す
• 等値比較は精度に依存して結果が変わる（float32では 0.1+0.2 == 0.3 が True）。np.isclose / math.isclose を使い、ゼロと比較する可能性があるなら絶対許容誤差を必ず指定する

浮動小数点の誤差は「避けられない不運」ではありません。$u$ と条件数という2つの数字さえ手元にあれば、どこで誤差が増幅するかは計算前に予測でき、式を書き換えることで回避できます。数値計算とは、同じ数学的な値に至る無数の式の中から、条件数の小さい経路を選ぶ作業なのだと考えると、見通しがよくなるはずです。

次のステップとして、以下の記事も参考にしてください。

• 条件数と数値的安定性の理論 — 本記事の条件数を、行列と連立一次方程式に拡張します
• 【Python】数値計算で微分可能な関数の微分値を求める — 刻み幅 $h$ の最適値が $\sqrt{u} \approx 1.5\times10^{-8}$ で決まる理由を、本記事の丸め誤差モデルから導けます
• 数値線形代数の基礎 — LU分解の部分ピボット選択や、QR分解が数値的に好まれる理由につながります

参考文献

• IEEE Std 754-2019, IEEE Standard for Floating-Point Arithmetic
• David Goldberg, “What Every Computer Scientist Should Know About Floating-Point Arithmetic,” ACM Computing Surveys, 23(1), 1991
• Nicholas J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002
• B. P. Welford, “Note on a Method for Calculating Corrected Sums of Squares and Products,” Technometrics, 4(3), 1962
• W. Kahan, “Further remarks on reducing truncation errors,” Communications of the ACM, 8(1), 1965