# 浮動小数点数と丸め誤差：IEEE 754・標準モデル・和と Horner 法の誤差解析 > IEEE 754 の 2 進形式と正しい丸めの規約から四則の標準モデルを導き、Sterbenz の補題と桁落ち、逐次和・対ごとの和の誤差評価、Horner 法の後退誤差を証明する。binary32 の和と 1−cos x の数値例で評価を確かめる。 > https://rikai.mugen-giken.com/mathematics/numerical-analysis/floating-point-arithmetic ## 0. この記事の要点 計算機は無限に長い小数を持てないので、実数を有限の桁で丸めて記憶し、計算のたびにまた丸めます。この記事は、その丸めが何度も重なると答えがどれだけずれるのかを、不等式で見積もる方法を学びます。読み終えると、足し算や多項式の計算で誤差がどの程度になるかを式で説明でき、`0.1 + 0.2` が `0.3` にならない理由も数字で追えるようになります。 - 計算機の実数は IEEE 754 の**浮動小数点数**（有限桁の 2 進数で実数を近似して表す方式）です。binary64（倍精度）は符号 1 ビット・指数 11 ビット・仮数 52 ビットで、精度 $p=53$ 桁の 2 進数を表します。 - 四則と平方根は**正しく丸められる**（厳密な結果を最も近い浮動小数点数に丸めたものが返る）と規格が定めています。ここから、オーバーフローとアンダーフローがなければ $\mathrm{fl}(x\circ y)=(x\circ y)(1+\delta)$, $|\delta|\le u=2^{-p}$ という**標準モデル**が従います。 - **Sterbenz の補題**: $y/2\le x\le 2y$ なら $x-y$ は誤差なく計算されます。桁落ちで失われる精度は、引き算そのものではなく、引く前の値に含まれていた誤差から来ます。 - $n$ 個の逐次和の誤差は $\gamma_{n-1}\sum|x_i|$（$\gamma_k=ku/(1-ku)$）以下、隣どうしを足していく**対ごとの和**は $\gamma_{\lceil\log_2n\rceil}\sum|x_i|$ 以下です。 - **Horner 法**で計算した多項式の値は、各係数を相対 $\gamma_{2n}$ 以内で動かした多項式の正確な値です（後退誤差）。この「計算結果は少し違う入力に対する正確な答え」という見方が、数値解析と数値線形代数の誤差解析の出発点になります。 ## 1. 動機 Python で `0.1 + 0.2 == 0.3` を評価すると `False` になります。数値計算を始めた人がまず出会う現象ですが、これは計算機の「誤り」ではありません。十進の $0.1$ は 2 進では有限桁で書けないので、読み込んだ時点で最も近い 2 進数に置き換わり、足し算の結果もまた最も近い 2 進数に置き換わります。どの段階でどれだけずれるかは、規格によって**正確に**決まっています。 数値解析の仕事は、この「置き換え」が積み重なったときに答えがどれだけずれるかを、定理として評価することです。ずれには 2 種類あります。 | 種類 | 原因 | 例 | |---|---|---| | 丸め誤差 | 実数を有限桁の浮動小数点数で表すこと | 和・内積・多項式の値の誤差 | | 打ち切り誤差 | 無限の操作（極限・級数・積分）を有限で止めること | 差分商、数値積分、反復の停止 | この記事は前者を扱います。丸めの基本的な誤差モデルと誤差因子の積の補題は、計算機援用証明のテーマの[浮動小数と区間演算](/mathematics/computer-assisted-proofs/floating-point-interval-arithmetic)の章で証明されているので、ここではその主張を引用して使い、IEEE 754 の形式と、和・桁落ち・多項式の誤差解析を積み上げます。 ## 2. IEEE 754 の 2 進形式 ### 2.1. 浮動小数点数の集合 精度 $p$ と指数の範囲 $e_{\min}\lt e_{\max}$ を固定した浮動小数点数の集合 $\mathbb{F}$ は、既存の定義（）を使います。要点を再掲します。 - **正規化数**は $\pm M\cdot2^{e-p+1}$（$M$ は整数で $2^{p-1}\le M\lt2^p$、$e_{\min}\le e\le e_{\max}$）、**非正規化数**は $\pm M\cdot2^{e_{\min}-p+1}$（$0\lt M\lt2^{p-1}$）で、これらと $0$ を合わせたものが $\mathbb{F}$ です。 - 最小の正規化数は $x_{\min}=2^{e_{\min}}$、最大の数は $x_{\max}=(2-2^{1-p})2^{e_{\max}}$、**単位丸め**は $u=2^{-p}$ です。 - $[2^e,2^{e+1}]$ の中の浮動小数点数は間隔 $2u\cdot2^e$ で等間隔に並びます（）。間隔は $2$ の冪をまたぐたびに 2 倍になり、相対的な間隔は $u$ から $2u$ の間にあります。
1/2 1 2 間隔 2⁻⁴ 間隔 2⁻³ 間隔 2⁻²
### 2.2. ビット配置 IEEE 754 は、この集合を固定長のビット列で表す方法を定めています。次の定義は記号が多いですが、考え方は「符号 1 ビット、指数 $E$、仮数の小数部 $f$ の 3 つの整数から、$\pm2^{(\text{指数})}\times1.(\text{小数部})$ の形の数を作る」というだけです。たとえば binary64 では、$E=1023$, $f=0$ が $1$、$E=1024$, $f=0$ が $2$、$E=1023$, $f=2^{51}$ が $1+2^{51}\cdot2^{-52}=1.5$ を表します。
s E f 1 ビット 11 ビット 52 ビット 符号 指数部 仮数部
正の整数 $w$（指数部の幅）と $t$（仮数部の幅）を固定し、$\mathrm{bias}=2^{w-1}-1$ とする。$1+w+t$ ビットの列を、符号ビット $s\in\{0,1\}$、指数部 $E\in\{0,1,\ldots,2^w-1\}$（$w$ ビットの符号なし整数）、仮数部 $f\in\{0,1,\ldots,2^t-1\}$（$t$ ビットの符号なし整数）に分け、次の値を表すものとする。 1. $1\le E\le2^w-2$ のとき **正規化数** $(-1)^s\,2^{E-\mathrm{bias}}\,(1+f\,2^{-t})$。 2. $E=0$, $f\ne0$ のとき **非正規化数** $(-1)^s\,2^{1-\mathrm{bias}}\,(f\,2^{-t})$。$E=0$, $f=0$ のとき $(-1)^s\,0$（$+0$ と $-0$）。 3. $E=2^w-1$, $f=0$ のとき $(-1)^s\infty$。$E=2^w-1$, $f\ne0$ のとき **NaN**（数でない値）。 binary32 は $(w,t)=(8,23)$、binary64 は $(w,t)=(11,52)$ である。 1 と 2 の値の全体は、$p=t+1$, $e_{\min}=1-\mathrm{bias}$, $e_{\max}=\mathrm{bias}$ とした $\mathbb{F}$ に一致します。実際、1 の値は $M=2^t+f$（$2^{p-1}\le M\lt2^p$）、$e=E-\mathrm{bias}$ として $M\cdot2^{e-p+1}$ と書け、2 の値は $M=f$（$0\lt M\lt2^{p-1}$）として $M\cdot2^{e_{\min}-p+1}$ と書けます。正規化数では先頭の 1 を格納しない（**隠れビット**）ので、$t$ ビットの仮数部で $p=t+1$ 桁の精度が得られます。 | 形式 | 全体 | $w$ | $t$ | $p$ | bias | $e_{\min}$ | $e_{\max}$ | $u=2^{-p}$ | |---|---|---|---|---|---|---|---|---| | binary32（単精度） | 32 | 8 | 23 | 24 | 127 | $-126$ | 127 | $\approx5.96\times10^{-8}$ | | binary64（倍精度） | 64 | 11 | 52 | 53 | 1023 | $-1022$ | 1023 | $\approx1.11\times10^{-16}$ | 非正規化数があるおかげで、$0$ と $x_{\min}$ の間も等間隔 $2^{e_{\min}-p+1}$ で埋まります。$x\ne y$ の浮動小数点数について $x-y$ が $0$ に丸められることはありません（**漸進アンダーフロー**）。非正規化数がなければ、$x_{\min}$ と $1.5\,x_{\min}$ の差は表せず $0$ になり、`if x != y: z = 1/(x - y)` のような計算が 0 除算を起こします。 $\pm0$ は比較では等しく、$1/(+0)=+\infty$, $1/(-0)=-\infty$ のように符号が意味をもつ演算があります。NaN は $0/0$, $\infty-\infty$, $\sqrt{-1}$ などの結果で、自分自身とも等しくない（`x != x` が真）ことで判定できます。 ### 2.3. 丸めの方向と正しい丸め $|x|\le x_{\max}$ の実数 $x$ について、次の 4 つの丸めを考える。 - **最近接偶数丸め** $\mathrm{fl}(x)$: $x$ に最も近い $\mathbb{F}$ の元。2 つあれば $M$ が偶数の方（ の $\mathrm{fl}$）。 - **下向き丸め** $\nabla(x)$、**上向き丸め** $\Delta(x)$（同じ定義の $\nabla$, $\Delta$）、**0 方向への丸め**（$x\ge0$ なら $\nabla(x)$、$x\lt0$ なら $\Delta(x)$）。 演算 $\circ\in\{+,-,\times,\div\}$ と $x,y\in\mathbb{F}$ について、計算機が返す値が厳密な実数 $x\circ y$ をいずれかの丸め（既定は最近接偶数丸め）で丸めたものであるとき、その演算は**正しく丸められる**という。平方根 $\sqrt{x}$（$x\ge0$）と融合積和 $xy+z$ についても同様に定める。 IEEE 754 は、四則・平方根・融合積和（fma）が正しく丸められることを**要求**します。以下、計算機の演算の結果を $\mathrm{fl}(x\circ y)$ と書きます。これは「厳密な $x\circ y$ を計算してから丸めた値」という意味で、どのように計算するかには立ち入りません。一方、`exp`・`sin` などの初等関数については、規格は正しい丸めを推奨するにとどめているので、以下の評価をそのまま適用できません。 $1$ と、$1$ より大きい最小の浮動小数点数との差 $\varepsilon_M=2^{1-p}$ を**機械イプシロン**という。その半分 $u=2^{-p}$ を**単位丸め**という。 $\varepsilon_M$ は $1$ の右側の間隔で、$u$ は最近接丸めの相対誤差の上界です（）。$1$ の左側の間隔は $u$ です（$[1/2,1)$ の間隔は $2u\cdot2^{-1}$）。困ったことに、文献ごとに同じ記号が違う量を指します。 | 文献・環境 | $2^{1-p}$ の呼び名 | $2^{-p}$ の呼び名 | |---|---|---| | Higham | machine epsilon $\epsilon_M$ | unit roundoff $u$ | | Golub–Van Loan | — | unit roundoff $\mathbf{u}$ | | Trefethen–Bau | — | $\epsilon_{\text{machine}}$ | | C の `DBL_EPSILON`、NumPy の `finfo.eps`、MATLAB の `eps` | eps | — | 論文やライブラリの誤差評価を読むときは、どちらの量かを確かめてください。この記事では $u=2^{-p}$ を使います。 `np.finfo` と `np.nextafter` で、binary32 と binary64 の精度・指数の範囲・1 の両側の間隔・最小の正規化数と非正規化数・最大値を出します。 ```python # Python 3.12.10 / NumPy 2.5.2（x86-64） import numpy as np for t in (np.float32, np.float64): fi = np.finfo(t) one = t(1) up = np.nextafter(one, t(2)) - one # 1 の右の間隔 dn = one - np.nextafter(one, t(0)) # 1 の左の間隔 print(fi.dtype, fi.nmant + 1, fi.minexp, fi.maxexp - 1) print(f" eps={fi.eps:.3e} up={up:.3e} dn={dn:.3e}") print(f" tiny={fi.tiny:.3e} sub={fi.smallest_subnormal:.3e}" f" max={fi.max:.3e}") ``` ```text float32 24 -126 127 eps=1.192e-07 up=1.192e-07 dn=5.960e-08 tiny=1.175e-38 sub=1.401e-45 max=3.403e+38 float64 53 -1022 1023 eps=2.220e-16 up=2.220e-16 dn=1.110e-16 tiny=2.225e-308 sub=4.941e-324 max=1.798e+308 ``` **検算.** (a) 既知値: $p=24,53$、$e_{\min}=-126,-1022$、$e_{\max}=127,1023$ は表と一致し、`eps` $=2^{-23}\approx1.192\times10^{-7}$, $2^{-52}\approx2.220\times10^{-16}$、`tiny` $=2^{-126}$, $2^{-1022}$、`sub` $=2^{-149}$, $2^{-1074}$ です。(b) 別の求め方: `nextafter` で測った右の間隔 `up` は `finfo.eps` と一致し、左の間隔 `dn` はその半分の $u$ です。(c) 次数: 左右の間隔の比が 2 であることは、間隔が $2$ の冪をまたぐと 2 倍になること（）のとおりです。 ## 3. 四則の標準モデル 正しい丸めの規約と、丸めの誤差モデル（）を組み合わせると、すべての誤差解析の出発点になる次の形が得られます。 $x,y\in\mathbb{F}$、$\circ\in\{+,-,\times,\div\}$（$\div$ では $y\ne0$）とし、厳密な値 $z=x\circ y$ が $z=0$ または $x_{\min}\le|z|\le x_{\max}$ を満たす（オーバーフローもアンダーフローも起きない）とする。演算が最近接偶数丸めで正しく丸められるとき、 $$ \mathrm{fl}(x\circ y)=(x\circ y)(1+\delta),\qquad\mathrm{fl}(x\circ y)=\frac{x\circ y}{1+\delta'},\qquad|\delta|,|\delta'|\le u $$ を満たす $\delta,\delta'$ がある。$x\circ y\in\mathbb{F}$ ならば $\delta=\delta'=0$ にとれる。$x\ge0$ の平方根 $\mathrm{fl}(\sqrt{x})$ についても同じことが成り立つ。 正しい丸めの規約（）により、計算結果は $\mathrm{fl}(z)$ です。$z\in\mathbb{F}$（$z=0$ を含む）なら丸めの定義から $\mathrm{fl}(z)=z$ で、$\delta=\delta'=0$ とします。そうでなければ $x_{\min}\le|z|\le x_{\max}$ なので、丸めの誤差モデル（ の 1）により、$\mathrm{fl}(z)=z(1+\delta)$, $|\delta|\le u$ であり、さらに $|\mathrm{fl}(z)-z|\le u\,|\mathrm{fl}(z)|$ です。後者から $\mathrm{fl}(z)\ne0$（$\mathrm{fl}(z)=0$ なら $z=0$ となり仮定に反する）なので、$\delta'=(z-\mathrm{fl}(z))/\mathrm{fl}(z)$ とおくと $|\delta'|\le u$ かつ $z=\mathrm{fl}(z)(1+\delta')$、すなわち $\mathrm{fl}(z)=z/(1+\delta')$ です。平方根も $z=\sqrt{x}$ として同じです。 2 つ目の形 $z/(1+\delta')$ は、誤差因子を分母に置きたいとき（割り算や、計算値の側から評価したいとき）に便利です。以下、この 2 つの形を単に「標準モデル」と呼び、途中でオーバーフローとアンダーフローが起きないことを常に仮定します。 丸めが何回も重なったときは、誤差因子をまとめる記号 $\gamma_n=nu/(1-nu)$（$nu\lt1$。）と、誤差因子の積の補題を使います。 > **誤差因子の積の補題**（）: $|\delta_i|\le u$, $\rho_i=\pm1$, $nu\lt1$ なら $\prod_{i=1}^n(1+\delta_i)^{\rho_i}=1+\theta_n$, $|\theta_n|\le\gamma_n$。 言い換えると、相対誤差 $u$ 以下の丸めが $n$ 回重なっても、全体の相対誤差は $\gamma_n$（$nu$ よりわずかに大きい程度）以下に収まる、ということです。以降の誤差解析は、数えた丸めの回数 $n$ を $\gamma_n$ に置き換える作業の繰り返しです。 この補題と組み合わせて使う $\gamma_n$ の計算規則を、補題として証明しておきます。数値線形代数の後退誤差解析で繰り返し使います。 $m,n\ge1$ を整数とする。 1. $m\le n$, $nu\lt1$ なら $\gamma_m\le\gamma_n$。また整数 $k\ge1$ について $knu\lt1$ なら $k\gamma_n\le\gamma_{kn}$。 2. $(m+n)u\lt1$ なら $\gamma_m+\gamma_n+\gamma_m\gamma_n\le\gamma_{m+n}$。とくに $\gamma_m+\gamma_n\le\gamma_{m+n}$。 3. $(m+n)u\lt1$ とする。$|\theta_m|\le\gamma_m$, $|\theta_n|\le\gamma_n$ なら $(1+\theta_m)(1+\theta_n)=1+\theta$ と書いたとき $|\theta|\le\gamma_{m+n}$。 1. $\gamma_k=1/(1-ku)-1$ は $k$ について増加です。また $k\gamma_n=knu/(1-nu)\le knu/(1-knu)=\gamma_{kn}$ です（分母 $1-nu\ge1-knu\gt0$）。 2. $1+\gamma_k=1/(1-ku)$ と、$(1-mu)(1-nu)=1-(m+n)u+mnu^2\ge1-(m+n)u\gt0$ から $$ (1+\gamma_m)(1+\gamma_n)=\frac{1}{(1-mu)(1-nu)}\le\frac{1}{1-(m+n)u}=1+\gamma_{m+n} $$ で、両辺から $1$ を引けば主張です。$\gamma_m\gamma_n\ge0$ なので後半も従います。 3. $\theta=\theta_m+\theta_n+\theta_m\theta_n$ なので、2 から $|\theta|\le\gamma_m+\gamma_n+\gamma_m\gamma_n\le\gamma_{m+n}$ です。 binary64 で $\mathrm{fl}(0.1)=M\cdot2^{-55}$, $\mathrm{fl}(0.2)=M\cdot2^{-54}$ です（$M=3602879701896397$。十進の $0.1$ の 2 進展開 $0.0\overline{0011}$ を 53 桁で丸めたもの）。厳密な和は $3M\cdot2^{-55}=10808639105689191\cdot2^{-55}$ で、$[1/4,1/2)$ の間隔 $2^{-54}$ で測ると $5404319552844595.5$ 目盛り、ちょうど 2 つの浮動小数点数の中点です。最近接偶数丸めは偶数の $5404319552844596\cdot2^{-54}$ を選びます。一方、十進の $0.3$ を丸めた $\mathrm{fl}(0.3)$ は $5404319552844595\cdot2^{-54}$ なので、両者は 1 目盛りずれます。有理数で計算すると $\mathrm{fl}(0.3)-3/10\approx-1.11\times10^{-17}$ で、これは 1 回の丸めの範囲 $u\cdot0.3\approx3.3\times10^{-17}$ に入ります。一方 $\mathrm{fl}(0.1+0.2)-3/10\approx4.44\times10^{-17}$ はこの範囲を超えています。$0.1$ と $0.2$ を読み込むときの 2 回の丸めの誤差（和は $\approx1.67\times10^{-17}$）が加わっているからです。足し算そのものの誤差は、厳密な和 $3M\cdot2^{-55}$ との差 $2^{-55}\approx2.8\times10^{-17}$ で、 の範囲 $u\,|3M\cdot2^{-55}|\approx3.3\times10^{-17}$ に入っています。 ## 4. 引き算：Sterbenz の補題と桁落ち ### 4.1. 誤差なく計算される引き算 $x,y\in\mathbb{F}$ が $y/2\le x\le2y$ を満たすならば $x-y\in\mathbb{F}$ であり、したがって $\mathrm{fl}(x-y)=x-y$（引き算は誤差なく計算される）。 証明の方針は、$x-y$ が $y$ の最下位桁の位置を単位とする整数倍で、しかも絶対値が $y$ 以下だから、$p$ 桁に収まる、というものです。非正規化数の場合分けが入るために長く見えますが、中心の考えはこれだけです。 仮定から $x,y\ge0$ です。$x=0$ なら $y\le2x=0$ で $y=0$、$x-y=0\in\mathbb{F}$ なので、$x,y\gt0$ とします。$\mathbb{F}$ は符号について対称なので、$x\ge y$ の場合に示せば十分です（$x\lt y$ なら $y-x\in\mathbb{F}$ を示して符号を変えます。仮定は $x$ と $y$ について対称です）。 正の浮動小数点数 $z$ に対し、$e_z=\max(\lfloor\log_2z\rfloor,e_{\min})$ とおきます。 から、$z$ は $z=M_z\cdot2^{e_z-p+1}$（$M_z$ は整数で $0\lt M_z\lt2^p$）と書けます（正規化数なら $e_z$ はその指数、非正規化数なら $e_z=e_{\min}$）。$e_z$ は $z$ について単調非減少なので、$x\ge y$ から $e_x\ge e_y$ です。よって $x=M_x2^{e_x-e_y}\cdot2^{e_y-p+1}$ は $2^{e_y-p+1}$ の整数倍で、 $$ x-y=k\cdot2^{e_y-p+1}\qquad(k\ \text{は整数}) $$ です。$0\le x-y\le2y-y=y$ から $0\le k\le M_y\lt2^p$ です。 あとは、$0\le k\lt2^p$ の整数 $k$ と $e_{\min}\le e\le e_{\max}$ について $k\cdot2^{e-p+1}\in\mathbb{F}$ を示せば十分です。$k=0$ なら $0\in\mathbb{F}$ です。$k\ge2^{p-1}$ なら指数 $e$ の正規化数です。$k\lt2^{p-1}$ のときは、$k$ を 2 倍して $e$ を 1 減らす操作を、$k\ge2^{p-1}$ になるか $e=e_{\min}$ になるまで繰り返します（値は変わりません）。前者で止まれば正規化数、後者で止まれば $0\lt k\lt2^{p-1}$ で非正規化数です。 この補題は「近い数どうしの引き算は**それ自体は**正確」と言っています。たとえば $x\in[1,2]$ の浮動小数点数について $\mathrm{fl}(x-1)=x-1$ です。 ### 4.2. 桁落ち それでも近い数の引き算が危険なのは、引く前の値がすでに誤差を含んでいるからです。 $a,b$ を $a\ne b$ の実数とし、$\hat a=a(1+\alpha)$, $\hat b=b(1+\beta)$, $|\alpha|,|\beta|\le\varepsilon$ とする。このとき $$ \frac{|(\hat a-\hat b)-(a-b)|}{|a-b|}\le\varepsilon\,\frac{|a|+|b|}{|a-b|} $$ であり、$a,b\gt0$ のとき等号は $\alpha=\varepsilon$, $\beta=-\varepsilon$ でとられる。 $(\hat a-\hat b)-(a-b)=a\alpha-b\beta$ なので、三角不等式から $|a\alpha-b\beta|\le\varepsilon(|a|+|b|)$ です。$a,b\gt0$, $\alpha=\varepsilon$, $\beta=-\varepsilon$ なら $a\alpha-b\beta=\varepsilon(a+b)$ で等号です。 $a\approx b$ なら倍率 $(|a|+|b|)/|a-b|$ は大きく、入力の小さな相対誤差 $\varepsilon$ が大きな相対誤差に拡大されます。これを**桁落ち**（cancellation）といいます。 と合わせると、桁落ちの誤差は**引き算の前に**入ったもので、引き算はそれを露わにするだけだとわかります。対策は、引き算の前に誤差の入らない式に書き換えることです。 $x\in\mathbb{F}$, $0\lt x\lt1$ とし、$y=\mathrm{fl}(\mathrm{fl}(1+x)-1)$ を計算します。標準モデルから $\mathrm{fl}(1+x)=(1+x)(1+\delta)$, $|\delta|\le u$ で、$\mathrm{fl}(1+x)\in[1,2]$ なので により引き算は厳密です。したがって $y=(1+x)(1+\delta)-1=x+(1+x)\delta$ で、相対誤差は $$ \frac{|y-x|}{x}=\frac{(1+x)|\delta|}{x}\le\frac{(1+x)u}{x} $$ です。$x=10^{-10}$ なら $u/x\approx10^{-6}$ で、16 桁のうち 10 桁が失われます。これは の倍率 $(|1+x|+1)/x$ と同じ大きさです。 $x$ が小さいとき $1-\cos x\approx x^2/2$ は $1$ よりずっと小さく、$\cos x$ の計算値に含まれる $u$ 程度の誤差が、相対誤差 $u/(x^2/2)$ に拡大されます。$1-\cos x=2\sin^2(x/2)$ と書き換えると、引き算がなくなり、どの段も相対誤差 $u$ 程度で済みます。$(1-\cos x)/x^2$ を両方の式で計算し、mpmath の 50 桁の値と比べます。 ```python # Python 3.12.10 / NumPy 2.5.2 / mpmath 1.3.0（x86-64） import numpy as np import mpmath as mp mp.mp.dps = 50 for k in range(1, 9): x = 10.0 ** -k exact = (1 - mp.cos(mp.mpf(x))) / mp.mpf(x) ** 2 # x は binary64 の値 naive = (1 - np.cos(x)) / x**2 stable = 2 * np.sin(x / 2) ** 2 / x**2 e1 = abs(naive - exact) / exact e2 = abs(stable - exact) / exact print(f"1e-{k} {float(e1):.2e} {float(e2):.2e}") ``` ```text 1e-1 1.11e-14 1.24e-16 1e-2 2.88e-13 1.14e-16 1e-3 1.57e-11 1.39e-17 1e-4 5.24e-09 1.53e-16 1e-5 8.27e-08 1.12e-16 1e-6 8.89e-05 4.44e-17 1e-7 7.99e-04 5.62e-17 1e-8 1.00e+00 8.33e-18 ``` **検算.** (a) 既知値: $(1-\cos x)/x^2\to1/2$ で、参照値は同じ binary64 の $x$ に対する 50 桁の値です。(b) 別精度・別の式: 書き換えた式の誤差は全行で $u\approx1.1\times10^{-16}$ 程度で、参照値と 15〜16 桁一致します。(c) 次数: 素朴な式の誤差は上の見積もり $u/(x^2/2)=2u/x^2$ 以下で、$x$ を $1/10$ にするとほぼ $100$ 倍に増えます（$10^{-4}$ から $10^{-6}$ で $5.2\times10^{-9}\to8.9\times10^{-5}$）。$x=10^{-8}$ では $x^2/2\lt u$ なので $\mathrm{fl}(\cos x)=1$ となり、答えは $0$、相対誤差は $1$ です。なお NumPy の `cos`, `sin` は正しく丸められる保証はありませんが、誤差は数 $u$ 程度で、上の見積もりの大きさを変えません。 ## 5. 和の誤差 ### 5.1. 逐次和 $n$ 個の浮動小数点数 $x_1,\ldots,x_n$ を前から順に足すアルゴリズム $$ \hat s_1=x_1,\qquad\hat s_k=\mathrm{fl}(\hat s_{k-1}+x_k)\quad(k=2,\ldots,n) $$ を**逐次和**といいます。 $x_1,\ldots,x_n\in\mathbb{F}$、$(n-1)u\lt1$ とし、逐次和の各段でオーバーフローもアンダーフローも起きないとする。このとき $$ \hat s_n=\sum_{i=1}^nx_i(1+\theta^{(i)}),\qquad|\theta^{(i)}|\le\gamma_{n-1}, $$ したがって $s=\sum_ix_i$ について $|\hat s_n-s|\le\gamma_{n-1}\sum_{i=1}^n|x_i|$ である。 標準モデル（）から $\hat s_k=(\hat s_{k-1}+x_k)(1+\delta_k)$, $|\delta_k|\le u$（$k=2,\ldots,n$）です。これを $k=n$ から順に代入して展開すると $$ \hat s_n=x_1\prod_{j=2}^n(1+\delta_j)+\sum_{i=2}^nx_i\prod_{j=i}^n(1+\delta_j) $$ です。$x_1$ に掛かる因子は $n-1$ 個、$x_i$（$i\ge2$）に掛かる因子は $n-i+1\le n-1$ 個です。誤差因子の積の補題（）から $k$ 個の積は $1+\theta_k$, $|\theta_k|\le\gamma_k$ で、 の 1 から $\gamma_k\le\gamma_{n-1}$ です。最後の不等式は $\hat s_n-s=\sum_ix_i\theta^{(i)}$ と三角不等式からです。 この評価は 2 つのことを言っています。第一に、計算した和は**入力を相対 $\gamma_{n-1}$ だけ動かした数の正確な和**です（後退誤差が小さい）。第二に、正の数の和なら $\sum|x_i|=s$ なので相対誤差は $\gamma_{n-1}\approx(n-1)u$ 以下ですが、正負が打ち消し合う和では $\sum|x_i|/|s|$ 倍に拡大されます（ と同じ構造）。 証明の展開をよく見ると、$x_i$ に掛かる因子は $n-i+1$ 個で、先に足す項ほど多くの丸めを受けます。正の数を足すときは小さい順に足すと評価が小さくなります（演習で確かめます）。 ### 5.2. 対ごとの和 **対ごとの和**（pairwise summation）は、隣どうしを足して長さを半分にする操作を繰り返します。長さ $m$ の列 $(y_1,\ldots,y_m)$ から、$m$ が奇数なら末尾に $0$ を足して偶数にし、$(\mathrm{fl}(y_1+y_2),\mathrm{fl}(y_3+y_4),\ldots)$ を作ります。これを長さが $1$ になるまで繰り返し、残った数を $\hat s$ とします。 $x_1,\ldots,x_n\in\mathbb{F}$（$n\ge2$）、$L=\lceil\log_2n\rceil$、$Lu\lt1$ とし、途中でオーバーフローもアンダーフローも起きないとする。対ごとの和の結果 $\hat s$ は $$ \hat s=\sum_{i=1}^nx_i(1+\theta^{(i)}),\qquad|\theta^{(i)}|\le\gamma_L $$ を満たし、$|\hat s-s|\le\gamma_{\lceil\log_2n\rceil}\sum_i|x_i|$ である。 **段数.** $\ell$ 段の後の列の長さは $\lceil n/2^\ell\rceil$ です。$\ell$ についての帰納法で、長さ $m=\lceil a\rceil$（$a=n/2^\ell$）の列が次の段で $\lceil m/2\rceil$ になるので、$\lceil m/2\rceil=\lceil a/2\rceil$ を示せばよいです。$a\le m$ から $\lceil a/2\rceil\le\lceil m/2\rceil$。逆に $2\lceil a/2\rceil$ は $a$ 以上の整数なので $m=\lceil a\rceil$ 以上で、$\lceil a/2\rceil\ge m/2$、整数なので $\lceil a/2\rceil\ge\lceil m/2\rceil$ です。長さが $1$ になるのは $n\le2^\ell$ のとき、すなわち $\ell\ge\log_2n$ のときなので、段数は $L$ です。 **誤差.** $\ell$ 段の後の各成分が、元の連続したいくつかの $x_i$ について $\sum x_i(1+\theta^{(i)}_\ell)$ の形（$1+\theta^{(i)}_\ell$ は $\ell$ 個の $(1+\delta)$ の積）であることを $\ell$ の帰納法で示します。$\ell=0$ は $x_i$ 自身です。次の段で 2 つの成分を足すと、標準モデル（）から和全体に 1 つの因子 $(1+\delta)$ が掛かり、各 $x_i$ の因子は 1 個増えます。末尾に補った $0$ との和は厳密で、$\delta=0$ とすれば同じ形です。$L$ 段の後、各 $x_i$ には $L$ 個の因子が掛かり、誤差因子の積の補題（）から $|\theta^{(i)}|\le\gamma_L$ です。最後の不等式は三角不等式です。 $n=10^6$ で、逐次和の評価の係数は $n-1\approx10^6$、対ごとの和は $\lceil\log_2 10^6\rceil=20$ です。演算回数はどちらも $n-1$ 回の足し算です。NumPy の `np.sum`（メモリ上で連続な軸に沿った和）は、変形した対ごとの和を使っています。列を（前半の長さを 8 の倍数にそろえて）半分に分けることを長さが 128 以下になるまで再帰的に繰り返し、128 個以下のかたまりの中では 8 本の逐次和を並行に走らせてから 8 本を対ごとに足し、かたまりどうしは分けたときと逆の順に対ごとに足します。かたまりの中に逐次和が残るので、上の定理の評価はそのままは当てはまりません。 乱数種を固定した $[0,1)$ の binary32 の一様乱数 $10^6$ 個について、先頭 $n$ 個の和を 3 通りで計算し、正しく丸めた和（`math.fsum` を binary64 で）と比べます。表の列は、$n$、逐次和の相対誤差とその評価 $\gamma_{n-1}$、自作の対ごとの和の相対誤差とその評価 $\gamma_{\lceil\log_2n\rceil}$、`np.sum` の相対誤差です（$u=2^{-24}$）。 ```python # Python 3.12.10 / NumPy 2.5.2（x86-64）／乱数種 20261006 import math from fractions import Fraction import numpy as np def pairwise(x): # 隣どうしを足して長さを半分にする while len(x) > 1: if len(x) % 2: x = np.append(x, np.float32(0)) # 0 を足すのは厳密 x = x[0::2] + x[1::2] return x[0] u = 2.0 ** -24 # binary32 の単位丸め g = lambda k: k * u / (1 - k * u) rng = np.random.default_rng(20261006) xs = rng.random(10**6, dtype=np.float32) for n in (10**3, 10**4, 10**5, 10**6): x = xs[:n] s = math.fsum(x.astype(np.float64)) # 正しく丸めた和 rec = np.cumsum(x, dtype=np.float32)[-1] # 逐次和 e = [abs(float(t) - s) / s for t in (rec, pairwise(x), np.sum(x))] k = math.ceil(math.log2(n)) print(f"{n:>8} {e[0]:.2e} {g(n - 1):.2e} {e[1]:.2e} {g(k):.2e}" f" {e[2]:.2e}") ex = sum(map(Fraction, xs.tolist())) # 厳密な和（有理数） print(f"fsum と厳密な和の相対差 {abs(s - ex) / ex:.2e}") ``` ```text 1000 6.22e-07 5.95e-05 7.56e-09 5.96e-07 7.56e-09 10000 1.16e-06 5.96e-04 9.95e-09 8.34e-07 1.08e-07 100000 4.52e-06 6.00e-03 1.47e-08 1.01e-06 1.47e-08 1000000 1.54e-05 6.34e-02 1.71e-08 1.19e-06 1.71e-08 fsum と厳密な和の相対差 0.00e+00 ``` **検算.** (a) 既知値: 各項は $2^{-24}$ の整数倍で和は $2^{20}$ 未満なので、厳密な和は 44 ビットに収まり binary64 で正確に表せます。最終行のとおり `fsum` は有理数で計算した厳密な和と一致し、参照値に丸め誤差はありません。(b) 別手法: 逐次和と自作の対ごとの和の誤差は対応する評価の内側にあり、自作の対ごとの和は逐次和より 2〜3 桁小さい誤差です。`np.sum` は $n=10^4$ を除いて自作の対ごとの和と同じ誤差で、$n=10^4$ でも逐次和より 1 桁小さく、$\gamma_{\lceil\log_2n\rceil}$ の内側にあります（`np.sum` は足す順序が違うので、この表の値は定理の対象外の参考値です）。(c) 次数: 評価 $\gamma_{n-1}$ は $n$ を 10 倍するごとに 10 倍、$\gamma_{\lceil\log_2n\rceil}$ は $\log_2 10\approx3.3$ ずつ係数が増えるだけです。実際の逐次和の誤差は 10 倍ごとに 2〜4 倍で、最悪評価の 10 倍よりゆるやかです。丸め誤差の符号がそろわず打ち消し合うためで、定理はこの典型的な大きさまでは主張しません。 調和級数 $\sum1/k$ は発散しますが、binary32 で前から足すと和はある所で増えなくなります。和 $\hat s$ が $[8,16)$ にあるとき間隔は $2u\cdot8=2^{-20}$ なので、足す項 $1/k$ が間隔の半分 $2^{-21}$ 以下になると、$\hat s+1/k$ は $\hat s$ に丸められます（ちょうど半分のときは偶数丸めで、仮数が偶数ならそのまま）。和が $16$ に届く前に $k=2^{21}$ に達するなら、そこで止まるはずです。$\ln 2^{21}+0.5772\approx15.13\lt16$ なので、そう見積もれます。 ```python # Python 3.12.10 / NumPy 2.5.2（x86-64） import numpy as np k = np.arange(1, 2**23 + 1, dtype=np.float32) t = np.float32(1) / k # 1/k を binary32 に丸めた項 s = np.cumsum(t, dtype=np.float32) # 前から逐次に足す m = np.flatnonzero(s[1:] == s[:-1])[0] + 2 # 和が初めて増えない項 print("止まる項", m, "= 2^21:", m == 2**21) print("binary32 の和", f"{s[m - 1]:.6f}", "最後まで", f"{s[-1]:.6f}") ref = np.sum(1.0 / np.arange(1, m + 1, dtype=np.float64)) print("binary64 の和", f"{ref:.6f}") ``` ```text 止まる項 2097152 = 2^21: True binary32 の和 15.403683 最後まで 15.403683 binary64 の和 15.133307 ``` **検算.** (a) 既知値: 止まる項は見積もりどおり $2^{21}=2097152$ で、以後 $2^{23}$ 項まで和は変わりません。(b) 別精度: binary64 で同じ $2^{21}$ 項を足すと $15.133307$ で、$\ln2^{21}+0.577216=15.133306$（$0.577216$ は Euler の定数）とほぼ一致します（差は $1/(2\cdot2^{21})\approx2.4\times10^{-7}$ 程度の補正）。binary32 の和 $15.4037$ はこれより $0.27$ 大きくなります。主な原因は最後の段階 $2^{20}\lt k\le2^{21}$ です。ここでは和が $[8,16)$ にあり、項 $1/k$ は間隔 $2^{-20}$ の半分と 1 間隔の間にあるので、どの項も 1 目盛り $2^{-20}$ に切り上げられて足されます。過大分は $\sum_{k=2^{20}+1}^{2^{21}}(2^{-20}-1/k)\approx1-\ln2\approx0.307$ です。それより前の段階では切り上げと切り捨てが混じり、差し引きで約 $0.037$ 小さくなるので（各段の丸め誤差を足し上げて確かめられます）、合計が $0.27$ です。(c) 誤差の大きさ: 差 $0.27$ は の評価 $\gamma_{n-1}\sum|x_i|$（$n=2^{21}$ では $(n-1)u\approx0.125$ で、$\gamma\approx0.14$、評価は約 $2.2$）の内側です。 ### 5.3. 補償和と融合積和 和の誤差をさらに小さくする方法として、各段の丸め誤差を別の変数に拾って次の段で足し戻す**補償和**（Kahan の方法）があります。 Kahan の補償和で計算した和 $\hat s$ は $|\hat s-s|\le(2u+O(nu^2))\sum_i|x_i|$ を満たすことが知られています。係数が $n$ によらないのが要点です。証明は各段の誤差の場合分けが長いので、この記事では扱いません（Higham『Accuracy and Stability of Numerical Algorithms』第 2 版 §4.3、Knuth『The Art of Computer Programming』第 2 巻 §4.2.2）。Python 3.12 以降の組み込み `sum` は浮動小数点数に対して補償つきの加算を使い、`math.fsum` は正しく丸めた和を返します。以下の証明ではこの評価を使いません。 融合積和 $\mathrm{fma}(a,b,c)=\mathrm{fl}(ab+c)$ は積と和を 1 回の丸めで計算します。正しく丸められることは IEEE 754 が要求するので、標準モデル（）と同じ形 $(ab+c)(1+\delta)$ が成り立ち、内積や Horner 法の 1 段の丸めが 2 回から 1 回に減ります。さらに $ab-\mathrm{fl}(ab)=\mathrm{fma}(a,b,-\mathrm{fl}(ab))$ が（アンダーフローがなければ）厳密に計算できるので、積の丸め誤差を正確に取り出せます。これを使う誤差なし変換と倍々精度の算法は Muller ほか『Handbook of Floating-Point Arithmetic』第 2 版の該当章を参照してください。この記事の定理は fma を仮定しません。 ## 6. Horner 法の後退誤差 多項式 $p(x)=\sum_{i=0}^na_ix^i$ の値を計算する標準的な方法は **Horner 法** $$ q_n=a_n,\qquad q_k=x\,q_{k+1}+a_k\quad(k=n-1,\ldots,0),\qquad p(x)=q_0 $$ です。乗算 $n$ 回・加算 $n$ 回で済みます。浮動小数点数では各段が $\hat q_k=\mathrm{fl}(\mathrm{fl}(x\hat q_{k+1})+a_k)$ になります。 $x,a_0,\ldots,a_n\in\mathbb{F}$、$2nu\lt1$ とし、Horner 法の途中でオーバーフローもアンダーフローも起きないとする。このとき計算値 $\hat q_0$ は $$ \hat q_0=\sum_{i=0}^na_i(1+\theta_i)\,x^i,\qquad|\theta_i|\le\gamma_{2n} $$ を満たす。すなわち $\hat q_0$ は、各係数を相対 $\gamma_{2n}$ 以内で動かした多項式 $\tilde p(x)=\sum_i\tilde a_ix^i$（$\tilde a_i=a_i(1+\theta_i)$）の、$x$ における正確な値である。 証明の方針は、各段の乗算と加算の丸めを、それぞれ $(1+\varepsilon_k)$, $(1+\delta_k)$ という因子として書き、最後に係数 $a_i$ へ掛かる因子の個数を数えることです。添字は多いですが、$a_i$ に掛かる因子が高々 $2n$ 個であることを示すための帳簿にすぎません。 標準モデル（）から、$k=n-1,\ldots,0$ について $$ \hat q_k=\bigl(x\hat q_{k+1}(1+\varepsilon_k)+a_k\bigr)(1+\delta_k),\qquad|\varepsilon_k|,|\delta_k|\le u $$ です。$k$ についての下向きの帰納法で $$ \hat q_k=\sum_{i=k}^na_ix^{i-k}\,\pi_{k,i},\qquad\pi_{k,i}=(1+\delta_i)^{[i\lt n]}\prod_{j=k}^{i-1}(1+\varepsilon_j)(1+\delta_j) $$ を示します（$[i\lt n]$ は $i\lt n$ なら $1$、$i=n$ なら $0$）。$k=n$ では $\hat q_n=a_n$ で $\pi_{n,n}=1$ です。$k+1$ で成り立てば、上の式に代入して $$ \hat q_k=\sum_{i=k+1}^na_ix^{i-k}\,\pi_{k+1,i}(1+\varepsilon_k)(1+\delta_k)+a_k(1+\delta_k) $$ で、$\pi_{k+1,i}(1+\varepsilon_k)(1+\delta_k)=\pi_{k,i}$, $1+\delta_k=\pi_{k,k}$ なので $k$ でも成り立ちます。$\pi_{0,i}$ の因子の個数は、$i\lt n$ なら $1+2i\le2n-1$、$i=n$ なら $2n$ です。誤差因子の積の補題（）と の 1 から $\pi_{0,i}=1+\theta_i$, $|\theta_i|\le\gamma_{2n}$ です。 この定理は、Horner 法の答えが「正しい問題の近似解」ではなく「近い問題の正確な解」であることを言っています。係数そのものが測定値や前の計算の結果で相対 $u$ 程度の誤差をもつなら、Horner 法の誤差はそれと同じ程度で、アルゴリズムとしてはこれ以上望めません。答え $p(x)$ そのものの誤差は、ここから三角不等式で $$ |\hat q_0-p(x)|\le\gamma_{2n}\sum_{i=0}^n|a_i||x|^i $$ と評価できます。右辺を $|p(x)|$ で割った倍率 $\sum|a_i||x|^i/|p(x)|$ は、$p(x)$ が根の近くで小さいときに大きくなります。根の近くで多項式の値を正確に計算できないのは、Horner 法が悪いのではなく、問題そのものが係数の摂動に敏感だからです。この「問題の感度（条件数）」と「アルゴリズムの後退誤差」の区別が、このテーマの[条件数と安定性](/mathematics/numerical-analysis/conditioning-and-stability)の章の主題です。 ## 7. 演習 binary64 の $[1,2)$ に含まれる浮動小数点数の個数を求めよ。また、正規化数の範囲で隣り合う浮動小数点数の相対的な間隔 $(\mathrm{succ}(x)-x)/x$ の最大値と最小値を、$u$ を使って表せ。 $[1,2)$ の浮動小数点数は $M\cdot2^{-52}$（$2^{52}\le M\lt2^{53}$）なので $2^{52}\approx4.5\times10^{15}$ 個です。$x\in[2^e,2^{e+1})$ の正規化数の右の間隔は $2u\cdot2^e$ なので、相対的な間隔は $2u\cdot2^e/x$ で、$x=2^e$ のとき最大値 $2u$、$x$ が $2^{e+1}$ の直前（$x=(2-2u)2^e$）のとき最小値 $2u/(2-2u)=u/(1-u)$ です。相対的な間隔が 2 倍近く揺れることが、$\varepsilon_M=2u$ と $u$ の 2 つの量が使われる理由の 1 つです。 Sterbenz の補題の条件 $y/2\le x\le2y$ はゆるめられないことを、$x$ を $2$ の次の浮動小数点数、$y$ を $1$ の前の浮動小数点数として確かめよ（精度 $p$ の一般の形式で、$e_{\min}\le-1$, $e_{\max}\ge1$ とする）。 $[2,4)$ の間隔は $4u$、$[1/2,1)$ の間隔は $u$ なので、$x=2+4u$, $y=1-u$ です。$x/y=(2+4u)/(1-u)\gt2$ で条件をわずかに外れます。$x-y=1+5u$ ですが、$[1,2)$ の浮動小数点数は $2u$ の整数倍なので、奇数倍の $5u$ を含む $1+5u$ は $\mathbb{F}$ に属しません。したがって $\mathrm{fl}(x-y)\ne x-y$ で、引き算に丸め誤差が入ります。 の証明の展開を使って、$|\hat s_n-s|\le\sum_{i=1}^n\gamma_{c_i}|x_i|$（$c_1=n-1$、$c_i=n-i+1$（$i\ge2$））を示せ。これを使って、正の数 $x_1,\ldots,x_n$ を足すとき、小さい順に並べて足すとこの評価が最小になることを説明せよ。 展開の式で $x_i$ に掛かる因子は $c_i$ 個なので、誤差因子の積の補題から $|\theta^{(i)}|\le\gamma_{c_i}$ で、三角不等式から主張が出ます。係数の列 $\gamma_{c_1},\gamma_{c_2},\ldots,\gamma_{c_n}$ は $\gamma_{n-1},\gamma_{n-1},\gamma_{n-2},\ldots,\gamma_1$ で、添字が後ろほど小さい非増加列です。$\sum\gamma_{c_i}x_i$ を最小にする並べ方は、並べ替え不等式（大きい係数に小さい数を対応させる）により $x_i$ を昇順に並べたものです。ただしこれは**評価**を最小にするだけで、実際の誤差が必ず最小になるとは限りません。 $p(x)=(x-1)^{10}$ を展開した係数で Horner 法により $x=1.01$ での値を計算するとき、 から得られる相対誤差の評価 $\gamma_{20}\sum|a_i||x|^i/|p(x)|$ を見積もれ。binary64 で何桁が信頼できるか。 $a_i=\binom{10}{i}(-1)^{10-i}$ なので $\sum|a_i||x|^i=(1+x)^{10}=2.01^{10}\approx1.08\times10^{3}$、$|p(x)|=0.01^{10}=10^{-20}$ です。倍率は約 $1.1\times10^{23}$ で、$\gamma_{20}\approx20u\approx2.2\times10^{-15}$ を掛けると評価は $2\times10^{8}$ になり、1 桁も保証されません。展開せずに $(x-1)^{10}$ と計算すれば、$x-1$ は により厳密で、相対誤差は $\gamma_{10}$ 程度です（$x$ を浮動小数点数とした場合）。同じ関数でも、計算する式によって答えの精度がまったく違います。 ## 参考文献 - N. J. Higham, *Accuracy and Stability of Numerical Algorithms*, 2nd ed., SIAM, 2002 — 第 2 章（浮動小数点数と標準モデル）、第 3 章（$\gamma_n$）、第 4 章（和の誤差と補償和）、第 5 章（多項式と Horner 法）。 - IEEE Computer Society, *IEEE Standard for Floating-Point Arithmetic* (IEEE Std 754-2019), 2019 — 第 3 章（形式）、第 4 章（丸めの属性）、第 5 章（演算）。 - D. Goldberg, "What every computer scientist should know about floating-point arithmetic", *ACM Computing Surveys* 23 (1991), 5–48. [doi:10.1145/103162.103163](https://doi.org/10.1145/103162.103163) - J.-M. Muller ほか, *Handbook of Floating-Point Arithmetic*, 2nd ed., Birkhäuser, 2018 — 浮動小数点数の形式、fma、誤差なし変換。 - P. H. Sterbenz, *Floating-Point Computation*, Prentice-Hall, 1974 — 引き算が厳密になる条件。 - L. N. Trefethen, D. Bau III, *Numerical Linear Algebra*, SIAM, 1997 — Lecture 13（浮動小数点数）、Lecture 14（安定性）。