# 精度の選び方：倍精度・任意精度・球算術と誤差予算

> 計算値の誤差を丸め・打ち切り・データの 3 項に分ける誤差予算を定式化し、条件数から必要な作業桁数を決める定理、桁落ちで失う桁数、精度を倍々にする Newton 法の費用、区間・球算術と 2 精度比較の保証の違いを証明し、mpmath の実行例で確かめる。
> https://rikai.mugen-giken.com/computer-science/research-computing/precision-and-error-budget

## 0. この記事の要点

計算機で出した数値を「小数点以下 10 桁まで正しい」と言うには、どこで誤差が入りうるかを数えて、足し合わせておく必要があります。この記事は、誤差の入り口を整理し、目標の桁数を得るには計算を何桁で行えばよいかを決め、得られた桁が本当に正しいかを確かめる方法を扱います。読み終えると、「倍精度で足りるのか、桁数を増やすべきか」を勘ではなく見積もりで判断できるようになります。

- 数値の主張に桁を付けるには、誤差を原因ごとに分けて上から抑えた**誤差予算**を作ります。計算値と真値の差は、丸め・打ち切り・入力データの 3 項の和で抑えられます。
- [相対条件数](/mathematics/numerical-analysis/conditioning-and-stability#def-condition-number) $\kappa$ の問題を後退誤差（入力をわずかにずらすと計算値がちょうど正しい答えになる、というときのずれの大きさ）$cu$ の算法で解くと、相対誤差はほぼ $c\kappa u$ です。目標の $d$ 桁を得るのに必要な作業桁数は $d+\log_{10}\kappa+\log_{10}c$ で足ります。**条件数の桁だけ余分に持つ**、が精度の選び方の基本です。
- 打ち消し合う和 $\sum x_i$ は $\log_{10}(\sum|x_i|/|\sum x_i|)$ 桁を失います。交代級数で $e^{-x}$ を計算すると約 $2x/\ln10\approx0.87x$ 桁を失い、$x=20$ で倍精度の 16 桁をすべて使い切ります。
- 2 次収束する反復では、作業精度を反復ごとに倍にしていけば、全体の費用は最後の 1 回の 2 倍以下で済みます。
- 区間演算・球算術（答えを 1 つの数ではなく、範囲として返す計算法）の結果は真値を含みます。浮動小数の「2 つの精度で一致した桁」は丸め誤差の推定にはなりますが、包含の保証にはならず、打ち切り誤差は見えません。
- 高精度計算を台無しにする最も多い原因は、定数や入力を倍精度で作ってから渡すことです。高精度の数は文字列か有理数から作ります。

<Aside type="tip" title="この記事で使う事前知識">

- **後退誤差・前進誤差**: 前進誤差は答えのずれ、後退誤差は「入力をどれだけずらせば計算値が正確な答えになるか」というずれ。→ [条件数と安定性](/mathematics/numerical-analysis/conditioning-and-stability)
- **相対条件数**: 入力の相対的な変化が、出力の相対的な変化を何倍にするかを表す数。→ [相対条件数](/mathematics/numerical-analysis/conditioning-and-stability#def-condition-number)
- **単位丸め $u$**: 浮動小数点数で 1 回の演算をしたときの、相対的な丸め誤差の上限。→ [浮動小数点数と丸め誤差](/mathematics/numerical-analysis/floating-point-arithmetic)
- **binary64 / binary32**: 仮数部がそれぞれ 53 ビット・24 ビットの浮動小数点数の形式（倍精度・単精度）。→ [浮動小数点数と丸め誤差](/mathematics/numerical-analysis/floating-point-arithmetic)
- **桁落ち**: 近い数の引き算で、有効な桁が減ること。→ [浮動小数点数と丸め誤差](/mathematics/numerical-analysis/floating-point-arithmetic)

</Aside>

## 1. 動機

「この積分の値は $0.4872$ である」と書くとき、最後の $2$ は正しいでしょうか。倍精度で計算したから 16 桁正しい、とは言えません。刻み幅による打ち切りで 5 桁目から違うかもしれませんし、近い数の引き算で 10 桁を失っているかもしれませんし、入力の定数が 3 桁しか正しくないかもしれません。

逆の誤りもあります。念のため任意精度で 50 桁計算したのに、入力の $0.1$ を倍精度の数として渡していたために、17 桁目から先は「$0.1$ ではない別の数」についての正確な計算になっている、という場合です。

精度の選び方は、次の 3 つの問いに順に答えることです。

1. 誤差はどこから入るか（誤差予算）。
2. 目標の桁数を得るには、作業精度を何桁にすればよいか（条件数）。
3. 得られた桁を、どうやって確かめるか（別精度での再計算、区間・球算術）。

条件数・前進誤差・後退誤差は[条件数と安定性](/mathematics/numerical-analysis/conditioning-and-stability)の章の定義を、浮動小数点数の誤差モデルは[浮動小数点数と丸め誤差](/mathematics/numerical-analysis/floating-point-arithmetic)の章のものを使います。区間演算と球算術の定義と包含の定理は[浮動小数と区間演算](/mathematics/computer-assisted-proofs/floating-point-interval-arithmetic)の章にあり、ここではその使い方の側を扱います。

## 2. 誤差予算

### 2.1. 用語

<Definition id="def-precision-terms" title="作業精度・目標精度・有効桁">
- 計算に使う浮動小数点数の精度 $p$（ビット）を**作業精度**といい、$D=-\log_{10}u=p\log_{10}2$（$u=2^{-p}$ は単位丸め。<Ref to="mathematics/numerical-analysis/floating-point-arithmetic#def-machine-epsilon" />）を**作業桁数**という。binary64 では $D\approx15.95$、binary32 では $D\approx7.22$ である。
- 最終結果に求める相対誤差の上限 $10^{-d}$ を**目標精度**、$d$ を目標の桁数という。
- 真値 $y\ne0$ に対する計算値 $\hat y$ の**有効桁**を $-\log_{10}(|\hat y-y|/|y|)$ と定める。
</Definition>

<Definition id="def-error-budget" title="誤差予算">
最終結果の誤差 $|\hat y-y|$ を、原因ごとの項の和で上から抑えた表を**誤差予算**という。項は次の 4 種類に分ける。

- **打ち切り誤差**: 極限・級数・積分・微分方程式を有限の手続きで置き換えたことによる誤差。
- **丸め誤差**: 有限の精度で演算したことによる誤差。
- **データの誤差**: 入力（測定値、前段の計算結果、10 進の定数を 2 進に直したもの）が真の値と違うことによる誤差。
- **統計誤差**: 乱数を使う推定の標本による揺らぎ。
</Definition>

はじめの 3 つは決定的な上界で抑えられます。統計誤差は確率的な量で、標準誤差（[再現できる数値実験](/computer-science/research-computing/reproducible-experiments)の <Ref to="computer-science/research-computing/reproducible-experiments#prop-seed-averaging" />）で見積もり、表では別の行にします。

### 2.2. 3 項への分解

真の問題を $y=f(x)$、打ち切りを入れた（刻み $h$ などで離散化した）問題を $f_h$、実際に手にしている入力を $\hat x$、計算機が返した値を $\hat y$ とします。

<Proposition id="prop-error-budget" title="誤差の 3 項への分解">
$I\subset\mathbb{R}$ を開区間、$f,f_h:I\to\mathbb{R}$、$x,\hat x\in I$、$\hat y\in\mathbb{R}$ とする。このとき
$$
|\hat y-f(x)|\le\underbrace{|\hat y-f_h(\hat x)|}_{\text{丸め}}+\underbrace{|f_h(\hat x)-f(\hat x)|}_{\text{打ち切り}}+\underbrace{|f(\hat x)-f(x)|}_{\text{データ}}
$$
である。さらに、各項は次のように評価できる。

1. （丸め）$f_h$ が $I$ 上微分可能で $|f_h'|\le L_h$ とする。$\hat y=f_h(\hat x+\Delta)$, $\hat x+\Delta\in I$, $|\Delta|\le\eta|\hat x|$ となる $\Delta$ がある（$f_h$ の計算の後退誤差が $\eta$ 以下）ならば、$|\hat y-f_h(\hat x)|\le L_h\,|\hat x|\,\eta$。
2. （打ち切り）定数 $C\ge0$, $q\gt0$ があってすべての $z\in I$ で $|f_h(z)-f(z)|\le Ch^q$ ならば、$|f_h(\hat x)-f(\hat x)|\le Ch^q$。
3. （データ）$f$ が $I$ 上微分可能で $|f'|\le L$、$|\hat x-x|\le\varepsilon|x|$ ならば、$|f(\hat x)-f(x)|\le L\,|x|\,\varepsilon$。
</Proposition>

<Proof of="prop-error-budget">
$\hat y-f(x)=(\hat y-f_h(\hat x))+(f_h(\hat x)-f(\hat x))+(f(\hat x)-f(x))$ に三角不等式を使えば最初の不等式です。

1. $\Delta=0$ なら左辺は $0$ です。$\Delta\ne0$ なら、$\hat x$ と $\hat x+\Delta$ を端点とする閉区間は $I$ に含まれ、$f_h$ はその上で微分可能（したがって連続）なので、平均値の定理（<Ref to="mathematics/calculus/mean-value-and-taylor#thm-mvt" />）から、両者の間のある $c$ について $f_h(\hat x+\Delta)-f_h(\hat x)=f_h'(c)\Delta$ です。よって $|\hat y-f_h(\hat x)|\le L_h|\Delta|\le L_h|\hat x|\eta$ です。
2. 仮定を $z=\hat x$ に適用します。
3. 1 と同じく、$\hat x=x$ なら $0$、そうでなければ平均値の定理から $|f(\hat x)-f(x)|=|f'(c)||\hat x-x|\le L|x|\varepsilon$ です。
</Proof>

1 と 3 の右辺を $|f(x)|$ で割ると、$L|x|/|f(x)|$ という量が現れます。これは相対条件数 $\kappa(x)=|xf'(x)/f(x)|$（<Ref to="mathematics/numerical-analysis/conditioning-and-stability#prop-condition-differentiable" />）を、$x$ の近くでの上限に置き換えたものです。つまり、丸めの項は「後退誤差 × 条件数」、データの項は「入力の相対誤差 × 条件数」、打ち切りの項は離散化の次数 $q$ で決まります。3 つは**別々の手段**で小さくします。丸めは作業精度を上げる、打ち切りは $h$ を小さくするか次数の高い公式に替える、データは入力を作り直す、です。作業精度をいくら上げても、打ち切りとデータの項は変わりません。

次の例は、この分解を実際の数に当てはめたものです。$\zeta(3)$ の近似では、打ち切り（途中で項を数えるのをやめたこと）と丸め（足し算ごとの誤差）のどちらが大きいかを先に見積もり、足す順序によって丸めの項が大きく変わることを確かめます。途中の丸めの回数の数え方は細かいので、結論の表（予算と実際の比較）だけを先に見てから戻っても構いません。

<Example id="ex-zeta3-budget" title="ζ(3) の部分和の誤差予算">
$\zeta(3)=\sum_{k\ge1}k^{-3}$ を、$N=10^6$ 項の部分和 $S_N$ を倍精度で足して近似します。予算を先に立てます。

- **打ち切り**: $f(x)=x^{-3}$ は $[1,\infty)$ で正かつ単調減少なので、積分判定法（<Ref to="mathematics/calculus/series-and-convergence#thm-integral-test" /> の 2）から $\dfrac{1}{2(N+1)^2}\le\zeta(3)-S_N\le\dfrac{1}{2N^2}=5\times10^{-13}$。
- **丸め**: $k^2\le10^{12}\lt2^{53}$ なので `k*k` は厳密です。$k\le208063$ では $k^3\lt2^{53}$ で `k*k*k` も厳密、各項の丸めは除算の 1 回です。それより大きい $k$ では $k^3$ が $2^{53}$ を超えるので、乗算でもう 1 回丸めが入りえます（各項の丸めは高々 2 回）。大きい項から逐次に足すと、第 $k$ 項が通る加算は $k\le2$ で $N-1$ 回、$k\ge3$ で $N-k+1$ 回なので、各項に掛かる誤差因子は $k\le2$ で $N$ 個、$k\ge3$ で高々 $N-k+3\le N$ 個です。誤差因子の積の補題（<Ref to="mathematics/computer-assisted-proofs/floating-point-interval-arithmetic#lem-gamma" />。除数の側の丸めは $(1+\delta)^{-1}$ の形で入ります）から誤差は $\gamma_NS_N\approx1.3\times10^{-10}$ 以下です（逐次和の誤差 <Ref to="mathematics/numerical-analysis/floating-point-arithmetic#thm-recursive-sum" /> と同じ議論）。小さい項から足す順序では、最初に足す項が加算 $N-1$ 回と丸め 2 回を受けるので因子は $N+1$ 個以下、上界は $\gamma_{N+1}S_N$ で、数値は同じ $1.3\times10^{-10}$ です。
- **データ**: 入力は整数 $k$ だけなので $0$ です。

```python
# Python 3.9.6 / NumPy 2.0.2 / mpmath 1.3.0（arm64）
import math
import numpy as np
import mpmath as mp

N = 10**6
u = 2.0 ** -53
k = np.arange(1, N + 1, dtype=np.float64)
terms = 1.0 / (k * k * k)             # k <= 208063 では k^3 < 2^53 で厳密。各項の丸めは高々 2 回
S_seq = float(np.cumsum(terms)[-1])   # 大きい項から逐次に足す
S_rev = float(np.cumsum(terms[::-1])[-1])   # 小さい項から逐次に足す
S_fsum = math.fsum(terms)

mp.mp.dps = 40
zeta3 = mp.zeta(3)
S_exact = mp.fsum(1 / mp.mpf(j) ** 3 for j in range(1, N + 1))   # 部分和（40 桁）
trunc = zeta3 - S_exact
gamma = N * u / (1 - N * u)
print(f"打ち切り  真値 {mp.nstr(trunc, 6)}  上界 1/(2N^2) = {1 / (2 * N**2):.6e}"
      f"  下界 1/(2(N+1)^2) = {1 / (2 * (N + 1)**2):.6e}")
print(f"丸め      上界 γ_N·S = {gamma * S_fsum:.3e}")
for name, S in (("大→小", S_seq), ("小→大", S_rev), ("fsum", S_fsum)):
    print(f"  {name}: 丸め誤差 {mp.nstr(abs(mp.mpf(S) - S_exact), 3)}"
          f"  全誤差 {mp.nstr(abs(mp.mpf(S) - zeta3), 6)}")
print(f"補正つき  fsum + 1/(2N^2) の誤差 {mp.nstr(abs(mp.mpf(S_fsum) + mp.mpf(1) / (2 * N**2) - zeta3), 3)}")
```

```text
打ち切り  真値 4.99999e-13  上界 1/(2N^2) = 5.000000e-13  下界 1/(2(N+1)^2) = 4.999990e-13
丸め      上界 γ_N·S = 1.335e-10
  大→小: 丸め誤差 8.77e-12  全誤差 9.2733e-12
  小→大: 丸め誤差 9.37e-17  全誤差 5.00093e-13
  fsum: 丸め誤差 9.37e-17  全誤差 5.00093e-13
補正つき  fsum + 1/(2N^2) の誤差 9.32e-17
```

| 項 | 予算（上界） | 実際 |
|---|---|---|
| 打ち切り | $5.0\times10^{-13}$ | $5.0\times10^{-13}$ |
| 丸め（逐次和） | $1.3\times10^{-10}$ | $8.8\times10^{-12}$（大→小）、$9.4\times10^{-17}$（小→大） |
| データ | $0$ | $0$ |

**検算.** (a) 既知値: $\zeta(3)$ と部分和の参照値は mpmath の 40 桁です。打ち切りの真値は上界と下界の間にあります（7 桁目まで書くと $4.999995\times10^{-13}$ です。mpmath 1.4.1 では `mp.nstr` の丸め方が違い、同じ値が `5.0e-13` と表示されます）。(b) 別手法: 足す順序を変えた 3 通りで、丸め誤差はどれも予算 $1.3\times10^{-10}$ の内側です。(c) 予算の読み方: 予算の合計は丸めの項が支配していて、保証できるのは $10^{-10}$ までです。実際、大きい項から足すと丸め誤差 $8.8\times10^{-12}$ が打ち切り誤差を上回りました。「$N$ を増やせば精度が上がる」と考えて $N$ だけを増やしても、この足し方では丸めの項が増えて改善しません。丸めの項を減らすには、正しく丸めた和 `math.fsum` を使います。その誤差は、各項の丸め（除算の分が $uS_N$ 以下。$k\gt208063$ の項の乗算の分は、それらの項の和が $10^{-11}$ 程度なので無視できます）と和の丸め（$u|S_N|$ 程度）の合計で約 $2uS_N\approx2.7\times10^{-16}$ 以下です。すると打ち切りが支配的になり、それは上界と下界にはさまれた値 $1/(2N^2)$ を足して補正できます。補正後の誤差は $9\times10^{-17}$ でした。
</Example>

この例の教訓は、**どの項が支配的かによって打つ手が違う**ことです。予算を立てずに「倍精度だから 16 桁」と考えると、打ち切りの $5\times10^{-13}$ を見落とし、「項数を増やせばよい」と考えると丸めの項を見落とします。

## 3. 必要な作業桁数

### 3.1. 条件数の桁だけ失う

丸めの項について、作業精度をどれだけにすればよいかを定理にします。既存の「前進誤差 ≤ 条件数 × 後退誤差」（<Ref to="mathematics/numerical-analysis/conditioning-and-stability#thm-forward-backward" />）を、桁数の形に直したものです。

次の定理は仮定が長いですが、言いたいことは 1 つです。後退安定な算法（後退誤差が $cu$ 以下の算法）で計算した結果の相対誤差は、ほぼ「条件数 × $cu$」で、条件数が $10^k$ なら $k$ 桁を失う、ということです。1 は入力の近くで線形近似できる場合、2 は近似を使わず、入力の近くでの変化率の上限で評価する場合、3 は桁数の形に言い換えたものです。

<Theorem id="thm-working-precision" title="必要な作業桁数">
$U\subset\mathbb{R}^n$ を開集合、$f:U\to\mathbb{R}^m$、$\boldsymbol{x}\in U$, $\boldsymbol{x}\ne\boldsymbol{0}$, $f(\boldsymbol{x})\ne\boldsymbol{0}$ とする。単位丸め $u$ の浮動小数点演算で $f(\boldsymbol{x})$ を計算する算法が、ある定数 $c\gt0$ について次を満たすとする（後退誤差が $cu$ 以下。<Ref to="mathematics/numerical-analysis/conditioning-and-stability#def-backward-stable" />）。

> 計算値 $\hat{\boldsymbol{y}}_u$ は、$\boldsymbol{x}+\Delta\boldsymbol{x}_u\in U$, $\|\Delta\boldsymbol{x}_u\|\le cu\|\boldsymbol{x}\|$ を満たすある $\Delta\boldsymbol{x}_u$ について $\hat{\boldsymbol{y}}_u=f(\boldsymbol{x}+\Delta\boldsymbol{x}_u)$ と書ける。

相対前進誤差を $E_u=\|\hat{\boldsymbol{y}}_u-f(\boldsymbol{x})\|/\|f(\boldsymbol{x})\|$ とおく。

1. $f$ が $\boldsymbol{x}$ で全微分可能ならば、相対条件数 $\kappa(\boldsymbol{x})$（<Ref to="mathematics/numerical-analysis/conditioning-and-stability#def-condition-number" />）について
   $$
   E_u\le c\,\kappa(\boldsymbol{x})\,u+r(u),\qquad\lim_{u\to+0}\frac{r(u)}{u}=0
   $$
   を満たす $r(u)\ge0$ がある。
2. $u_0\gt0$ とし、閉球 $\bar B=\{\boldsymbol{z}:\|\boldsymbol{z}-\boldsymbol{x}\|\le cu_0\|\boldsymbol{x}\|\}$ が $U$ に含まれ、$f$ が $\bar B$ の各点で全微分可能で、$K=\sup_{\boldsymbol{z}\in\bar B}\|J_f(\boldsymbol{z})\|\,\|\boldsymbol{x}\|/\|f(\boldsymbol{x})\|\lt\infty$ とする。このとき $u\le u_0$ ならば $E_u\le c\,K\,u$。
3. 2 の状況で、作業桁数 $D=-\log_{10}u$ が
   $$
   D\ge d+\log_{10}K+\log_{10}c
   $$
   を満たせば、$E_u\le10^{-d}$ である。すなわち、目標の $d$ 桁に、条件数の桁 $\log_{10}K$ と算法の定数の桁 $\log_{10}c$ を足した作業桁数で足りる。
</Theorem>

<Proof of="thm-working-precision">
**1.** <Ref to="mathematics/numerical-analysis/conditioning-and-stability#thm-forward-backward" /> の 1 から、$\boldsymbol{x}$ と $f$ だけで決まる関数 $R\ge0$ があって、$\Delta\boldsymbol{x}_u\ne\boldsymbol{0}$ のとき
$$
E_u\le\kappa(\boldsymbol{x})\frac{\|\Delta\boldsymbol{x}_u\|}{\|\boldsymbol{x}\|}+R(\Delta\boldsymbol{x}_u)\le c\,\kappa(\boldsymbol{x})\,u+R(\Delta\boldsymbol{x}_u),\qquad\lim_{\Delta\boldsymbol{x}\to\boldsymbol{0}}\frac{R(\Delta\boldsymbol{x})}{\|\Delta\boldsymbol{x}\|}=0
$$
です。$r(u)=R(\Delta\boldsymbol{x}_u)$（$\Delta\boldsymbol{x}_u=\boldsymbol{0}$ のときは $E_u=0$ なので $r(u)=0$）とおきます。$\Delta\boldsymbol{x}_u\ne\boldsymbol{0}$ のとき
$$
\frac{r(u)}{u}=\frac{R(\Delta\boldsymbol{x}_u)}{\|\Delta\boldsymbol{x}_u\|}\cdot\frac{\|\Delta\boldsymbol{x}_u\|}{u}\le\frac{R(\Delta\boldsymbol{x}_u)}{\|\Delta\boldsymbol{x}_u\|}\cdot c\|\boldsymbol{x}\|
$$
で、$u\to+0$ のとき $\|\Delta\boldsymbol{x}_u\|\le cu\|\boldsymbol{x}\|\to0$ なので右辺は $0$ に収束します。

**2.** $u\le u_0$ のとき $\|\Delta\boldsymbol{x}_u\|\le cu_0\|\boldsymbol{x}\|$ なので、線分 $\{\boldsymbol{x}+t\Delta\boldsymbol{x}_u:0\le t\le1\}$ は $\bar B$ に含まれます（球は凸）。<Ref to="mathematics/numerical-analysis/conditioning-and-stability#thm-forward-backward" /> の 2 から、$E_u$ は「線分上での $\|J_f(\boldsymbol{z})\|\|\boldsymbol{x}\|/\|f(\boldsymbol{x})\|$ の上限」と $\|\Delta\boldsymbol{x}_u\|/\|\boldsymbol{x}\|\le cu$ の積以下で、線分上の上限は $K$ 以下です。

**3.** $E_u\le cKu$ の右辺の常用対数は $\log_{10}c+\log_{10}K-D$ で、仮定からこれは $-d$ 以下です。
</Proof>

3 が実務で使う形です。binary64 の作業桁数は $15.95$ なので、$\kappa\approx10^{10}$ の問題では、後退安定な算法を使っても得られるのは 5〜6 桁です。10 桁欲しければ、作業桁数を $10+10+\log_{10}c$、つまり 21 桁以上にします。$f$ が 2 回連続微分可能なら 1 の $r(u)$ は $O(u^2)$ で、$\kappa u\ll1$ の範囲では 1 と 2 の差は無視できます。

連立 1 次方程式 $A\boldsymbol{x}=\boldsymbol{b}$ では、後退誤差 $\eta$ の解の相対誤差は $2\kappa(A)\eta/(1-\kappa(A)\eta)$ 以下です（<Ref to="mathematics/numerical-linear-algebra/norms-and-conditioning#cor-forward-error-bound" />。$\kappa(A)=\|A\|\|A^{-1}\|$ は <Ref to="mathematics/numerical-linear-algebra/norms-and-conditioning#def-condition-number" />）。$\kappa(A)\eta\le1/2$ なら $4\kappa(A)\eta$ 以下なので、桁数の勘定は同じで、「$\log_{10}\kappa(A)$ 桁を失う」と読めます。

<Example id="ex-hilbert-digits" title="Hilbert 行列：得られる桁数は作業桁数 − log₁₀κ">
Hilbert 行列 $H_{ij}=1/(i+j-1)$ は条件数が $n$ とともに急に大きくなることで知られています。binary64 に丸めた $H$ と $\boldsymbol{b}=H\boldsymbol{1}$ を問題として固定し、倍精度（`np.linalg.solve`）と mpmath（`mp.lu_solve`）の 30 桁・50 桁で解いて、80 桁の参照解と比べます。表の各組は「得られた有効桁, 予測（作業桁数 $-\log_{10}\kappa_\infty$）」です。

```python
# Python 3.9.6 / NumPy 2.0.2 / SciPy 1.13.1 / mpmath 1.3.0（arm64）
import math
import numpy as np
import mpmath as mp
from scipy.linalg import hilbert

def digits(x, ref):                   # 相対誤差（∞ ノルム）の -log10
    err = max(abs(a - b) for a, b in zip(x, ref)) / max(abs(b) for b in ref)
    return float(-mp.log10(err))

def work_digits(dps):                 # mp.lu_solve は内部で精度を 10 ビット上げる
    mp.mp.dps = dps
    return (mp.mp.prec + 10) * math.log10(2)

D64 = 53 * math.log10(2)              # binary64 の作業桁数 15.95
print(" n log10κ | float64 予測 | dps=30 予測 | dps=50 予測")
for n in range(6, 15, 2):
    H = hilbert(n)                    # binary64 に丸めた Hilbert 行列（これを問題とする）
    b = H @ np.ones(n)
    mp.mp.dps = 80                    # 参照解: 同じ H, b を 80 桁で解く
    Hm, bm = mp.matrix(H.tolist()), mp.matrix(b.tolist())
    ref = mp.lu_solve(Hm, bm)
    lk = float(mp.log10(mp.norm(Hm, mp.inf) * mp.norm(mp.inverse(Hm), mp.inf)))
    row = [(digits(np.linalg.solve(H, b).tolist(), ref), D64 - lk)]
    for dps in (30, 50):
        wd = work_digits(dps)
        x = mp.lu_solve(mp.matrix(H.tolist()), mp.matrix(b.tolist()))
        mp.mp.dps = 80
        row.append((digits(x, ref), wd - lk))
    print(f"{n:>2} {lk:6.1f} | " + " | ".join(f"{r:6.1f} {p:5.1f}" for r, p in row))
```

```text
 n log10κ | float64 予測 | dps=30 予測 | dps=50 予測
 6    6.6 |    9.8   9.3 |   28.3  27.4 |   47.8  47.2
 8    9.6 |    6.4   6.3 |   25.3  24.4 |   45.4  44.3
10   12.5 |    3.6   3.4 |   22.8  21.5 |   41.8  41.3
12   15.6 |    0.8   0.4 |   19.0  18.5 |   39.0  38.3
14   16.7 |   -0.4  -0.8 |   19.1  17.3 |   38.6  37.2
```

**検算.** (a) 参照: 参照解は同じ $H$, $\boldsymbol{b}$ を 80 桁で解いたもので、$\kappa\le10^{17}$ なので <Ref to="thm-working-precision" /> から 60 桁以上正しく、表の比較には十分です。(b) 別精度: 3 つの作業精度のすべての行で、得られた桁数は予測と 2 桁以内で一致し、予測より少し良い側にあります（予測は最悪の場合の評価なので、良い側にずれるのは定理と矛盾しません）。(c) 次数: $n$ を 2 増やすごとに $\log_{10}\kappa$ は約 3 増え、得られる桁数はどの精度でも約 3 ずつ減ります。倍精度は $n=12$ で 1 桁を割り、$n=14$ では有効桁が負、つまり相対誤差が 1 を超えています。mpmath の予測には、`mp.lu_solve` が内部で作業精度を 10 ビット（約 3 桁）上げていること（ソースコードで確認）を入れてあります。$n=14$ の $\log_{10}\kappa$ が $16.7$ にとどまるのは、binary64 に丸めた行列を問題にしているからで、真の Hilbert 行列の条件数ではありません。
</Example>

### 3.2. 作業桁数を決める手順

<Ref to="thm-working-precision" /> を手順に直すと次のようになります。

1. **条件数を見積もる。** 解析的にわかるならその値を使います。わからなければ、入力を相対 $10^{-8}$ 程度ずらして出力の相対変化を見る（比が条件数の目安）か、行列なら条件数の推定ルーチンを使います。
2. **桁数を決める。** 目標の $d$ 桁に $\log_{10}\kappa$ を足し、算法の定数と見積もりの粗さのために数桁（3〜5 桁）の余裕を足します。mpmath なら `mp.dps` にその値を設定します。
3. **倍の桁で再計算して確かめる。** 作業桁数を 2 倍にして同じ計算を繰り返し、2 つの結果が目標の桁数まで一致することを見ます（7 節。これは丸め誤差の検査で、打ち切り誤差の検査にはなりません）。

## 4. 桁落ちで失う桁数

条件数が大きくなる最もありふれた原因は、近い数の引き算と、正負が打ち消し合う和です。失う桁数を式にしておきます。

次の補題は、$e^x$ が級数で書けるという基本の事実です。交代級数 $e^{-x}=\sum(-x)^k/k!$ の項の絶対値を足すと $e^x$ になり、答えの $e^{-x}$ より桁違いに大きい値を途中で扱うことになります。この差が、後の命題 3 で失う桁数の原因です。

<Lemma id="lem-exp-series" title="指数級数">
すべての実数 $x$ について $\displaystyle\sum_{k=0}^{\infty}\frac{x^k}{k!}=e^x$ である。
</Lemma>

<Proof of="lem-exp-series">
$x=0$ なら両辺は $1$ です。$x\ne0$ とします。$\exp$ は $\mathbb{R}$ 上で何回でも微分可能で、導関数はすべて $\exp$ です（<Ref to="mathematics/calculus/derivatives#thm-elementary-derivatives" /> の 4）。テイラーの定理（<Ref to="mathematics/calculus/mean-value-and-taylor#thm-taylor" />）を $a=0$ で使うと、各 $n\ge0$ について、$0$ と $x$ の間のある $c_n$ で
$$
e^x-\sum_{k=0}^n\frac{x^k}{k!}=\frac{e^{c_n}}{(n+1)!}x^{n+1}
$$
です。$\exp'=\exp\gt0$ なので $\exp$ は狭義単調増加で（<Ref to="mathematics/calculus/mean-value-and-taylor#cor-monotonicity" /> の 3）、$c_n\lt|x|$ から $e^{c_n}\le e^{|x|}$ です。$a_n=|x|^{n+1}/(n+1)!$ とおくと、右辺の絶対値は $e^{|x|}a_n$ 以下です。$N\ge2|x|$ となる整数 $N$ をとると、$n\ge N$ で $a_{n+1}/a_n=|x|/(n+2)\lt1/2$ なので $a_n\le a_N2^{-(n-N)}\to0$ です。よって部分和は $e^x$ に収束します。
</Proof>

<Proposition id="prop-digits-lost" title="差と和で失う桁数">
1. （差）$a\ne b$ を実数、$\hat a=a(1+\alpha)$, $\hat b=b(1+\beta)$, $|\alpha|,|\beta|\le\varepsilon$ とする。$\hat a-\hat b$ の相対誤差は $\varepsilon\cdot\dfrac{|a|+|b|}{|a-b|}$ 以下で、$a,b\gt0$ のときこの値をとりうる。したがって、有効桁は入力に比べて高々
   $$
   \log_{10}\frac{|a|+|b|}{|a-b|}
   $$
   桁減る。
2. （和）$x_1,\ldots,x_n$ を実数、$s=\sum_ix_i\ne0$、$\hat x_i=x_i(1+\varepsilon_i)$, $|\varepsilon_i|\le\varepsilon$ とする。$\sum_i\hat x_i$ の相対誤差は $\varepsilon\cdot\dfrac{\sum_i|x_i|}{|s|}$ 以下で、$\varepsilon_i=\varepsilon\,\mathrm{sgn}(x_i)$ のとき等号が成り立つ。失う桁数は高々 $\log_{10}\bigl(\sum_i|x_i|/|s|\bigr)$ である。
3. （交代級数）$x\gt0$ とし、$e^{-x}$ を級数 $\sum_{k\ge0}t_k$, $t_k=(-x)^k/k!$ の部分和 $S_N=\sum_{k=0}^Nt_k$ で計算する。2 の倍率 $A_N=\sum_{k=0}^N|t_k|/|S_N|$ は $N\to\infty$ で $e^{2x}$ に収束する。失う桁数は
   $$
   \log_{10}e^{2x}=\frac{2x}{\ln10}\approx0.869\,x
   $$
   に近づく。
</Proposition>

<Proof of="prop-digits-lost">
1. 桁落ちの相対誤差の評価（<Ref to="mathematics/numerical-analysis/floating-point-arithmetic#prop-cancellation" />）そのものです。相対誤差が $\varepsilon$ から $\varepsilon M$（$M=(|a|+|b|)/|a-b|\ge1$）に増えると、有効桁 $-\log_{10}(\text{相対誤差})$ は $\log_{10}M$ だけ減ります。
2. $\sum_i\hat x_i-s=\sum_ix_i\varepsilon_i$ なので、三角不等式から $|\sum_i\hat x_i-s|\le\varepsilon\sum_i|x_i|$ です。$\varepsilon_i=\varepsilon\,\mathrm{sgn}(x_i)$ なら $x_i\varepsilon_i=\varepsilon|x_i|$ で、和は $\varepsilon\sum_i|x_i|$ となり等号です。$|s|$ で割れば主張です。
3. $|t_k|=x^k/k!$ なので、<Ref to="lem-exp-series" /> から $\sum_{k=0}^N|t_k|\to e^x$、また $x$ を $-x$ に置き換えて $S_N\to e^{-x}$ です。$e^{-x}\ne0$ なので $A_N\to e^x/e^{-x}=e^{2x}$ です。$\log_{10}e^{2x}=2x\log_{10}e=2x/\ln10$ で、$2/\ln10=0.8685\ldots$ です。
</Proof>

3 の直感はこうです。$e^{-x}$ は $x=20$ で約 $2\times10^{-9}$ と小さいのに、級数の途中には $x^k/k!$ という大きな項（$k\approx x$ で $10^7$ 台）が現れ、それらが打ち消し合って小さな答えになります。足す項の大きさの合計は約 $e^{x}$、答えは $e^{-x}$ なので、比が $e^{2x}$ になります。

2 の倍率 $\sum|x_i|/|\sum x_i|$ は、和という問題の成分ごとの条件数です。3 は「答えが小さいのに、途中の項が大きい」計算の典型で、$x=20$ では $17.4$ 桁を失います。倍精度の作業桁数は $15.95$ なので、1 桁も残りません。

<Example id="ex-exp-alternating" title="交代級数で e⁻ˣ を計算する">
$x=20,30,40$ について、交代級数を倍精度と mpmath（`dps` が 30 と 60）で足し、有効桁を調べます。参照値は 100 桁の `mp.exp(-x)` です。

```python
# Python 3.9.6 / mpmath 1.3.0（arm64）
import math
import mpmath as mp

def exp_neg_series(x, one):           # Σ (-x)^k / k! を項が効かなくなるまで足す
    s = term = one
    k = 0
    while True:
        k += 1
        term = term * (-x) / k
        new = s + term
        if new == s and k > x:
            return s
        s = new

def digits(approx, exact):            # 正しい 10 進桁数（相対誤差の -log10）
    err = abs((mp.mpf(approx) - exact) / exact)
    return float(-mp.log10(err)) if err else float("inf")

print(" x  失う桁   float64   dps=30   dps=60")
for x in (20, 30, 40):
    mp.mp.dps = 100
    exact = mp.exp(-x)
    row = [digits(exp_neg_series(float(x), 1.0), exact)]
    for dps in (30, 60):
        mp.mp.dps = dps
        val = exp_neg_series(mp.mpf(x), mp.mpf(1))
        mp.mp.dps = 100
        row.append(digits(val, exact))
    lost = 2 * x / math.log(10)
    print(f"{x:>2} {lost:7.1f} " + " ".join(f"{d:8.1f}" for d in row))
mp.mp.dps = 15
print(exp_neg_series(20.0, 1.0), math.exp(-20), 1 / exp_neg_series(-20.0, 1.0))
```

```text
 x  失う桁   float64   dps=30   dps=60
20    17.4     -0.2     14.7     45.0
30    26.1     -8.5      6.0     36.6
40    34.7    -17.9     -2.6     27.5
5.621884472130418e-09 2.061153622438558e-09 2.0611536224385583e-09
```

**検算.** (a) 既知値: 参照値は 100 桁の $e^{-x}$ で、最後の行の `math.exp(-20)` $=2.0611\ldots\times10^{-9}$ とも一致します。(b) 命題との照合: mpmath の作業桁数は `dps=30` で 103 ビット（$31.0$ 桁）、`dps=60` で 203 ビット（$61.1$ 桁）です。作業桁数から <Ref to="prop-digits-lost" /> の「失う桁」を引いた予測は、倍精度で $-1.4,-10.1,-18.8$、`dps=30` で $13.6,4.9,-3.7$、`dps=60` で $43.7,35.0,26.4$ で、得られた桁数はどれもこれより 1〜2 桁良い側にあります（負の桁数は、相対誤差が 1 を超えていることを表します）。命題の倍率は最悪の場合の上界なので、良い側にずれるのは矛盾ではありません。(c) 別の式: 最後の行の 3 つ目は $e^{-20}=1/e^{20}$ として、すべての項が正の級数で $e^{20}$ を計算してから逆数をとったものです。打ち消しがないので、倍精度のままで 16 桁合います。交代級数の値 $5.6\times10^{-9}$ は 1 桁目から違います。
</Example>

精度を上げるのは**最後の手段**です。この例では、式を $1/e^{x}$ に書き換えるだけで倍精度で足ります。書き換えが見つからないとき、あるいは書き換えた式の正しさを確かめたいときに、<Ref to="thm-working-precision" /> と <Ref to="prop-digits-lost" /> で必要な桁数を見積もって任意精度を使います。

## 5. 精度の道具

### 5.1. 何が使えるか

<Definition id="def-arbitrary-precision" title="任意精度浮動小数と厳密な有理数">
精度 $p$ を実行時に選べる浮動小数点演算を**任意精度浮動小数**という。誤差モデルは固定精度と同じで、単位丸めが $u=2^{-p}$ になる。Python の mpmath では、`mp.prec` が 2 進の精度 $p$、`mp.dps` が対応する 10 進の桁数で、一方を設定すると他方も決まる（`mp.dps = 30` で `mp.prec` は $103$）。

分子と分母を任意の桁数の整数でもつ有理数（Python の `fractions.Fraction`）の四則は丸めを伴わない。これを**厳密な有理数演算**という。
</Definition>

<Definition id="def-mixed-precision" title="混合精度">
1 つの計算の中で、部分ごとに異なる精度を使うことを**混合精度**という。費用の大半を占める部分を低い精度で、結果の精度を決める部分を高い精度で行う。
</Definition>

混合精度の代表例は、連立 1 次方程式の**反復改良**です。LU 分解と補正方程式の求解を低い精度で行い、残差 $\boldsymbol{b}-A\boldsymbol{x}_k$ だけを高い精度で計算すると、誤差は 1 回ごとに縮みます（<Ref to="mathematics/numerical-linear-algebra/lu-and-pivoting#thm-iterative-refinement" />）。$O(n^3)$ の分解は 1 回で済み、高い精度が要るのは $O(n^2)$ の残差の計算だけです。

Python で使える数の型を、精度の観点で並べます。速度の欄は桁の目安で、計測値ではありません。

| 型 | 精度 | 速度の目安 | 配列演算 | 包含の保証 |
|---|---|---|---|---|
| `np.float32` | 24 ビット（7.2 桁） | ハードウェア。倍精度と同程度かやや速い | NumPy | なし |
| `float`, `np.float64` | 53 ビット（16.0 桁） | ハードウェア | NumPy | なし |
| `np.longdouble` | **環境依存**。x86-64 の Linux では仮数 64 ビット、手元の arm64 の macOS では 53 ビット（`float64` と同じ） | 環境依存 | NumPy | なし |
| mpmath の `mpf` | 任意 | ソフトウェア実装。倍精度より 2〜3 桁遅い | `mp.matrix`（要素ごとに Python の演算） | なし |
| mpmath の `iv`（区間） | 任意 | `mpf` と同程度 | `iv.matrix` | **あり** |
| python-flint の `arb`（球） | 任意 | C 実装。mpmath より速い | 行列型あり | **あり** |
| `fractions.Fraction` | 厳密（四則のみ） | 分母の桁数が増えると急に遅くなる | なし | 誤差なし |

`np.longdouble` は、同じコードが環境によって別の精度で動くので、再現性（[再現できる数値実験](/computer-science/research-computing/reproducible-experiments)）の観点から、結果を左右する場所には使わないのがよいと思います。「倍精度より少し多く」が必要なら、桁数を明示できる任意精度を使います。

### 5.2. 高精度計算を台無しにする誤り

任意精度の計算で最も多い誤りは、**入力や定数を倍精度で作ってから渡す**ことです。`mp.dps = 30` のもとで、手元の mpmath 1.3.0 は次を返します。

- `mp.mpf(0.1)` は `0.100000000000000005551115123126`。Python が `0.1` を binary64 に丸めたあとの値を 30 桁で表したもので、17 桁目から $1/10$ と違います。
- `mp.mpf("0.1")` は `0.1`。文字列から作れば、10 進の $0.1$ を作業精度で丸めた数になります。
- `mp.mpf(math.pi) - mp.pi` は $-1.2246\times10^{-16}$。`math.pi` は倍精度の定数です。高精度の $\pi$ は `mp.pi` を使います。

<Ref to="prop-error-budget" /> の言葉では、これは**データの誤差**です。作業精度を 30 桁にしても 1000 桁にしても、入力の相対誤差 $\varepsilon\approx10^{-16}$ は変わらず、結果の誤差は $\kappa\cdot10^{-16}$ 程度のままです。高精度の数は、文字列（`mp.mpf("0.1")`）、整数の比（`mp.mpf(1) / 10`）、または `Fraction` から作ります。式の途中に書いた `2.5 * x` のような Python の浮動小数リテラルも、$2.5$ のように 2 進で正確に表せる数なら無害ですが、`0.1 * x` は同じ誤りになります。

一方、入力がもともと倍精度の数として与えられている場合（測定値や、前段の倍精度計算の出力）は、`mp.mpf(x)` で正確にその値を取り込むのが正しい扱いです。4 節までの例で「binary64 に丸めた行列を問題とする」としたのはこの立場です。どちらの立場で計算しているかを、コードのコメントに書いてください。

### 5.3. 精度を倍々にする

高精度の計算は、桁数が多いほど 1 回の演算が高くつきます。2 次収束する反復では、はじめから目標の精度で計算する必要はありません。

この命題の発想は、Newton 法のような 2 次収束する反復は、1 回で正しい桁数がほぼ倍になる、という点にあります。たとえば今 50 桁正しければ、次の反復のあとは約 100 桁正しくなります。だから反復の途中は、これから得る桁数より少し多い程度の作業精度で足り、はじめから目標の桁数で計算する必要はありません。1 は誤差が 2 乗ずつ減ること（記号 $\tau_j$ はその上界）、2 は費用の合計が最後の 1 回の 2 倍未満で済むことを述べています。

<Proposition id="prop-precision-doubling" title="精度を倍にしながらの反復">
1. （精度）$J=[x^*-\delta,x^*+\delta]$、$\varphi:J\to\mathbb{R}$ とし、定数 $K\gt0$ についてすべての $x\in J$ で $|\varphi(x)-x^*|\le K|x-x^*|^2$ とする。$0\lt\tau_0\le\delta$, $2K\tau_0\lt1$ とし、$\tau_{j+1}=2K\tau_j^2$ と定める。$x_0\in J$, $|x_0-x^*|\le\tau_0$ から出発して
   $$
   x_{j+1}=\varphi(x_j)+\rho_j,\qquad|\rho_j|\le K\tau_j^2
   $$
   （$\rho_j$ は $j$ 回目の反復を有限の精度で行ったことによる誤差）で列を作ると、すべての $j$ で $x_j\in J$, $|x_j-x^*|\le\tau_j$ であり、$2K\tau_j=(2K\tau_0)^{2^j}$ である。
2. （費用）関数 $M:\{1,2,\ldots\}\to(0,\infty)$ が、すべての $q$ で $M(2q)\ge2M(q)$ を満たすとする。目標の精度を $p=2^Jp_0$（$p_0\ge1$, $J\ge0$ は整数）とし、精度 $q$ での反復 1 回の費用が $aM(q)$ 以下（$a\gt0$ は定数）とする。反復を精度 $p_0,2p_0,4p_0,\ldots,2^Jp_0=p$ で 1 回ずつ行うと、費用の合計は
   $$
   \sum_{i=0}^{J}aM\bigl(p/2^i\bigr)\le aM(p)\sum_{i=0}^{J}2^{-i}\lt2aM(p)
   $$
   である。すべての反復を精度 $p$ で行う場合の上界 $(J+1)aM(p)$ に対し、最後の 1 回の 2 倍未満で済む。
</Proposition>

<Proof of="prop-precision-doubling">
**1.** まず $2K\tau_{j+1}=(2K\tau_j)^2$ なので、帰納法で $2K\tau_j=(2K\tau_0)^{2^j}\lt1$ です。すると $\tau_{j+1}=(2K\tau_j)\tau_j\lt\tau_j$ で、$\tau_j\le\tau_0\le\delta$ です。$|x_j-x^*|\le\tau_j$ を $j$ の帰納法で示します。$j=0$ は仮定です。$j$ で成り立てば $x_j\in J$ なので
$$
|x_{j+1}-x^*|\le|\varphi(x_j)-x^*|+|\rho_j|\le K|x_j-x^*|^2+K\tau_j^2\le2K\tau_j^2=\tau_{j+1}
$$
で、$\tau_{j+1}\le\delta$ から $x_{j+1}\in J$ です。

**2.** $M(2q)\ge2M(q)$ を $q=p/2^i$（$1\le i\le J$）に繰り返し使うと $M(p)\ge2^iM(p/2^i)$、すなわち $M(p/2^i)\le2^{-i}M(p)$ です。$\sum_{i=0}^J2^{-i}=2-2^{-J}\lt2$ です。
</Proof>

Newton 法は 1 の仮定を満たします。$f$ が $C^2$ 級で $f(x^*)=0$, $f'(x^*)\ne0$ なら、$x^*$ の近くで $|x_{k+1}-x^*|\le K|x_k-x^*|^2$ です（<Ref to="mathematics/numerical-analysis/nonlinear-equations#thm-newton-local" />）。1 が言うのは、$j$ 回目の反復の丸め誤差は $K\tau_j^2=\tau_{j+1}/2$ まで許される、つまり**その反復のあとで達成できる精度と同じ程度の作業精度で足りる**ということです。達成できる精度は 1 回ごとに桁数が倍になるので、作業精度も倍々にすればよく、丸め誤差があっても誤差は $(2K\tau_0)^{2^j}$ の形で 2 乗ずつ減ります。2 の仮定 $M(2q)\ge2M(q)$ は、桁数を倍にすると費用が倍以上になることで、$p$ ビットの乗算の費用 $M(p)$ について、筆算（$p^2$ に比例）でも高速な乗算法でも成り立ちます。

<Example id="ex-newton-doubling" title="√2 を 3200 ビットまで">
$f(x)=x^2-2$ の Newton 法 $x\leftarrow(x+2/x)/2$ を、倍精度の $\sqrt2$ から出発し、作業精度を 100, 200, …, 3200 ビットと倍にしながら 1 回ずつ適用します。

```python
# Python 3.9.6 / mpmath 1.3.0（arm64）
import math
import mpmath as mp

target = 3200                          # 目標の精度（ビット）
precs = [100, 200, 400, 800, 1600, 3200]
mp.mp.prec = target + 64
exact = mp.sqrt(2)

def correct_bits(x):
    mp.mp.prec = target + 64
    return float(-mp.log(abs(x - exact) / exact, 2))

x = mp.mpf(math.sqrt(2))               # 初期値: binary64 の平方根
print(f"初期値      正しいビット数 {correct_bits(x):7.1f}")
for p in precs:                        # 精度を倍にしながら x ← (x + 2/x)/2
    mp.mp.prec = p
    x = (x + 2 / x) / 2
    print(f"prec={p:>4}   正しいビット数 {correct_bits(x):7.1f}")
ratio = sum(p * p for p in precs) / target**2     # 乗除算の費用を p^2 として数える
print(f"費用の比: 倍々 {ratio:.3f}  すべて目標精度 {len(precs):.3f}")
```

```text
初期値      正しいビット数    53.7
prec= 100   正しいビット数   100.1
prec= 200   正しいビット数   199.6
prec= 400   正しいビット数   401.0
prec= 800   正しいビット数   800.6
prec=1600   正しいビット数  1602.9
prec=3200   正しいビット数  3202.2
費用の比: 倍々 1.333  すべて目標精度 6.000
```

**検算.** (a) 既知値: 参照値は 3264 ビットの `mp.sqrt(2)` です。(b) 命題との照合: 各回のあとで正しいビット数はその回の作業精度にほぼ等しく、6 回で目標の 3200 ビットに達しました。すべてを 3200 ビットで行っても反復回数は 6 回で変わりません（正しいビット数は 1 回で高々倍にしかならないため）。(c) 費用: 乗除算の費用を $p^2$ とすると、倍々の場合の合計は最後の 1 回の $1+\tfrac14+\tfrac1{16}+\cdots\approx1.333$ 倍で、<Ref to="prop-precision-doubling" /> の上界 2 倍の内側です。すべてを目標精度で行うと 6 倍です。ここでの費用は演算量のモデルによる比で、時間の計測値ではありません。
</Example>

## 6. 区間演算と球算術：包含の保証

### 6.1. 保証があるものとないもの

浮動小数の計算は、誤差がどれだけあるかを結果に含めません。区間演算と球算術は、真値を含む集合（区間 $[\ell,r]$、あるいは中心と半径の組 $B(m,r)=[m-r,m+r]$。<Ref to="mathematics/computer-assisted-proofs/floating-point-interval-arithmetic#def-ball" />）を結果として返します。

次の命題の証明（4）では、記号 $\nabla$ と $\Delta$ を使います。$\nabla x$ は $x$ を下向きに丸めた値（$x$ 以下で最も近い浮動小数点数）、$\Delta x$ は上向きに丸めた値（$x$ 以上で最も近い浮動小数点数）を表します。球算術は、このように丸めの向きを意識して範囲を見積もります。

<Proposition id="prop-ball-vs-float" title="球算術の保証と限界">
球 $B(m,r)$ と演算 $\oplus,\otimes$ は <Ref to="mathematics/computer-assisted-proofs/floating-point-interval-arithmetic#def-ball" /> のものとし、$B_1\ominus B_2=B_1\oplus B(-m_2,r_2)$ と定める。

1. （包含）$m_1+m_2$, $m_1m_2$ が $0$ か正規化数の範囲にあるとき、$x\in B_1$, $y\in B_2$ ならば $x+y\in B_1\oplus B_2$, $xy\in B_1\otimes B_2$。
2. （依存性）$B=B(m,r)$ について $B\ominus B=B(0,r')$, $r'\ge2r$ である。$\{x-x:x\in B\}=\{0\}$ なのに、結果の半径は $2r$ 以上になる。
3. （精度を上げると縮む）$m_1,m_2$ を精度 $p_0$ の浮動小数点数とし、精度 $p\ge p_0$（単位丸め $u=2^{-p}$）で半径 $0$ の球 $B(m_1,0)$, $B(m_2,0)$ の和と積を計算する。半径を定義の右辺の $(1+\lambda)$ 倍以下（$\lambda\ge0$ は実装で決まる定数）にとるならば、オーバーフローもアンダーフローも起きないとき、和の球の半径は $(1+\lambda)(1+u)\,u\,|m_1+m_2|$ 以下、積の球の半径は $(1+\lambda)(1+u)\,u\,|m_1m_2|$ 以下である。
4. （2 つの精度の一致は保証にならない）$q\ge2$ を整数、$x=2^{-q}$ とする。精度 $p\lt q$ の浮動小数点演算（$x$ が表せる指数範囲をもつもの）で $(1+x)-1$ を計算すると、結果は $0$ である。したがって $p_1\lt p_2\lt q$ の 2 つの精度の結果は完全に一致する（ともに $0$）が、真値は $x\ne0$ で、相対誤差は $1$ である。同じ計算を球算術で行うと、結果の球の半径は $u$ 以上で、真値 $x$ を含む。
</Proposition>

<Proof of="prop-ball-vs-float">
**1.** 球算術の包含の定理（<Ref to="mathematics/computer-assisted-proofs/floating-point-interval-arithmetic#thm-ball-inclusion" />）そのものです。

**2.** $B\ominus B=B(m,r)\oplus B(-m,r)$ の中心は $\mathrm{fl}(m+(-m))=\mathrm{fl}(0)=0$、半径は定義から $r'\ge r+r+u\cdot|0|=2r$ です。

**3.** 精度 $p_0$ の浮動小数点数は、仮数の末尾に $0$ を補えば精度 $p$ の浮動小数点数でもあります。和の球の中心は $m=\mathrm{fl}(m_1+m_2)$ で、四則の標準モデル（<Ref to="mathematics/numerical-analysis/floating-point-arithmetic#cor-standard-model-ops" />）から $|m|\le(1+u)|m_1+m_2|$ です。定義の右辺は $0+0+u|m|$ なので、半径は $(1+\lambda)u|m|\le(1+\lambda)(1+u)u|m_1+m_2|$ 以下です。積も、中心 $m=\mathrm{fl}(m_1m_2)$, $|m|\le(1+u)|m_1m_2|$、右辺 $|m_1|\cdot0+|m_2|\cdot0+0\cdot0+u|m|$ から同じです。

**4.** $0\lt x=2^{-q}\le2^{-p-1}=u/2$ です。$\mathbb{F}\cap[1,2]$ は間隔 $2u$ の等間隔な点列（<Ref to="mathematics/computer-assisted-proofs/floating-point-interval-arithmetic#lem-spacing" />）なので、$\nabla(1+x)=1$, $\Delta(1+x)=1+2u$ で、$1+x$ からの距離はそれぞれ $x\le u/2$ と $2u-x\ge3u/2$ です。近いのは $1$ なので $\mathrm{fl}(1+x)=1$、続く引き算は $\mathrm{fl}(1-1)=0$ です。球算術では、$B(1,0)\oplus B(x,0)$ の中心は $\mathrm{fl}(1+x)=1$、半径は $r\ge0+0+u\cdot1=u$。そこから $B(1,0)$ を引くと中心 $\mathrm{fl}(1-1)=0$、半径 $r'\ge r+0+u\cdot0\ge u$ で、$|x|\le u/2\lt r'$ から $x\in B(0,r')$ です。
</Proof>

4 つの項目の読み方をまとめます。

- 1 は、球算術の結果が**証明として使える**ことを言っています。式全体についての包含は、区間演算の包含原理（<Ref to="mathematics/computer-assisted-proofs/floating-point-interval-arithmetic#thm-inclusion" />）と同じく、演算ごとの包含を積み重ねて得られます。
- 2 は代償です。球算術は 2 つの引数が同じ変数から来たことを知らないので、半径は過大になりえます。包含は保たれますが、精度についての情報が薄まります。式を書き換えて同じ変数の出現を減らすのが対策です（演習）。
- 3 は、丸めから来る半径が $u$ に比例することを言っています。1 回の演算についての主張ですが、演算を重ねても「半径が $u$ の定数倍」という形は、定義の式を通じて次の演算に受け継がれます。作業精度を倍にすれば $u$ は 2 乗になり、丸めに由来する半径は一気に縮みます。入力がもともと幅をもつ場合（データの誤差）は、その分の半径は精度を上げても縮みません。
- 4 は、「2 つの精度で結果が一致したから正しい」という推論の反例です。両方の精度で同じ向きに外れることがあります。球算術はこのとき半径 $u$ 以上の球を返し、「この精度では何もわからない」ことを正直に伝えます。

次の例の出力は、区間の下端と上端 $[\ell,r]$ を、中心 $(\ell+r)/2$ と半径 $(r-\ell)/2$ に直して表示しています。半径が小さいほど、真値が狭い範囲に絞られています。半径が中心と同じ大きさなら、真値が $0$ から中心の 2 倍の間のどこかにある、としかわからず、その精度では桁が 1 つも決まっていません。

<Example id="ex-interval-one-minus-cos" title="(1 − cos x)/x² を区間演算で評価する">
$(1-\cos x)/x^2$ は $x\to0$ で $1/2$ に近づきますが、素朴に計算すると桁落ちします（<Ref to="prop-digits-lost" /> の 1 で $a=1$, $b=\cos x$）。mpmath の区間演算 `iv` で、作業精度を 53, 106, 212 ビットに変えて評価し、得られた区間の中心と半径を、150 桁の参照値と比べます。この例を実行した環境（Python 3.9.6）には球算術ライブラリ python-flint を入れていないため、ここでは端点で持つ区間演算を使い、中心と半径に直して表示しています（`arb` による同じ計算は、例のあとに別の環境での実行結果を載せます）。包含の保証は同じです（<Ref to="mathematics/computer-assisted-proofs/floating-point-interval-arithmetic#thm-inclusion" />）。

```python
# Python 3.9.6 / mpmath 1.3.0（arm64）
import math
import mpmath as mp
from mpmath import iv

mp.mp.dps = 150                                               # 参照値は 150 桁
for xs in ("1e-3", "1e-6", "1e-9"):
    exact = (1 - mp.cos(mp.mpf(xs))) / mp.mpf(xs) ** 2
    xf = float(xs)
    print(f"x={xs} float64 の素朴な計算 {(1 - math.cos(xf)) / xf**2!r}")
    for prec in (53, 106, 212):
        iv.prec = prec
        x = iv.mpf(xs)                                        # 10 進の x を含む区間
        y = (1 - iv.cos(x)) / x**2
        mid, rad = y.mid, y.delta / 2                         # 中心と半径
        err = abs(mp.mpf(mid.a) - exact)
        inside = y.a <= exact <= y.b
        print(f"  prec={prec:>3} 中心={mp.nstr(mp.mpf(mid.a), 12)}"
              f" 半径={mp.nstr(mp.mpf(rad.b), 3)} 真の誤差={mp.nstr(err, 3)} 包含={inside}")
iv.prec = 53
x = iv.mpf(["0.1", "0.2"])
print(x - x)
print(x * (1 - x), x - x * x)
```

```text
x=1e-3 float64 の素朴な計算 0.49999995832550326
  prec= 53 中心=0.499999958381 半径=5.55e-11 真の誤差=4.77e-11 包含=True
  prec=106 中心=0.499999958333 半径=6.16e-27 真の誤差=2.76e-27 包含=True
  prec=212 中心=0.499999958333 半径=7.6e-59 真の誤差=6.68e-59 包含=True
x=1e-6 float64 の素朴な計算 0.5000444502911705
  prec= 53 中心=0.49998893914 半径=5.55e-5 真の誤差=1.11e-5 包含=True
  prec=106 中心=0.5 半径=6.16e-21 真の誤差=2.54e-21 包含=True
  prec=212 中心=0.5 半径=7.6e-53 真の誤差=6.99e-53 包含=True
x=1e-9 float64 の素朴な計算 0.0
  prec= 53 中心=55.5111512313 半径=55.5 真の誤差=55.0 包含=True
  prec=106 中心=0.5 半径=6.16e-15 真の誤差=1.96e-15 包含=True
  prec=212 中心=0.5 半径=7.6e-47 真の誤差=4.95e-47 包含=True
[-0.10000000000000001943, 0.10000000000000001943]
[0.07999999999999997391, 0.18000000000000002109] [0.059999999999999983902, 0.19000000000000002998]
```

**検算.** (a) 包含: 9 通りすべてで、参照値は区間の中にあり、中心の真の誤差は半径以下です。(b) 別精度: 作業精度を 53 から 106 ビットに上げると、半径は $5.55\times10^{-11}$ から $6.16\times10^{-27}$ へ、ほぼ $2^{-53}\approx1.1\times10^{-16}$ 倍に縮みます。<Ref to="prop-ball-vs-float" /> の 3 が 1 回の演算について述べた「半径は $u$ に比例する」という挙動です。(c) 桁落ちとの対応: $x=10^{-9}$ の倍精度の素朴な計算は $0.0$ を返し、これが誤りであることは結果からはわかりません。同じ 53 ビットの区間演算は中心 $55.5$・半径 $55.5$、すなわち区間 $[0,111]$ を返し、この精度では何も言えないことを示します。106 ビットでは半径 $6\times10^{-15}$ で、15 桁が保証されます。

最後の 2 行は依存性の例です。$X=[0.1,0.2]$ について $X-X$ は $\{0\}$ ではなく $[-0.1,0.1]$ になり、同じ関数 $x(1-x)=x-x^2$ の 2 つの書き方が違う幅の区間を返します（演習）。

なお、参照値をはじめ 60 桁で計算したところ、212 ビット（約 64 桁）の行で「包含 = False」と表示されました。区間演算が誤ったのではなく、**参照値の方が精度不足**だったのです。参照値は、比べる相手より十分に高い精度で用意します。
</Example>

python-flint が使える環境では、同じ計算を `arb` 型（中心と半径を持つ球）で書けます。球は区間より格納が軽く（高精度の中心と、低精度の半径）、高い精度での計算に向きます。$x=10^{-9}$ の場合を、python-flint を入れた別の環境（このテーマの後半の章と同じ Python 3.12 の環境。版はコードの先頭）で実行した結果を載せます。

```python
# Python 3.12.14 / python-flint 0.9.0（arm64）
from flint import arb, ctx

def g(x):
    return (1 - x.cos()) / (x * x)

ctx.prec = 600                                    # 参照の球（600 ビット）
ref = g(arb("1e-9"))
for prec in (53, 106, 212):
    ctx.prec = prec
    y = g(arb("1e-9"))                            # 10 進の 1e-9 を含む球から出発する
    err = abs(y.mid() - ref)                      # 中心の真の誤差（を含む球）
    print(f"prec={prec:>3} 半径={float(y.rad()):.3e} 真の誤差={float(err.mid()):.3e}"
          f" 包含={y.contains(ref)}  {y}")
ctx.prec = 53
x = arb("0.15", "0.05")                           # 中心 0.15、半径 0.05 の球
print(x - x, x * (1 - x), x - x * x)
```

```text
prec= 53 半径=5.000e-01 真の誤差=5.000e-01 包含=True  [+/- 0.501]
prec=106 半径=1.237e-14 真の誤差=8.125e-15 包含=True  [0.5000000000000 +/- 2.05e-14]
prec=212 半径=1.525e-46 真の誤差=2.643e-47 包含=True  [0.499999999999999999958333333333333333334722222 +/- 4.02e-46]
[+/- 0.101] [0.1 +/- 0.0801] [0.1 +/- 0.0951]
```

**検算.** (a) 包含: 3 つの精度のすべてで、600 ビットの参照の球は結果の球に含まれ、中心の真の誤差は半径以下です。(b) 別精度: 精度を 106 から 212 ビットに倍にすると、半径は $1.2\times10^{-14}$ から $1.5\times10^{-46}$ へ、$2^{-106}\approx1.2\times10^{-32}$ 倍に縮みます（<Ref to="prop-ball-vs-float" /> の 3 の挙動）。(c) 別実装との比較: 53 ビットでは球は中心 $0$・半径 $0.5$ で、mpmath の区間 $[0,111]$ より狭いものの、やはりこの精度では 1 桁も決まりません。半径の値は実装（初等関数の誤差の見積もり方）によるので、mpmath の `iv` の半径とは一致しません。どちらも真値を含むことだけが共通の保証です。最後の行は依存性の例で、球 $B(0.15,0.05)$ について $x-x$ は半径 $0.1$ の球になり（同 2）、$x(1-x)$ と $x-x^2$ は違う半径を返します。

数値の主張を証明の一部として使うときの手順は、[浮動小数と区間演算](/mathematics/computer-assisted-proofs/floating-point-interval-arithmetic)から始まる計算機援用証明のテーマで扱います。

## 7. 2 つの精度で比べる

包含を保証するライブラリを使わないとき、丸め誤差を見積もる手軽な方法は、同じ計算を 2 つの精度で行って差をとることです。これが何を教え、何を教えないかを命題にします。

<Proposition id="prop-two-precision-estimate" title="2 精度の差は丸め誤差の推定値">
ある式（打ち切りを含んでよい）の厳密な値を $y_h$、真に求めたい値を $y$ とする。同じ式を精度 $p$ と精度 $p'\gt p$ で計算した値を $\hat y_p=y_h+\rho_p$, $\hat y_{p'}=y_h+\rho_{p'}$ と書く（$\rho$ は丸め誤差）。

1. $0\le\lambda\lt1$ について $|\rho_{p'}|\le\lambda|\rho_p|$ ならば
   $$
   (1-\lambda)\,|\rho_p|\le|\hat y_p-\hat y_{p'}|\le(1+\lambda)\,|\rho_p|.
   $$
   すなわち、2 つの結果の差は、低い方の精度の丸め誤差に相対 $\lambda$ 以内で一致する。
2. 丸め誤差が単位丸めに比例する範囲、すなわち定数 $0\lt c_1\le c_2$ があって両方の精度で $c_1u\le|\rho|\le c_2u$（$u$ はそれぞれの単位丸め $2^{-p}$, $2^{-p'}$）のとき、$p'=2p$ ならば 1 は $\lambda=(c_2/c_1)\,2^{-p}$ で成り立つ。
3. 差 $\hat y_p-\hat y_{p'}=\rho_p-\rho_{p'}$ は $y$ にも $y_h-y$ にもよらない。とくに、2 つの結果が何桁一致していても、打ち切り誤差 $|y_h-y|$ については何もわからず、$|\hat y_p-y|\ge|y_h-y|-|\rho_p|$ である。
</Proposition>

<Proof of="prop-two-precision-estimate">
1. $\hat y_p-\hat y_{p'}=\rho_p-\rho_{p'}$ なので、三角不等式から $|\rho_p|-|\rho_{p'}|\le|\hat y_p-\hat y_{p'}|\le|\rho_p|+|\rho_{p'}|$ で、$|\rho_{p'}|\le\lambda|\rho_p|$ を代入します。
2. $|\rho_{2p}|\le c_2\,2^{-2p}$ と $|\rho_p|\ge c_1\,2^{-p}$ から $|\rho_{2p}|\le(c_2/c_1)\,2^{-p}\,|\rho_p|$ です。
3. 前半は 1 の最初の等式です。後半は $\hat y_p-y=(y_h-y)+\rho_p$ と三角不等式からです。
</Proof>

2 の仮定のうち下界 $c_1u\le|\rho|$ は、低い精度の丸め誤差がたまたま $0$ になる場合を除くためのものです。<Ref to="prop-ball-vs-float" /> の 4 は、仮定「$|\rho_{p'}|\le\lambda|\rho_p|$」が破れる例になっています。そこでは両方の精度で $\rho=-x$ と同じ誤差が生じ、$\lambda\lt1$ がとれません。丸め誤差が $u$ に比例するのは、作業精度が「条件数の桁」を超えてからです（<Ref to="thm-working-precision" />）。2 精度比較は、**少なくとも一方の精度が十分であるとき**にだけ、誤差の推定になります。

<Example id="ex-two-precision" title="2 精度比較が当たる例と外れる例">
前半は <Ref to="prop-ball-vs-float" /> の 4 の例（$x=10^{-40}$ を 53 ビットと 106 ビットで）、後半は交代級数の $e^{-20}$ を 60 ビットと 120 ビットで計算して差をとる例です。

```python
# Python 3.9.6 / mpmath 1.3.0（arm64）
import mpmath as mp
from mpmath import iv

x_str = "1e-40"
for prec in (53, 106):
    mp.mp.prec = prec
    x = mp.mpf(x_str)
    print(prec, (1 + x) - 1)                     # 浮動小数: どちらの精度でも 0
    iv.prec = prec
    y = (1 + iv.mpf(x_str)) - 1                  # 区間演算
    print("   区間", mp.nstr(mp.mpf(y.a), 3), mp.nstr(mp.mpf(y.b), 3))
mp.mp.prec = 160
print(160, mp.nstr((1 + mp.mpf(x_str)) - 1, 5))

# 2 つの精度の差による丸め誤差の推定（交代級数の e^{-20}）
def series(x):
    s = term = mp.mpf(1)
    k = 0
    while True:
        k += 1
        term = term * (-x) / k
        if s + term == s and k > x:
            return s
        s += term
vals = {}
for prec in (60, 120, 400):
    mp.mp.prec = prec
    vals[prec] = series(mp.mpf(20))
mp.mp.prec = 400
est = abs(vals[60] - vals[120])
true = abs(vals[60] - mp.exp(-20))
print("推定", mp.nstr(est, 6), "真の誤差", mp.nstr(true, 6), "比 - 1 =", mp.nstr(est / true - 1, 3))
```

```text
53 0.0
   区間 0.0 2.22e-16
106 0.0
   区間 0.0 2.47e-32
160 1.0e-40
推定 5.34785e-11 真の誤差 5.34785e-11 比 - 1 = -5.79e-19
```

**検算.** (a) 外れる例: $(1+10^{-40})-1$ は 53 ビットでも 106 ビットでも $0.0$ で、2 つの結果は完全に一致しますが、真値は $10^{-40}$ です（$10^{-40}\lt2^{-107}$ なので、命題の $x=2^{-q}$ の場合と同じ理由で $\mathrm{fl}(1+x)=1$ になります）。区間演算は $[0,2.2\times10^{-16}]$, $[0,2.5\times10^{-32}]$ を返し、真値を含みます。160 ビットにして初めて $10^{-40}$ が見えます。(b) 当たる例: 60 ビットで計算した $e^{-20}$ の誤差 $5.3\times10^{-11}$（真値 $2.06\times10^{-9}$ に対して 1〜2 桁しか合っていない）を、120 ビットの結果との差が 18 桁の精度で言い当てています。比のずれ $5.8\times10^{-19}$ は $2^{-60}\approx8.7\times10^{-19}$ と同じ大きさで、<Ref to="prop-two-precision-estimate" /> の 2 の $\lambda$ に対応します。(c) 打ち切り: <Ref to="ex-zeta3-budget" /> では、倍精度の `fsum` と 40 桁の部分和の差は $9\times10^{-17}$ でしたが、$\zeta(3)$ との差は $5\times10^{-13}$ でした。精度を変えても同じ部分和 $S_N$ を計算しているので、打ち切り誤差は 2 精度比較では見えません（命題の 3）。打ち切りは、$N$ や刻みを変えて調べます。
</Example>

刻みを変えたときの誤差の次数を使う検算と、数値実験の結果を証明に昇格させる手順は、このテーマの検算と検証の章で扱います。

## 8. Julia との対比

同じ道具立ては Julia にもあります。対応を表にします。Julia のコードはこの記事の実行環境では実行していないので、出力は載せません。

| 用途 | Python | Julia |
|---|---|---|
| 任意精度浮動小数 | mpmath の `mpf`（`mp.prec`, `mp.dps`） | `BigFloat`（`setprecision` で精度を指定） |
| 10 進の数を正確に取り込む | `mp.mpf("0.1")` | `big"0.1"`、`parse(BigFloat, "0.1")` |
| 厳密な有理数 | `fractions.Fraction(1, 10)` | `1//10`（`Rational`） |
| 高精度の円周率 | `mp.pi` | `BigFloat(π)` |
| 球算術 | python-flint の `arb` | Arblib.jl の `Arb` |

落とし穴も同じです。Julia でも、浮動小数リテラル `0.1` は binary64 の数なので、`BigFloat(0.1)` は 5.2 節の `mp.mpf(0.1)` と同じく「binary64 の $0.1$」を高精度に広げたものになります。

```julia
setprecision(BigFloat, 200) do
    a = parse(BigFloat, "0.1")   # 10 進の 0.1 を 200 ビットに丸めた数
    b = BigFloat(0.1)            # binary64 の 0.1 を 200 ビットに広げた数
    a - b                        # 0 にはならない（データの誤差）
end
```

ここで `big"0.1"` ではなく `parse` を使っているのには理由があります。`big"..."` はマクロなので、文字列を数に直すのはコードが読み込まれた時点で、そのときの大域的な精度（既定は 256 ビット）が使われます。`setprecision(...) do ... end` のブロックの中に書いても、ブロックの精度にはなりません（Julia の公式ドキュメントの `@big_str` の注意書き）。実行時の精度で作りたい数は `parse(BigFloat, "0.1")` で作ります。

言語が変わっても、誤差予算の立て方（2 節）と作業桁数の決め方（3 節）は変わりません。変わるのは、精度の指定が Python では大域的な設定（`mp.dps`）、Julia ではブロック単位（`setprecision(...) do ... end`）や型で行える、といった書き方の部分です。

## 9. 演習

<Exercise id="exr-working-digits" difficulty="易">
相対条件数（の近傍での上限）が $K=10^{8}$ の問題を、後退誤差が $10u$ 以下の算法で解く。(1) binary64 と binary32 で得られる桁数を <Ref to="thm-working-precision" /> から見積もれ。(2) 12 桁の精度が必要なとき、mpmath の `mp.dps` をいくつ以上にすればよいか。
<Solution>
(1) 定理の 2 から相対誤差は $cKu=10\cdot10^8\cdot u$ 以下です。得られる桁数は少なくとも $D-\log_{10}K-\log_{10}c=D-9$ で、binary64（$D=15.95$）では約 7 桁、binary32（$D=7.22$）では $-1.8$ 桁、すなわち 1 桁も保証されません。

(2) 定理の 3 から $D\ge12+8+1=21$ で足ります。`mp.dps = 21` 以上とします。条件数の見積もりの粗さを考えて数桁の余裕を足し、`mp.dps = 25` 程度にして、さらに `mp.dps = 50` で再計算して 12 桁が一致することを確かめるのが 3.2 節の手順です。
</Solution>
</Exercise>

<Exercise id="exr-digits-lost-cos" difficulty="標準">
$x=10^{-4}$ で $1-\cos x$ を倍精度で計算するとき、<Ref to="prop-digits-lost" /> の 1 から失う桁数を見積もり、得られる有効桁を予測せよ。また、同じ量を倍精度で 15 桁以上得る方法を 1 つ挙げよ。
<Solution>
$a=1$, $b=\cos x$ として、$|a|+|b|=1+\cos x\approx2$、$|a-b|=1-\cos x\approx x^2/2=5\times10^{-9}$ です。倍率は $2/(5\times10^{-9})=4\times10^{8}$ で、失う桁数は $\log_{10}(4\times10^8)\approx8.6$ 桁です。$\cos x$ の計算値の相対誤差が $u$ 程度とすると、有効桁は約 $15.95-8.6\approx7.3$ 桁、相対誤差にして $4\times10^8\times1.1\times10^{-16}\approx4\times10^{-8}$ 以下と予測できます。[浮動小数点数と丸め誤差](/mathematics/numerical-analysis/floating-point-arithmetic)の章の数値例では、$x=10^{-4}$ での実際の相対誤差は $5.2\times10^{-9}$ で、この範囲に入っています。

15 桁以上得るには、$1-\cos x=2\sin^2(x/2)$ と書き換えます。引き算がなくなり、各演算の相対誤差は $u$ 程度のままです。精度を上げる（$16+9=25$ 桁以上で計算する）方法もありますが、書き換えの方が安くつきます。
</Solution>
</Exercise>

<Exercise id="exr-doubling-cost" difficulty="標準">
<Ref to="prop-precision-doubling" /> の 2 で、(1) $M(q)=q$（費用が桁数に比例）、(2) $M(q)=q^{\log_23}$（Karatsuba 法による乗算。$\log_23\approx1.585$）のそれぞれについて、$J\to\infty$ での費用の合計と最後の 1 回の費用の比 $\sum_{i\ge0}M(p/2^i)/M(p)$ を求めよ。どちらも $M(2q)\ge2M(q)$ を満たすことも確かめよ。
<Solution>
(1) $M(2q)=2q=2M(q)$ で、仮定を等号で満たします。$M(p/2^i)/M(p)=2^{-i}$ なので、比は $\sum_{i\ge0}2^{-i}=2$ です。命題の上界 2 は、この場合にちょうど達成されます。

(2) $M(2q)=(2q)^{\log_23}=2^{\log_23}q^{\log_23}=3M(q)\ge2M(q)$ です。$M(p/2^i)/M(p)=(2^{-i})^{\log_23}=3^{-i}$ なので、比は $\sum_{i\ge0}3^{-i}=3/2$ です。

乗算の費用が桁数に対して急に増えるほど、低い精度での反復は相対的に安くなり、比は 1 に近づきます。本文の例の筆算モデル $M(q)=q^2$ では $\sum4^{-i}=4/3$ でした。
</Solution>
</Exercise>

<Exercise id="exr-dependency" difficulty="標準">
区間 $X=[0.1,0.2]$ について、厳密な区間演算（丸めなし）で (1) $X\cdot(1-X)$、(2) $X-X\cdot X$ を計算し、関数 $g(x)=x(1-x)$ の $X$ 上の真の値域と比べよ。どちらの書き方が狭い区間を与えるか。真の値域をちょうど与える書き方はあるか。
<Solution>
(1) $1-X=[0.8,0.9]$。正の区間どうしの積は端点の積で、$X\cdot(1-X)=[0.1\cdot0.8,\ 0.2\cdot0.9]=[0.08,0.18]$ です。

(2) $X\cdot X=[0.01,0.04]$（区間演算は 2 つの $X$ を独立に動かすので、2 乗としては過大ではありませんが、次の引き算で依存性が失われます）。$X-X\cdot X=[0.1-0.04,\ 0.2-0.01]=[0.06,0.19]$ です。

真の値域: $g'(x)=1-2x\gt0$（$x\lt1/2$）なので $g$ は $X$ 上で増加し、値域は $[g(0.1),g(0.2)]=[0.09,0.16]$ です。どちらの結果もこれを含みますが、幅は (1) が $0.10$、(2) が $0.13$、真の値域が $0.07$ で、どちらも過大です。本文の例の最後の行（外向きに丸めた `iv` の結果）は、この (1), (2) と末尾の桁を除いて一致しています。

$x$ が 1 回しか現れない形に書き換えると、依存性の問題は起きません。$g(x)=\tfrac14-(x-\tfrac12)^2$ とすると、$X-\tfrac12=[-0.4,-0.3]$、その 2 乗は $[0.09,0.16]$（負の区間の 2 乗は端点の 2 乗を入れ替えたもの）、$\tfrac14-[0.09,0.16]=[0.09,0.16]$ で、真の値域に一致します。
</Solution>
</Exercise>

## 参考文献

- N. J. Higham, *Accuracy and Stability of Numerical Algorithms*, 2nd ed., SIAM, 2002 — 第 1 章（誤差の種類、前進誤差・後退誤差・条件数、桁落ち）、第 12 章（反復改良）。
- R. P. Brent, P. Zimmermann, *Modern Computer Arithmetic*, Cambridge University Press, 2010 — 任意精度の乗算の費用と、Newton 法で精度を倍にしていく手法。
- W. Tucker, *Validated Numerics: A Short Introduction to Rigorous Computations*, Princeton University Press, 2011 — 区間演算と依存性の問題。
- F. Johansson, "Arb: efficient arbitrary-precision midpoint-radius interval arithmetic", *IEEE Transactions on Computers* 66 (2017), 1281–1292. [doi:10.1109/TC.2017.2690633](https://doi.org/10.1109/TC.2017.2690633) — 球算術の実装（python-flint の `arb` の基礎）。
- mpmath 公式ドキュメント, [Basic usage / Setting the precision](https://mpmath.org/doc/current/basics.html) — `mp.prec` と `mp.dps`、文字列からの数の生成、区間演算 `iv`。
- L. N. Trefethen, D. Bau III, *Numerical Linear Algebra*, SIAM, 1997 — Lecture 12（条件数）、Lecture 15（後退安定な算法の精度）。