浮動小数点数と丸め誤差：IEEE 754・標準モデル・和と Horner 法の誤差解析

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}$ 以内で動かした多項式の正確な値です（後退誤差）。この「計算結果は少し違う入力に対する正確な答え」という見方が、数値解析と数値線形代数の誤差解析の出発点になります。

この記事で使う事前知識

• 仮数と指数: 数を「仮数 × 2 の指数乗」の形で書くときの、仮数（有効数字にあたる部分）と指数のこと。→ 定義 2.1
• オーバーフローとアンダーフロー: 結果が表せる最大の値を超えること、または 0 でないのに表せる最小の正規化数より小さくなること。→ 定義 2.1［浮動小数と区間演算］
• 相対誤差: 誤差の大きさを、真の値の大きさで割ったもの。桁数の感覚に対応します。
• Horner 法: 多項式 $a_nx^n+\cdots+a_0$ を $(\cdots(a_nx+a_{n-1})x+\cdots)x+a_0$ の順に、掛け算と足し算だけで計算する方法（第 6 節で定義します）。
• 条件数: 入力が少し変わったとき、答えがどれだけ大きく変わりうるかを表す量（問題の側の性質）。→ 条件数（相対条件数）

1. 動機

Python で 0.1 + 0.2 == 0.3 を評価すると False になります。数値計算を始めた人がまず出会う現象ですが、これは計算機の「誤り」ではありません。十進の $0.1$ は 2 進では有限桁で書けないので、読み込んだ時点で最も近い 2 進数に置き換わり、足し算の結果もまた最も近い 2 進数に置き換わります。どの段階でどれだけずれるかは、規格によって正確に決まっています。

数値解析の仕事は、この「置き換え」が積み重なったときに答えがどれだけずれるかを、定理として評価することです。ずれには 2 種類あります。

種類	原因	例
丸め誤差	実数を有限桁の浮動小数点数で表すこと	和・内積・多項式の値の誤差
打ち切り誤差	無限の操作（極限・級数・積分）を有限で止めること	差分商、数値積分、反復の停止

この記事は前者を扱います。丸めの基本的な誤差モデルと誤差因子の積の補題は、計算機援用証明のテーマの浮動小数と区間演算の章で証明されているので、ここではその主張を引用して使い、IEEE 754 の形式と、和・桁落ち・多項式の誤差解析を積み上げます。

2. IEEE 754 の 2 進形式

2.1. 浮動小数点数の集合

精度 $p$ と指数の範囲 $e_{\min}\lt e_{\max}$ を固定した浮動小数点数の集合 $\mathbb{F}$ は、既存の定義（定義 2.1［浮動小数と区間演算］）を使います。要点を再掲します。

• 正規化数は $\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［浮動小数と区間演算］）。間隔は $2$ の冪をまたぐたびに 2 倍になり、相対的な間隔は $u$ から $2u$ の間にあります。

1 の近くの浮動小数点数（p = 4 の場合）。[1/2, 1) では間隔 2⁻⁴、[1, 2) では 2⁻³、[2, 4) では 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$ を表します。

binary64 のビット配置。符号 1 ビット、指数部 E が 11 ビット、仮数部 f が 52 ビット。

定義 2.1（IEEE 754 の 2 進交換形式）↗

正の整数 $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. 丸めの方向と正しい丸め

定義 2.2（丸めの方向と正しい丸め）↗

$|x|\le x_{\max}$ の実数 $x$ について、次の 4 つの丸めを考える。

• 最近接偶数丸め $\mathrm{fl}(x)$: $x$ に最も近い $\mathbb{F}$ の元。2 つあれば $M$ が偶数の方（定義 2.3［浮動小数と区間演算］ の $\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 などの初等関数については、規格は正しい丸めを推奨するにとどめているので、以下の評価をそのまま適用できません。

定義 2.3（機械イプシロンと単位丸め）↗

$1$ と、$1$ より大きい最小の浮動小数点数との差 $\varepsilon_M=2^{1-p}$ を機械イプシロンという。その半分 $u=2^{-p}$ を単位丸めという。

$\varepsilon_M$ は $1$ の右側の間隔で、$u$ は最近接丸めの相対誤差の上界です（定理 2.4［浮動小数と区間演算］）。$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}$ を使います。

例 2.4（NumPy で単位丸めと間隔を確かめる）

np.finfo と np.nextafter で、binary32 と binary64 の精度・指数の範囲・1 の両側の間隔・最小の正規化数と非正規化数・最大値を出します。

# 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}")

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 倍になること（補題 2.2［浮動小数と区間演算］）のとおりです。

3. 四則の標準モデル

正しい丸めの規約と、丸めの誤差モデル（定理 2.4［浮動小数と区間演算］）を組み合わせると、すべての誤差解析の出発点になる次の形が得られます。

系 3.1（四則の標準モデル）✓

$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})$ についても同じことが成り立つ。

証明（系 3.1）

正しい丸めの規約（定義 2.2）により、計算結果は $\mathrm{fl}(z)$ です。$z\in\mathbb{F}$（$z=0$ を含む）なら丸めの定義から $\mathrm{fl}(z)=z$ で、$\delta=\delta'=0$ とします。そうでなければ $x_{\min}\le|z|\le x_{\max}$ なので、丸めの誤差モデル（定理 2.4［浮動小数と区間演算］ の 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$。定義 3.1［浮動小数と区間演算］）と、誤差因子の積の補題を使います。

