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

0. この記事の要点

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

• 数値の主張に桁を付けるには、誤差を原因ごとに分けて上から抑えた誤差予算を作ります。計算値と真値の差は、丸め・打ち切り・入力データの 3 項の和で抑えられます。
• 相対条件数 $\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 つの精度で一致した桁」は丸め誤差の推定にはなりますが、包含の保証にはならず、打ち切り誤差は見えません。
• 高精度計算を台無しにする最も多い原因は、定数や入力を倍精度で作ってから渡すことです。高精度の数は文字列か有理数から作ります。

この記事で使う事前知識

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

1. 動機

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

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

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

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

条件数・前進誤差・後退誤差は条件数と安定性の章の定義を、浮動小数点数の誤差モデルは浮動小数点数と丸め誤差の章のものを使います。区間演算と球算術の定義と包含の定理は浮動小数と区間演算の章にあり、ここではその使い方の側を扱います。

2. 誤差予算

2.1. 用語

定義 2.1（作業精度・目標精度・有効桁）↗

• 計算に使う浮動小数点数の精度 $p$（ビット）を作業精度といい、$D=-\log_{10}u=p\log_{10}2$（$u=2^{-p}$ は単位丸め。定義 2.3［浮動小数点数と丸め誤差］）を作業桁数という。binary64 では $D\approx15.95$、binary32 では $D\approx7.22$ である。
• 最終結果に求める相対誤差の上限 $10^{-d}$ を目標精度、$d$ を目標の桁数という。
• 真値 $y\ne0$ に対する計算値 $\hat y$ の有効桁を $-\log_{10}(|\hat y-y|/|y|)$ と定める。

定義 2.2（誤差予算）↗

最終結果の誤差 $|\hat y-y|$ を、原因ごとの項の和で上から抑えた表を誤差予算という。項は次の 4 種類に分ける。

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

はじめの 3 つは決定的な上界で抑えられます。統計誤差は確率的な量で、標準誤差（再現できる数値実験の 命題 6.1［再現できる数値実験］）で見積もり、表では別の行にします。

2.2. 3 項への分解

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

命題 2.3（誤差の 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$。

証明（命題 2.3）

$\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$ はその上で微分可能（したがって連続）なので、平均値の定理（定理 3.3［平均値の定理とテイラーの定理］）から、両者の間のある $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$ です。

∎

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

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

例 2.4（ζ(3) の部分和の誤差予算）

$\zeta(3)=\sum_{k\ge1}k^{-3}$ を、$N=10^6$ 項の部分和 $S_N$ を倍精度で足して近似します。予算を先に立てます。

• 打ち切り: $f(x)=x^{-3}$ は $[1,\infty)$ で正かつ単調減少なので、積分判定法（定理 5.1［級数と収束判定］ の 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$ 個です。誤差因子の積の補題（補題 3.2［浮動小数と区間演算］。除数の側の丸めは $(1+\delta)^{-1}$ の形で入ります）から誤差は $\gamma_NS_N\approx1.3\times10^{-10}$ 以下です（逐次和の誤差 定理 5.1［浮動小数点数と丸め誤差］ と同じ議論）。小さい項から足す順序では、最初に足す項が加算 $N-1$ 回と丸め 2 回を受けるので因子は $N+1$ 個以下、上界は $\gamma_{N+1}S_N$ で、数値は同じ $1.3\times10^{-10}$ です。
• データ: 入力は整数 $k$ だけなので $0$ です。

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

打ち切り  真値 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}$ でした。

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

3. 必要な作業桁数

3.1. 条件数の桁だけ失う

丸めの項について、作業精度をどれだけにすればよいかを定理にします。既存の「前進誤差 ≤ 条件数 × 後退誤差」（定理 4.3［条件数と安定性］）を、桁数の形に直したものです。

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

定理 3.1（必要な作業桁数）✓

$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$ 以下。定義 4.2［条件数と安定性］）。

計算値 $\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})$（定義 2.1［条件数と安定性］）について $$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$ を足した作業桁数で足りる。

証明（定理 3.1）

1. 定理 4.3［条件数と安定性］ の 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$ に含まれます（球は凸）。定理 4.3［条件数と安定性］ の 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$ 以下です。

∎

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)$ 以下です（系 7.3［行列のノルムと線形方程式の条件数］。$\kappa(A)=\|A\|\|A^{-1}\|$ は 定義 6.2［行列のノルムと線形方程式の条件数］）。$\kappa(A)\eta\le1/2$ なら $4\kappa(A)\eta$ 以下なので、桁数の勘定は同じで、「$\log_{10}\kappa(A)$ 桁を失う」と読めます。

例 3.2（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 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))

 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}$ なので 定理 3.1 から 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 行列の条件数ではありません。

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

定理 3.1 を手順に直すと次のようになります。

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 で失う桁数の原因です。

補題 4.1（指数級数）✓

すべての実数 $x$ について $\displaystyle\sum_{k=0}^{\infty}\frac{x^k}{k!}=e^x$ である。

証明（補題 4.1）

$x=0$ なら両辺は $1$ です。$x\ne0$ とします。$\exp$ は $\mathbb{R}$ 上で何回でも微分可能で、導関数はすべて $\exp$ です（定理 5.5［導関数の定義と基本的な微分法］ の 4）。テイラーの定理（定理 5.3［平均値の定理とテイラーの定理］）を $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$ は狭義単調増加で（系 3.4［平均値の定理とテイラーの定理］ の 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$ に収束します。

∎

命題 4.2（差と和で失う桁数）✓

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$$ に近づく。

証明（命題 4.2）

1. 桁落ちの相対誤差の評価（命題 4.2［浮動小数点数と丸め誤差］）そのものです。相対誤差が $\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!$ なので、補題 4.1 から $\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$ です。

∎

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 桁も残りません。

例 4.3（交代級数で e⁻ˣ を計算する）

$x=20,30,40$ について、交代級数を倍精度と mpmath（dps が 30 と 60）で足し、有効桁を調べます。参照値は 100 桁の mp.exp(-x) です。

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

 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$ 桁）です。作業桁数から 命題 4.2 の「失う桁」を引いた予測は、倍精度で $-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 桁目から違います。

精度を上げるのは最後の手段です。この例では、式を $1/e^{x}$ に書き換えるだけで倍精度で足ります。書き換えが見つからないとき、あるいは書き換えた式の正しさを確かめたいときに、定理 3.1 と 命題 4.2 で必要な桁数を見積もって任意精度を使います。

5. 精度の道具

5.1. 何が使えるか

定義 5.1（任意精度浮動小数と厳密な有理数）↗

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

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

定義 5.2（混合精度）↗

1 つの計算の中で、部分ごとに異なる精度を使うことを混合精度という。費用の大半を占める部分を低い精度で、結果の精度を決める部分を高い精度で行う。

混合精度の代表例は、連立 1 次方程式の反復改良です。LU 分解と補正方程式の求解を低い精度で行い、残差 $\boldsymbol{b}-A\boldsymbol{x}_k$ だけを高い精度で計算すると、誤差は 1 回ごとに縮みます（定理 5.1［Gauss 消去と LU 分解］）。$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 は、同じコードが環境によって別の精度で動くので、再現性（再現できる数値実験）の観点から、結果を左右する場所には使わないのがよいと思います。「倍精度より少し多く」が必要なら、桁数を明示できる任意精度を使います。

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 を使います。

命題 2.3 の言葉では、これはデータの誤差です。作業精度を 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 倍未満で済むことを述べています。

命題 5.3（精度を倍にしながらの反復）✓

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 倍未満で済む。

証明（命題 5.3）

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

∎

Newton 法は 1 の仮定を満たします。$f$ が $C^2$ 級で $f(x^*)=0$, $f'(x^*)\ne0$ なら、$x^*$ の近くで $|x_{k+1}-x^*|\le K|x_k-x^*|^2$ です（定理 5.1［非線形方程式の反復解法］）。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$ に比例）でも高速な乗算法でも成り立ちます。

例 5.4（√2 を 3200 ビットまで）

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

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

初期値      正しいビット数    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$ 倍で、命題 5.3 の上界 2 倍の内側です。すべてを目標精度で行うと 6 倍です。ここでの費用は演算量のモデルによる比で、時間の計測値ではありません。

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

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

浮動小数の計算は、誤差がどれだけあるかを結果に含めません。区間演算と球算術は、真値を含む集合（区間 $[\ell,r]$、あるいは中心と半径の組 $B(m,r)=[m-r,m+r]$。定義 4.8［浮動小数と区間演算］）を結果として返します。

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

命題 6.1（球算術の保証と限界）✓

球 $B(m,r)$ と演算 $\oplus,\otimes$ は 定義 4.8［浮動小数と区間演算］ のものとし、$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$ を含む。