誤差因子の積の補題（補題 3.2［浮動小数と区間演算］）: $|\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$ の計算規則を、補題として証明しておきます。数値線形代数の後退誤差解析で繰り返し使います。

補題 3.2（γₙ の計算規則）✓

$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}$。

証明（補題 3.2）

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}$ です。

∎

例 3.3（0.1 + 0.2 が 0.3 にならない理由）

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}$ で、系 3.1 の範囲 $u\,|3M\cdot2^{-55}|\approx3.3\times10^{-17}$ に入っています。

4. 引き算：Sterbenz の補題と桁落ち

4.1. 誤差なく計算される引き算

補題 4.1（Sterbenz の補題）✓

$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$ 桁に収まる、というものです。非正規化数の場合分けが入るために長く見えますが、中心の考えはこれだけです。

証明（補題 4.1）

仮定から $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})$ とおきます。定義 2.1［浮動小数と区間演算］ から、$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. 桁落ち

それでも近い数の引き算が危険なのは、引く前の値がすでに誤差を含んでいるからです。

命題 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$ でとられる。

証明（命題 4.2）

$(\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）といいます。補題 4.1 と合わせると、桁落ちの誤差は引き算の前に入ったもので、引き算はそれを露わにするだけだとわかります。対策は、引き算の前に誤差の入らない式に書き換えることです。

例 4.3（(1 + x) − 1）

$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]$ なので 補題 4.1 により引き算は厳密です。したがって $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 桁が失われます。これは 命題 4.2 の倍率 $(|1+x|+1)/x$ と同じ大きさです。

例 4.4（1 − cos 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 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}")

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)$$

を逐次和といいます。

定理 5.1（逐次和の誤差）✓

$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|$ である。

証明（定理 5.1）

標準モデル（系 3.1）から $\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$ 個です。誤差因子の積の補題（補題 3.2［浮動小数と区間演算］）から $k$ 個の積は $1+\theta_k$, $|\theta_k|\le\gamma_k$ で、補題 3.2 の 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|$ 倍に拡大されます（命題 4.2 と同じ構造）。

証明の展開をよく見ると、$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$ とします。

定理 5.2（対ごとの和の誤差）✓

$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|$ である。

証明（定理 5.2）

段数. $\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 つの成分を足すと、標準モデル（系 3.1）から和全体に 1 つの因子 $(1+\delta)$ が掛かり、各 $x_i$ の因子は 1 個増えます。末尾に補った $0$ との和は厳密で、$\delta=0$ とすれば同じ形です。$L$ 段の後、各 $x_i$ には $L$ 個の因子が掛かり、誤差因子の積の補題（補題 3.2［浮動小数と区間演算］）から $|\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 本を対ごとに足し、かたまりどうしは分けたときと逆の順に対ごとに足します。かたまりの中に逐次和が残るので、上の定理の評価はそのままは当てはまりません。

例 5.3（binary32 の 10⁶ 個の和）

乱数種を固定した $[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 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}")

    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 倍よりゆるやかです。丸め誤差の符号がそろわず打ち消し合うためで、定理はこの典型的な大きさまでは主張しません。

例 5.4（binary32 の調和級数はどこで止まるか）

調和級数 $\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 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}")

止まる項 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$ は 定理 5.1 の評価 $\gamma_{n-1}\sum|x_i|$（$n=2^{21}$ では $(n-1)u\approx0.125$ で、$\gamma\approx0.14$、評価は約 $2.2$）の内側です。

5.3. 補償和と融合積和

和の誤差をさらに小さくする方法として、各段の丸め誤差を別の変数に拾って次の段で足し戻す補償和（Kahan の方法）があります。

注意 5.5（補償和の誤差（主張のみ・文献））

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 は正しく丸めた和を返します。以下の証明ではこの評価を使いません。

注意 5.6（融合積和の効用（文献））

融合積和 $\mathrm{fma}(a,b,c)=\mathrm{fl}(ab+c)$ は積と和を 1 回の丸めで計算します。正しく丸められることは IEEE 754 が要求するので、標準モデル（系 3.1）と同じ形 $(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)$ になります。

定理 6.1（Horner 法の後退誤差）✓

$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$ 個であることを示すための帳簿にすぎません。

証明（定理 6.1）

標準モデル（系 3.1）から、$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$ です。誤差因子の積の補題（補題 3.2［浮動小数と区間演算］）と 補題 3.2 の 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 法が悪いのではなく、問題そのものが係数の摂動に敏感だからです。この「問題の感度（条件数）」と「アルゴリズムの後退誤差」の区別が、このテーマの条件数と安定性の章の主題です。

7. 演習

演習 7.1易

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 つです。

演習 7.2標準

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$ で、引き算に丸め誤差が入ります。

演習 7.3標準

定理 5.1 の証明の展開を使って、$|\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$ を昇順に並べたものです。ただしこれは評価を最小にするだけで、実際の誤差が必ず最小になるとは限りません。

演習 7.4標準

$p(x)=(x-1)^{10}$ を展開した係数で Horner 法により $x=1.01$ での値を計算するとき、定理 6.1 から得られる相対誤差の評価 $\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$ は 補題 4.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
• 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（安定性）。