証明（命題 6.1）

1. 球算術の包含の定理（定理 4.9［浮動小数と区間演算］）そのものです。

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)$ で、四則の標準モデル（系 3.1［浮動小数点数と丸め誤差］）から $|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$ の等間隔な点列（補題 2.2［浮動小数と区間演算］）なので、$\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')$ です。

∎

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

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

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

例 6.2（(1 − cos x)/x² を区間演算で評価する）

$(1-\cos x)/x^2$ は $x\to0$ で $1/2$ に近づきますが、素朴に計算すると桁落ちします（命題 4.2 の 1 で $a=1$, $b=\cos x$）。mpmath の区間演算 iv で、作業精度を 53, 106, 212 ビットに変えて評価し、得られた区間の中心と半径を、150 桁の参照値と比べます。この例を実行した環境（Python 3.9.6）には球算術ライブラリ python-flint を入れていないため、ここでは端点で持つ区間演算を使い、中心と半径に直して表示しています（arb による同じ計算は、例のあとに別の環境での実行結果を載せます）。包含の保証は同じです（定理 4.4［浮動小数と区間演算］）。

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

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}$ 倍に縮みます。命題 6.1 の 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」と表示されました。区間演算が誤ったのではなく、参照値の方が精度不足だったのです。参照値は、比べる相手より十分に高い精度で用意します。

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

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

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}$ 倍に縮みます（命題 6.1 の 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$ は違う半径を返します。

数値の主張を証明の一部として使うときの手順は、浮動小数と区間演算から始まる計算機援用証明のテーマで扱います。

7. 2 つの精度で比べる

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

命題 7.1（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|$ である。

証明（命題 7.1）

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$ と三角不等式からです。

∎

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

例 7.2（2 精度比較が当たる例と外れる例）

前半は 命題 6.1 の 4 の例（$x=10^{-40}$ を 53 ビットと 106 ビットで）、後半は交代級数の $e^{-20}$ を 60 ビットと 120 ビットで計算して差をとる例です。

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

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}$ と同じ大きさで、命題 7.1 の 2 の $\lambda$ に対応します。(c) 打ち切り: 例 2.4 では、倍精度の fsum と 40 桁の部分和の差は $9\times10^{-17}$ でしたが、$\zeta(3)$ との差は $5\times10^{-13}$ でした。精度を変えても同じ部分和 $S_N$ を計算しているので、打ち切り誤差は 2 精度比較では見えません（命題の 3）。打ち切りは、$N$ や刻みを変えて調べます。

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

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$」を高精度に広げたものになります。

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. 演習

演習 9.1易

相対条件数（の近傍での上限）が $K=10^{8}$ の問題を、後退誤差が $10u$ 以下の算法で解く。(1) binary64 と binary32 で得られる桁数を 定理 3.1 から見積もれ。(2) 12 桁の精度が必要なとき、mpmath の mp.dps をいくつ以上にすればよいか。

解答

(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 節の手順です。

演習 9.2標準

$x=10^{-4}$ で $1-\cos x$ を倍精度で計算するとき、命題 4.2 の 1 から失う桁数を見積もり、得られる有効桁を予測せよ。また、同じ量を倍精度で 15 桁以上得る方法を 1 つ挙げよ。

解答

$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}$ 以下と予測できます。浮動小数点数と丸め誤差の章の数値例では、$x=10^{-4}$ での実際の相対誤差は $5.2\times10^{-9}$ で、この範囲に入っています。

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

演習 9.3標準

命題 5.3 の 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)$ を満たすことも確かめよ。

解答

(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$ でした。

演習 9.4標準

区間 $X=[0.1,0.2]$ について、厳密な区間演算（丸めなし）で (1) $X\cdot(1-X)$、(2) $X-X\cdot X$ を計算し、関数 $g(x)=x(1-x)$ の $X$ 上の真の値域と比べよ。どちらの書き方が狭い区間を与えるか。真の値域をちょうど与える書き方はあるか。

解答

(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]$ で、真の値域に一致します。

参考文献

• 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 — 球算術の実装（python-flint の arb の基礎）。
• mpmath 公式ドキュメント, Basic usage / Setting the precision — mp.prec と mp.dps、文字列からの数の生成、区間演算 iv。
• L. N. Trefethen, D. Bau III, Numerical Linear Algebra, SIAM, 1997 — Lecture 12（条件数）、Lecture 15（後退安定な算法の精度）。