浮動小数と区間演算：IEEE 754 の丸めモデルから区間 Newton 法まで

0. この記事の要点

計算機で 0.1 + 0.2 を計算すると、ぴったり 0.3 にはなりません。計算機は小数を有限の桁で丸めて持つからです。この記事では、その丸めの大きさを式で押さえ、「真の値は必ずこの範囲の中にある」と言い切れる計算のしかた（区間演算）を学びます。読み終えると、計算機の出力を、方程式の解が存在するという証明の一部として使う理由が分かります。

• 倍精度の浮動小数の四則と平方根は、IEEE 754 により正しく丸められます。そこから、正規化数の範囲では $\mathrm{fl}(x\circ y)=(x\circ y)(1+\delta)$, $|\delta|\le u=2^{-53}$ という誤差モデルが従います。
• 丸めが $n$ 回重なった誤差は Higham の $\gamma_n=nu/(1-nu)$ で抑えられ、内積の計算値は $|\widehat{s}-\boldsymbol{x}^{\mathsf{T}}\boldsymbol{y}|\le\gamma_n\sum_i|x_iy_i|$ を満たします。
• 区間演算は「真の値を含む区間」を演算の結果として返します。包含原理により、式を区間で評価した結果は値域を必ず含みます（ただし過大評価はあり得ます）。中心と半径で表す球算術（Arb の方式）も同じ原理で動きます。
• 丸めの向きを指定できない言語（Python・JavaScript など）では、最近点丸めの結果を nextafter で 1 つずつ外へ動かせば保守的な区間が得られます。これが正しいのは、IEEE 754 が正しい丸めを要求する演算（四則・平方根・fma）だけです。exp や sin などの数学関数には同じ保証がありません。
• 区間 Newton 法は平均値の定理と中間値の定理から、$N(X)\subset X$ なら $X$ に零点がちょうど 1 つあることを結論します。計算機の出力が、そのまま存在と一意性の証明になります。

この記事で使う事前知識

• 平均値の定理: 区間の両端の値の差は、途中のどこかの点での微分係数と区間の幅の積に等しい、という定理。→ 平均値の定理
• 中間値の定理: 連続関数は、2 点の値の間の値をすべてとる、という定理。→ 中間値の定理
• 凸集合: 集合の 2 点を結ぶ線分がいつも集合の中に収まる、という性質。→ 凸集合
• 連続写像: 近い点を近い点に写す写像。→ 連続写像

1. 動機

計算機援用証明は、証明の一部を計算機による数値計算に任せる方法です。Hales による Kepler 予想の証明、Tucker による Lorenz アトラクタの存在証明（Smale の第 14 問題）は、いずれも大量の区間演算を証明の一部として含みます。そのとき問題になるのは、浮動小数の計算結果は一般に真の値と一致しない、という点です。0.1 + 0.2 は倍精度で 0.30000000000000004 になります。計算結果を定理の一部として使うには、真の値がどこにあるかを数学的に保証しなければなりません。

このテーマ「計算機援用証明」の 4 章は、その保証の仕組みを下から積み上げます。

章	内容	保証されるもの
1（本章）	IEEE 754 の丸めモデル、誤差解析、区間演算、区間 Newton 法	個々の演算と方程式の解
2	Chebyshev 展開、Gauss–Legendre 求積の誤差、区間 Taylor 法	積分と関数近似
3	LDLᵀ・Cholesky と慣性による正定値性の認証（Rump の方法）、Sturm 列、Weyl・Bauer–Fike、LP・SDP の証明書	行列の固有値と最適値
4	Rayleigh–Ritz、Temple・Lehmann–Goerisch、有限要素法の誤差込みの下界	微分作用素の固有値

本章は全体の土台です。後の章の「計算機で確かめた」はすべて、本章の包含原理に帰着します。

2. 浮動小数と丸め

定義 2.1（2 進浮動小数の集合）↗

整数 $p\ge2$（精度）と $e_{\min}\lt e_{\max}$ を固定する。正規化数を $\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}$ を単位丸めという。

IEEE 754 の倍精度（binary64）は $p=53$, $e_{\min}=-1022$, $e_{\max}=1023$ で、$u=2^{-53}\approx1.11\times10^{-16}$ である。

正規化数 $M\cdot2^{e-p+1}$ は区間 $[2^e,2^{e+1})$ に入ります。$M$ が $1$ ずつ動くので、次がわかります。

補題 2.2（浮動小数の間隔）✓

$e_{\min}\le e\le e_{\max}$ のとき、$\mathbb{F}\cap[2^e,2^{e+1}]$ は間隔 $2^{e-p+1}=2u\cdot2^e$ の等間隔な点列で、両端 $2^e$, $2^{e+1}$ を含む（$e=e_{\max}$ のときは右端を除く）。$[0,2^{e_{\min}}]$ の中の $\mathbb{F}$ の元は間隔 $2^{e_{\min}-p+1}$ の等間隔な点列である。

証明（補題 2.2）

$[2^e,2^{e+1})$ に入る正規化数は、$2^{p-1}\le M\lt2^p$ に対する $M\cdot2^{e-p+1}$ で、$M=2^{p-1}$ のとき $2^e$ です。指数 $e'\ne e$ の正規化数は $[2^{e'},2^{e'+1})$ に入るのでこの区間には来ません。右端 $2^{e+1}$ は指数 $e+1$ の $M=2^{p-1}$ です。非正規化数と $0$ は $M\cdot2^{e_{\min}-p+1}$（$0\le M\lt2^{p-1}$）なので、指数 $e_{\min}$ の正規化数とつながって同じ間隔で並びます。

∎

定義 2.3（丸め）↗

$|x|\le x_{\max}$ の実数 $x$ に対し、

• $\nabla(x)=\max\{f\in\mathbb{F}:f\le x\}$（下向き丸め）、$\Delta(x)=\min\{f\in\mathbb{F}:f\ge x\}$（上向き丸め）、
• $\mathrm{fl}(x)$ は $\nabla(x)$, $\Delta(x)$ のうち $x$ に近い方（等距離なら $M$ が偶数の方。最近点丸め）

と定める。$\nabla(x)$ と $\Delta(x)$ のどちらかを返す丸めを忠実な丸めという。

IEEE 754 は、$+,-,\times,\div$ と $\sqrt{\ }$（および fma $a\times b+c$。積と和を丸め 1 回で行う演算）について、計算結果が厳密な実数の結果を指定の向きに丸めたものであることを要求します。既定の向きは最近点丸めで、$\mathrm{fl}(x\circ y)$ が返ります。この要求が、以下のすべての誤差評価の出発点です。

定理 2.4（丸めの誤差モデル）✓

$x_{\min}\le|x|\le x_{\max}$ とする。

1. $\mathrm{fl}(x)=x(1+\delta)$ となる $|\delta|\le u$ がある。また $|\mathrm{fl}(x)-x|\le u\,|\mathrm{fl}(x)|$ である。
2. $\nabla(x)$, $\Delta(x)$ はいずれも $x(1+\delta)$, $|\delta|\lt2u$ の形である。
3. $|x|\lt x_{\min}$ なら $|\mathrm{fl}(x)-x|\le2^{e_{\min}-p}$ である（倍精度で $2^{-1075}$）。

証明（定理 2.4）

対称性から $x\gt0$ とし、$2^e\le x\lt2^{e+1}$（$e_{\min}\le e\le e_{\max}$）とします。補題 2.2 より、$x$ をはさむ隣り合う浮動小数 $\nabla(x)\le x\le\Delta(x)$ の差は $2u\cdot2^e$ 以下です。

1. 最近点丸めは近い方を選ぶので $|\mathrm{fl}(x)-x|\le u\cdot2^e\le u\,x$ で、$\delta=(\mathrm{fl}(x)-x)/x$ とおけば $|\delta|\le u$ です。また $2^e\in\mathbb{F}$ かつ $2^e\le x$ なので $\nabla(x)\ge2^e$、したがって $\mathrm{fl}(x)\ge2^e$ で、$|\mathrm{fl}(x)-x|\le u\cdot2^e\le u\,\mathrm{fl}(x)$ です。
2. $x\notin\mathbb{F}$ なら $|\nabla(x)-x|$, $|\Delta(x)-x|$ はともに間隔 $2u\cdot2^e$ より真に小さく、$2u\cdot2^e\le2u\,x$ です。$x\in\mathbb{F}$ なら $\delta=0$ です。
3. $[0,x_{\min}]$ の浮動小数の間隔は $2^{e_{\min}-p+1}$ なので、最近点との距離はその半分以下です。

∎

3 はアンダーフローの場合で、相対誤差ではなく絶対誤差 $2^{-1075}$ しか保証されません。以下では簡単のため、途中の値が $0$ か正規化数の範囲に入る（アンダーフローもオーバーフローも起こらない）と仮定します。実際の検証プログラムでは絶対誤差の項を足すか、範囲外の値が出たら検証失敗として止めます。

例 2.5（0.1 は 1/10 ではない）

十進の 0.1 は、読み込みの時点で $\mathrm{fl}(1/10)$ に丸められます。Python の fractions.Fraction(0.1) で厳密な値を見ると $\mathrm{fl}(1/10)-1/10=1/180143985094819840\approx5.55\times10^{-18}$ で、定理 2.4 の 1 の通り $u\cdot\frac1{10}\approx1.1\times10^{-17}$ 以下です。したがって「$1/10$ を含む区間」が要るときは、0.1 そのものではなく、その前後の浮動小数で区間を作るか、有理数で扱います。

3. 丸め誤差解析と Higham の γₙ

丸めの $\delta$ が $n$ 個掛け合わさったとき、全体の相対誤差をまとめて抑える記号を用意します。足し算や掛け算を何十回も重ねた計算で、誤差が結局どれくらいになるのかを一つの式で言うための準備です。次の補題は、$n$ 回分の誤差因子の積が $1\pm\gamma_n$ に収まることを、$n$ についての帰納法で示します。

定義 3.1（γₙ）↗

$nu\lt1$ のとき $\gamma_n=\dfrac{nu}{1-nu}$ とおく。

補題 3.2（誤差因子の積）✓

$|\delta_i|\le u$, $\rho_i\in\{1,-1\}$（$i=1,\ldots,n$）、$nu\lt1$ とする。このとき

$$\prod_{i=1}^n(1+\delta_i)^{\rho_i}=1+\theta_n,\qquad|\theta_n|\le\gamma_n$$

である。

証明（補題 3.2）

$P_k=\prod_{i\le k}(1+\delta_i)^{\rho_i}$ とし、$|P_k-1|\le\gamma_k$ を $k$ の帰納法で示します。$k=1$ は、$\rho_1=1$ なら $|\delta_1|\le u\le\gamma_1$、$\rho_1=-1$ なら $|1/(1+\delta_1)-1|=|\delta_1|/|1+\delta_1|\le u/(1-u)=\gamma_1$ です。$k$ で成り立つとし、$\delta=\delta_{k+1}$ とします。

$\rho_{k+1}=1$ のとき、$P_k(1+\delta)-1=(P_k-1)(1+\delta)+\delta$ なので

$$|P_{k+1}-1|\le\gamma_k(1+u)+u=\frac{ku(1+u)+u(1-ku)}{1-ku}=\frac{(k+1)u}{1-ku}\le\gamma_{k+1}.$$

$\rho_{k+1}=-1$ のとき、$P_k/(1+\delta)-1=(P_k-1-\delta)/(1+\delta)$ なので

$$|P_{k+1}-1|\le\frac{\gamma_k+u}{1-u}=\frac{(k+1-ku)\,u}{(1-ku)(1-u)}.$$

これが $\gamma_{k+1}=(k+1)u/(1-(k+1)u)$ 以下であることは、分母を払った

$$(k+1)(1-ku)(1-u)-(k+1-ku)\bigl(1-(k+1)u\bigr)=ku\ge0$$

（両辺を展開すると $u^0$ と $u^2$ の項は打ち消し合い、$u^1$ の係数が $-(k+1)-k(k+1)+(k+1)^2+k=k$ になる）からわかります。

∎

定理 3.3（内積の丸め誤差）✓

$\boldsymbol{x},\boldsymbol{y}\in\mathbb{F}^n$、$nu\lt1$ とし、$\widehat{s}$ を $s_1=\mathrm{fl}(x_1y_1)$, $s_i=\mathrm{fl}\bigl(s_{i-1}+\mathrm{fl}(x_iy_i)\bigr)$ で計算した $s_n$ とする。途中でアンダーフロー・オーバーフローが起きなければ

$$\widehat{s}=\sum_{i=1}^nx_iy_i(1+\theta^{(i)}),\quad|\theta^{(i)}|\le\gamma_n,\qquad\text{したがって}\quad|\widehat{s}-\boldsymbol{x}^{\mathsf{T}}\boldsymbol{y}|\le\gamma_n\sum_{i=1}^n|x_iy_i|.$$

証明（定理 3.3）

定理 2.4 の 1 により、各乗算と加算の結果は $(1+\varepsilon_i)$, $(1+\delta_i)$（$|\varepsilon_i|,|\delta_i|\le u$）を掛けたものです（結果が $0$ なら丸めは厳密なので $\delta=0$ とします）。すなわち $s_1=x_1y_1(1+\varepsilon_1)$, $s_i=\bigl(s_{i-1}+x_iy_i(1+\varepsilon_i)\bigr)(1+\delta_i)$ です。これを展開すると

$$\widehat{s}=x_1y_1(1+\varepsilon_1)\prod_{j=2}^n(1+\delta_j)+\sum_{i=2}^nx_iy_i(1+\varepsilon_i)\prod_{j=i}^n(1+\delta_j)$$

です。第 1 項の因子は $n$ 個、第 $i$ 項の因子は $n-i+2\le n$ 個です。補題 3.2 により、$k$ 個の因子の積は $1+\theta_k$, $|\theta_k|\le\gamma_k\le\gamma_n$ と書けます（$\gamma_k$ は $k$ について増加）。最後の不等式は三角不等式です。

∎

評価が $\sum|x_iy_i|$ に比例することが要点です。打ち消し合いで $|\boldsymbol{x}^{\mathsf{T}}\boldsymbol{y}|$ が $\sum|x_iy_i|$ よりずっと小さいと、相対誤差は大きくなり得ます。また、足す順番を変えても各項が受ける丸めは $n$ 回以下なので、同じ評価が成り立ちます。この評価は第 3 章の Cholesky 分解の誤差解析で使います。

例 3.4（0.1 を 10 回足す）

s = 0.0 から s = s + 0.1 を 10 回繰り返すと 0.9999999999999999 になります。足しているのは $f=\mathrm{fl}(1/10)$ なので、厳密な和は $10f$ で、差は有理数で計算すると約 $1.67\times10^{-16}$ です。定理 3.3（$y_i=1$ とし、乗算は厳密）の評価は $\gamma_{10}\cdot10f\approx1.1\times10^{-15}$ で、確かにその範囲に入っています。なお Python 3.12 以降の組み込み sum は補償つきの加算を使い、NumPy の sum は対ごとの加算を使うので、同じ実験を sum で行うと結果が変わります。誤差解析の対象にするアルゴリズムを、コードの上で明示しておくことが大切です。

4. 区間演算と球算術

4.1. 区間演算の包含原理

実数の閉区間 $X=[\underline{x},\overline{x}]$（$\underline{x}\le\overline{x}$）の全体を $\mathbb{IR}$ と書きます。幅を $w(X)=\overline{x}-\underline{x}$ とします。

定義 4.1（区間の四則）↗

$X,Y\in\mathbb{IR}$ と $\circ\in\{+,-,\times,\div\}$ に対し $X\circ Y=\{x\circ y:x\in X,\ y\in Y\}$ と定める（$\div$ は $0\notin Y$ のときだけ定義する）。

命題 4.2（端点による公式）✓

$X=[a,b]$, $Y=[c,d]$ について $X+Y=[a+c,b+d]$, $X-Y=[a-d,b-c]$, $X\times Y=[\min P,\max P]$（$P=\{ac,ad,bc,bd\}$）、$0\notin Y$ なら $X\div Y=X\times[1/d,1/c]$ である。

証明（命題 4.2）

積について示します。$g(x,y)=xy$ は $y$ を固定すると $x$ の 1 次関数なので、$x\in[a,b]$ での最大・最小は端点でとられ、さらに $x$ を端点に固定すると $y$ の 1 次関数なので同様です。よって $X\times Y$ の上限・下限は $P$ の最大・最小で、それらはとられます。最小点 $\boldsymbol{p}$ と最大点 $\boldsymbol{q}$ を結ぶ線分は凸集合 $X\times Y$ に含まれ、$t\mapsto g((1-t)\boldsymbol{p}+t\boldsymbol{q})$ は連続なので、中間値の定理（定理 5.7［極限と連続性］）からその間の値はすべてとられます。和と差は各変数について単調なので同様です。商は $\{1/y:y\in Y\}=[1/d,1/c]$（$1/y$ は $0$ を含まない区間で単調減少）からです。

∎

定義 4.3（自然な区間拡張）↗

変数 $x_1,\ldots,x_m$ と定数から $+,-,\times,\div$ を有限回使って作った式を算術式という。算術式 $f$ の変数に区間 $X_1,\ldots,X_m$ を代入し、各演算を区間の演算で置き換えて評価したものを $F(X_1,\ldots,X_m)$ と書き、$f$ の自然な区間拡張という。単調な初等関数 $\varphi$（$\sqrt{\ }$, $\exp$ など）は $\varphi([a,b])=[\varphi(a),\varphi(b)]$（減少なら端点を入れ替える）として式に含めてよい。

定理 4.4（区間演算の包含原理）✓

算術式 $f$ の自然な区間拡張 $F(X_1,\ldots,X_m)$ が（$0$ を含む区間での割り算なしに）計算できたとする。このとき、すべての $x_i\in X_i$ について $f(x_1,\ldots,x_m)$ は定義されて $F(X_1,\ldots,X_m)$ に属する。さらに $X_i\subset Y_i$（すべての $i$）で $F(Y_1,\ldots,Y_m)$ も計算できるなら、$F(X)\subset F(Y)$ である（包含単調性）。

証明（定理 4.4）

式の構成に関する帰納法です。変数 $x_i$ と定数では明らかです。$f=g\circ h$ で $g$, $h$ について主張が成り立つとします。$x_i\in X_i$ なら $g(x)\in G(X)$, $h(x)\in H(X)$ で、定義 4.1 により $g(x)\circ h(x)\in G(X)\circ H(X)=F(X)$ です。割り算では $H(X)\not\ni0$ から $h(x)\ne0$ なので $f(x)$ は定義されます。単調な初等関数でも同様です。包含単調性も同じ帰納法で、$A\subset A'$, $B\subset B'$ なら定義から $A\circ B\subset A'\circ B'$ となることを使います。

∎

包含原理は値域を含むことしか言いません。同じ変数が式に何度も現れると、区間演算はそれらを独立に動かすので過大評価が生じます（依存性の問題）。たとえば $X=[0,1]$ で $X-X=[-1,1]$ ですが、$x-x$ の値域は $\{0\}$ です。

依存性の問題。X = [0,1] のとき、x − x の値域は {0} の 1 点だが、区間演算は 2 つの X を独立に動かすので X − X = [−1,1] になる。

4.2. 計算機の上の区間演算

端点を浮動小数に保つため、下端は下向き、上端は上向きに丸めます。

定義 4.5（外向き丸めの区間演算）↗

端点が $\mathbb{F}$ に属する区間 $X,Y$ に対し、$X\circ Y=[\ell,r]$ のとき $X\circledast Y=[\nabla(\ell),\Delta(r)]$ と定める。

命題 4.6（外向き丸めの包含）✓

$X\circ Y\subset X\circledast Y$ である。したがって 定理 4.4 は、各演算を $\circledast$ に置き換えても成り立つ。また $\nabla(\min_kp_k)=\min_k\nabla(p_k)$, $\Delta(\max_kp_k)=\max_k\Delta(p_k)$ なので、命題 4.2 の各端点を丸めてから最小・最大をとってよい。

証明（命題 4.6）

定義から $\nabla(\ell)\le\ell$, $\Delta(r)\ge r$ です。帰納法の各段で、真の区間 $G(X)\circ H(X)$ は、それを含む機械区間どうしの $\circ$ に含まれ（包含単調性）、さらにその $\circledast$ に含まれます。最後の等式は $\nabla$, $\Delta$ が単調非減少であることから従います。

∎

4.3. 丸めの向きを指定できないとき

C 言語には丸めの向きを切り替える fesetround がありますが、コンパイラの最適化と組み合わせると期待どおりに働かないことがあり、Python や JavaScript には切り替える手段がそもそもありません。使えるのは最近点丸めだけです。その場合は、最近点丸めの結果を隣の浮動小数へ 1 つ外側に動かします。次の補題は、これで真の値を取りこぼさないことの根拠です。「忠実な丸め」とは、真の値のすぐ下か、すぐ上の浮動小数のどちらかを返す丸めのことです。浮動小数 $f$ の 1 つ下・1 つ上の浮動小数を $\mathrm{pred}(f)$, $\mathrm{succ}(f)$ と書きます（Python の math.nextafter(f, -math.inf)、math.nextafter(f, math.inf)）。

補題 4.7（1 つ外側に動かせば含む）✓

$s$ を $|s|\le x_{\max}$ の実数とし、$f$ を $s$ の忠実な丸め（$f\in\{\nabla(s),\Delta(s)\}$）とする。このとき $\mathrm{pred}(f)\le s\le\mathrm{succ}(f)$ である。とくに最近点丸めの結果 $f=\mathrm{fl}(s)$ について成り立つ。

証明（補題 4.7）

$s\in\mathbb{F}$ なら $f=s$ で明らかです。$s\notin\mathbb{F}$ なら $\nabla(s)\lt s\lt\Delta(s)$ で、$\nabla(s)$ は $s$ 以下の最大の浮動小数、$\Delta(s)$ は $s$ 以上の最小の浮動小数なので、その間に浮動小数はなく $\Delta(s)=\mathrm{succ}(\nabla(s))$ です。$f=\nabla(s)$ なら $s\lt\Delta(s)=\mathrm{succ}(f)$ かつ $\mathrm{pred}(f)\lt f\lt s$、$f=\Delta(s)$ なら $\mathrm{pred}(f)=\nabla(s)\lt s\lt f$ です。

∎

nextafter による外向き丸め。最近点丸めの結果 fl(s) は真の値 s の隣にある浮動小数のどちらかなので、1 つ下と 1 つ上の浮動小数ではさめば必ず s を含む。

これで、四則演算については「最近点丸めの結果を pred と succ で外に広げる」区間演算が 命題 4.6 と同じ包含を満たします。幅は真の外向き丸めより最大で 1 ulp（隣り合う浮動小数の間隔）ずつ広くなりますが、正しさは失われません。注意すべき点を挙げます。

• 正しく丸められる演算に限る。 補題 4.7 の仮定は「結果が忠実な丸め」であることです。IEEE 754 はこれを四則・平方根・fma について要求しています。math.exp, math.sin, math.pow などの数学関数については、IEEE 754-2019 は正しい丸めを推奨するだけで要求しておらず、C や Python の数学ライブラリも誤差の上限を保証しません。したがって math.exp(x) を 1 つずつ外に広げても、$e^x$ を含む保証はありません。初等関数は、第 2 章のように四則だけで Taylor 展開と剰余を区間で評価するか、誤差が証明されたライブラリを使います。
• 定数の変換も丸めである。 例 2.5 のとおり、ソースに書いた 0.1 は $1/10$ ではありません。I(pred(0.1), succ(0.1)) とするか、整数の割り算 I(1)/I(10) で区間を作ります。
• オーバーフローと NaN。 途中で inf や nan が出たら、その結果を区間として信用せず検証失敗とします。

4.4. 球算術

区間を中心と半径で表す方法を球算術（中心半径演算）といいます。Johansson のライブラリ Arb（現在は FLINT に統合）は、中心を任意精度の浮動小数、半径を低精度の浮動小数で持つ球算術です。中心の計算は普通の浮動小数演算 1 回で済み、半径の計算は安価な上向き丸めの演算で足ります。

定義 4.8（球と球の演算）↗

$m\in\mathbb{F}$, $r\in\mathbb{F}$, $r\ge0$ に対し $B(m,r)=[m-r,m+r]$ と書く。2 つの球 $B_i=B(m_i,r_i)$ に対し

$$\begin{aligned} B_1\oplus B_2&=B(m,r),\quad m=\mathrm{fl}(m_1+m_2),\ r\ge r_1+r_2+u|m|,\\ B_1\otimes B_2&=B(m,r),\quad m=\mathrm{fl}(m_1m_2),\ r\ge|m_1|r_2+|m_2|r_1+r_1r_2+u|m| \end{aligned}$$

と定める。$r$ は右辺を上向き丸め（またはその代わりの succ）で計算した浮動小数とする。

定理 4.9（球算術の包含）✓

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

証明（定理 4.9）

定理 2.4 の 1 の後半から $|m_1+m_2-m|\le u|m|$ なので、

$$|x+y-m|\le|x-m_1|+|y-m_2|+|m_1+m_2-m|\le r_1+r_2+u|m|\le r.$$

積では $xy-m_1m_2=m_1(y-m_2)+m_2(x-m_1)+(x-m_1)(y-m_2)$ と $|m_1m_2-m|\le u|m|$ から

$$|xy-m|\le|m_1|r_2+|m_2|r_1+r_1r_2+u|m|\le r$$

です。

∎

球算術は区間の幅が小さいときに効率的ですが、幅が大きいと粗くなります。$[0,2]\times[0,2]=[0,4]$ ですが、球では $B(1,1)\otimes B(1,1)=B(1,3)=[-2,4]$ です。一方、高精度の計算では中心の演算だけが重く、半径の計算はほぼ無料なので、Arb は数千桁の検証計算にも使われています。

5. 区間 Newton 法

連続関数の零点の存在は中間値の定理で、一意性は単調性で示せます。区間 Newton 法は、その 2 つを同時に、しかも計算機で確かめられる形にしたものです。

定義 5.1（区間 Newton 作用素）↗

$X=[a,b]$ とし、$f:X\to\mathbb{R}$ は $X$ で連続、$(a,b)$ で微分可能とする。$m\in X$、区間 $Y\ni f(m)$、区間 $D$ が $f'((a,b))\subset D$ と $0\notin D$ を満たすとき、

$$N=m-Y\div D$$

を区間 Newton 作用素の値という（$D$ は $f'$ の区間拡張を $X$ で評価したもの、$Y$ は $f$ の区間拡張を $[m,m]$ で評価したものを想定する）。

定理 5.2（区間 Newton 法）✓

定義 5.1 の状況で、次が成り立つ。

1. $f$ の $X$ 内の零点はすべて $N$ に属する。
2. $N\subset X$ ならば、$f$ は $X$ にちょうど 1 つの零点をもつ。
3. $N\cap X=\emptyset$ ならば、$f$ は $X$ に零点をもたない。

証明（定理 5.2）

$D$ は $0$ を含まない区間なので、すべて正かすべて負です。$f$ を $-f$ に替えると $Y$, $D$ が $-Y$, $-D$ に替わり $N$ は変わらないので、$D\subset(0,\infty)$ と仮定してかまいません。

1. $x^\ast\in X$ を零点とします。$x^\ast=m$ なら $0=f(m)\in Y$ なので $m=m-0/d\in N$（$d\in D$）です。$x^\ast\ne m$ なら、平均値の定理（定理 3.3［平均値の定理とテイラーの定理］）により $x^\ast$ と $m$ の間の $\xi$ で $0=f(x^\ast)=f(m)+f'(\xi)(x^\ast-m)$ です。$f'(\xi)\in D$ は $0$ でないので $x^\ast=m-f(m)/f'(\xi)\in m-Y\div D=N$ です。

2. 一意性。 $(a,b)$ で $f'\gt0$ なので $f$ は $X$ で狭義単調増加です（系 3.4［平均値の定理とテイラーの定理］）。零点は高々 1 つです。

   存在。 $f(m)=0$ なら終わりです。$f(m)\gt0$ とします。$m=a$ だとすると、任意の $d\in D$ について $m-f(m)/d\lt a$ ですが、これは $N\subset X$ に反します。よって $a\lt m$ で、平均値の定理から $f(m)-f(a)=f'(\xi)(m-a)$ となる $\xi\in(a,m)$ があります。$m-f(m)/f'(\xi)\in N\subset X$ なので $m-f(m)/f'(\xi)\ge a$、すなわち $f(m)\le f'(\xi)(m-a)=f(m)-f(a)$ で、$f(a)\le0$ です。$f(a)\le0\lt f(m)$ と中間値の定理（定理 5.7［極限と連続性］）から $[a,m]$ に零点があります。$f(m)\lt0$ のときは同じ議論を $b$ の側で行い、$f(b)\ge0$ を得ます。

3. 1 から、零点があれば $N\cap X$ に属するはずです。

∎

3 つの主張は役割が違います。1 は解を落とさないこと、2 は存在と一意性の証明、3 は解がない領域の除外です。$X_{k+1}=N(X_k)\cap X_k$ と反復すると、零点を含む区間が単調に縮みます。$f$ が $C^2$ で零点が単純なら、幅は 2 次収束で縮みます（Moore–Kearfott–Cloud の区間 Newton 法の章）。多変数では、割り算の代わりに近似逆行列を使う Krawczyk 作用素が同じ役割を果たします。次にそれを証明します。考え方は、方程式 $f(x)=0$ を「$g(x)=x$ となる点（不動点）を探す問題」に書き直し、$g$ が直方体をそれ自身に写すことを区間演算で確かめる、というものです。存在の部分は平均値の定理ではなく、Brouwer の不動点定理（定理 6.1［胞体ホモロジーと写像度］）を使います。

定理 5.3（Krawczyk の定理）✓

$f:U\to\mathbb{R}^n$ は開集合 $U\subset\mathbb{R}^n$ で $C^1$ 級、$X=\prod_i[a_i,b_i]\subset U$（すべて $a_i\lt b_i$）は直方体、$\tilde x\in X$、$R\in\mathbb{R}^{n\times n}$ とする。区間行列 $\mathbf{J}$（各成分が区間の行列）が、すべての $\xi\in X$ について $f'(\xi)\in\mathbf{J}$（成分ごとに属する）を満たすとする。区間演算で

$$K(X)=\tilde x-Rf(\tilde x)+(I-R\mathbf{J})(X-\tilde x)$$

を計算する（区間行列と区間ベクトルの積は成分ごとの区間演算）。このとき次が成り立つ。

1. $f$ の $X$ 内の零点はすべて $K(X)$ に属する。
2. $K(X)\subset\operatorname{int}X$ ならば、$R$ と $\mathbf{J}$ の任意の行列は正則で、$f$ は $X$ にちょうど 1 つの零点をもつ。

証明（定理 5.3）

$g(x)=x-Rf(x)$ とおきます。$x\in X$ に対し $J_x=\int_0^1f'\bigl(\tilde x+t(x-\tilde x)\bigr)dt$（成分ごとの積分）とすると、$X$ は凸なので被積分関数の各成分は対応する $\mathbf{J}$ の成分の区間に入り、その平均 $J_x$ も入ります。すなわち $J_x\in\mathbf{J}$ です。微積分学の基本定理から $f(x)=f(\tilde x)+J_x(x-\tilde x)$ なので

$$g(x)=\tilde x-Rf(\tilde x)+(I-RJ_x)(x-\tilde x)\in K(X)\qquad(x\in X)\tag{*}$$

です（区間演算の定義から、各行列 $M\in I-R\mathbf{J}$ とベクトル $z\in X-\tilde x$ について $Mz$ は積の区間に入る）。

1. $x^\ast\in X$ が零点なら $g(x^\ast)=x^\ast$ で、$(\ast)$ から $x^\ast=g(x^\ast)\in K(X)$ です。

2. $K(X)\subset\operatorname{int}X$ とします。$Z=X-\tilde x$, $\boldsymbol{b}=-Rf(\tilde x)$ とおくと、$I-R\mathbf{J}$ に属するどの行列 $M$ についても、アフィン写像 $T_M(z)=\boldsymbol{b}+Mz$ は $Z$ を $\operatorname{int}Z$ へ写します。

   $R\mathbf{J}$ の正則性。 $J\in\mathbf{J}$ を任意にとり $M=I-RJ$ とします。$RJ$ が特異なら $My=y$ となる $\boldsymbol{y}\ne\boldsymbol{0}$ があります。$X$ は直方体で $Z$ は $\boldsymbol{0}$ を含む有界閉凸集合、$T_M(Z)\subset Z$ で、直方体は閉球と同相なので、Brouwer の不動点定理（定理 6.1［胞体ホモロジーと写像度］）により不動点 $z_0\in Z$ があります。このとき $T_M(z_0+t\boldsymbol{y})=\boldsymbol{b}+Mz_0+tM\boldsymbol{y}=z_0+t\boldsymbol{y}$ がすべての実数 $t$ で成り立ちます。$Z$ は有界なので、$t_2=\max\{t:z_0+t\boldsymbol{y}\in Z\}$ が存在し、$z_0+t_2\boldsymbol{y}\in\partial Z$ です。ところが $T_M$ はこの点を自分自身へ写し、一方で $T_M(Z)\subset\operatorname{int}Z$ なので、$\partial Z$ の点は像に入りません。矛盾です。よって $RJ$ は正則で、$R$ も $J$ も正則です。

   存在。 $(\ast)$ と $K(X)\subset X$ から $g$ は $X$ をそれ自身へ写す連続写像です。Brouwer の不動点定理により $g(x^\ast)=x^\ast$ となる $x^\ast\in X$ があり、$Rf(x^\ast)=\boldsymbol{0}$ です。$R$ は正則なので $f(x^\ast)=\boldsymbol{0}$ です。

   一意性。 $x_1\ne x_2\in X$ がともに零点なら、$f(x_1)-f(x_2)=J(x_1-x_2)$（$J=\int_0^1f'(x_2+t(x_1-x_2))dt\in\mathbf{J}$）より $J(x_1-x_2)=\boldsymbol{0}$ ですが、$J$ は正則なので $x_1=x_2$ です。

∎

1 変数のとき、$R=1/f'(m)$ ととれば $K(X)$ は区間 Newton 作用素に近い形になります。区間 Newton 法では存在が中間値の定理で出ますが、多変数では Brouwer の定理が必要なのが違いです。$R$ には $f'(\tilde x)$ の近似逆行列を使い、$\mathbf{J}$ は $f'$ の区間拡張を $X$ で評価して作ります。実際の検証では $I-R\mathbf{J}$ の成分が小さくなるので、$K(X)$ は $X$ の内側に収まりやすくなります。

例 5.4（√2 の包含）

$f(x)=x^2-2$, $X=[1,2]$ とし、nextafter による外向き丸めの区間で 定理 5.2 を反復します。区間の四則はすべて 補題 4.7 の方法で外に広げます。

import math
from fractions import Fraction

INF = math.inf
def down(x): return math.nextafter(x, -INF)   # 1 つ下の浮動小数
def up(x):   return math.nextafter(x, INF)    # 1 つ上の浮動小数

class I:
    """閉区間 [lo, hi]。+ - * / は最近点丸めの結果を 1 つずつ外へ広げる。"""
    def __init__(self, lo, hi=None):
        self.lo, self.hi = lo, (lo if hi is None else hi)
        assert self.lo <= self.hi
    def __add__(self, o):
        return I(down(self.lo + o.lo), up(self.hi + o.hi))
    def __sub__(self, o):
        return I(down(self.lo - o.hi), up(self.hi - o.lo))
    def __mul__(self, o):
        p = [self.lo * o.lo, self.lo * o.hi, self.hi * o.lo, self.hi * o.hi]
        return I(down(min(p)), up(max(p)))
    def __truediv__(self, o):
        if o.lo <= 0 <= o.hi:
            raise ZeroDivisionError("0 を含む区間で割った")
        q = [self.lo / o.lo, self.lo / o.hi, self.hi / o.lo, self.hi / o.hi]
        return I(down(min(q)), up(max(q)))
    def __repr__(self):
        return f"[{self.lo!r}, {self.hi!r}]"

def newton_step(f, df, X):
    m = min(max(X.lo + (X.hi - X.lo) / 2, X.lo), X.hi)   # X に属する浮動小数
    N = I(m) - f(I(m)) / df(X)
    lo, hi = max(N.lo, X.lo), min(N.hi, X.hi)
    return N, (I(lo, hi) if lo <= hi else None)

f  = lambda X: X * X - I(2.0)
df = lambda X: I(2.0) * X
X = I(1.0, 2.0)
for k in range(5):
    N, Xn = newton_step(f, df, X)
    print(k, N, X.lo <= N.lo and N.hi <= X.hi, f"{N.hi - N.lo:.3e}")
    X = Xn
print(Fraction(X.lo) ** 2 < 2 < Fraction(X.hi) ** 2)   # 有理数で確認

出力は次のとおりです（3 列目は $N\subset X$ の判定）。

0 [1.3749999999999996, 1.4375000000000004] True 6.250e-02
1 [1.4140624999999998, 1.414417613636364] True 3.551e-04
2 [1.414213559294524, 1.4142135659471793] True 6.653e-09
3 [1.4142135623730947, 1.4142135623730956] True 8.882e-16
4 [1.4142135623730947, 1.4142135623730954] True 6.661e-16
True

最初の反復で $N=[1.375\ldots,1.4375\ldots]\subset[1,2]$ となった時点で、定理 5.2 の 2 により $[1,2]$ に $x^2=2$ の解がちょうど 1 つあることが証明されています。その後の反復は幅を縮めるだけで、幅は $6\times10^{-2}$, $4\times10^{-4}$, $7\times10^{-9}$ と 2 次収束の様子を示し、隣り合う浮動小数 3〜4 個分で止まります。ここで使った演算は四則だけなので、補題 4.7 の仮定が満たされています。

6. 演習

演習 6.1易

（Sterbenz の補題）$x,y\in\mathbb{F}$ が正の正規化数で $y/2\le x\le2y$ を満たすなら、$x-y\in\mathbb{F}$（引き算は丸めなしで厳密）であることを示せ。

解答

$x\ge y$ としてかまいません（逆なら $y-x$ について示し、$\mathbb{F}$ が符号について対称であることを使います）。$y=M\cdot2^{e-p+1}$（$2^{p-1}\le M\lt2^p$）と書きます。$x\ge y$ より $x$ の指数は $e$ 以上なので、$x$ も $2^{e-p+1}$ の整数倍です。よって $x-y=k\cdot2^{e-p+1}$（$k$ は整数）で、$0\le x-y\le y$ から $0\le k\le M\lt2^p$ です。$k\lt2^p$ の整数に $2^{e-p+1}$（$e\ge e_{\min}$）を掛けた数は、$k$ を $2^{p-1}\le k'\lt2^p$ となるまで 2 倍して指数を下げれば正規化数として、指数が $e_{\min}$ まで下がれば非正規化数として表せるので、$\mathbb{F}$ に属します。

演習 6.2標準

$(m+n)u\lt1$ のとき $(1+\gamma_m)(1+\gamma_n)\le1+\gamma_{m+n}$ を示し、これから $(1+\theta_m)(1+\theta_n)=1+\theta_{m+n}$（$|\theta_k|\le\gamma_k$）と書けることを説明せよ。

解答

$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)=\frac1{(1-mu)(1-nu)}\le\frac1{1-(m+n)u}=1+\gamma_{m+n}.$$

$\theta=(1+\theta_m)(1+\theta_n)-1$ とすると、上からは $\theta\le(1+\gamma_m)(1+\gamma_n)-1\le\gamma_{m+n}$、下からは $\theta\ge(1-\gamma_m)(1-\gamma_n)-1=-\gamma_m-\gamma_n+\gamma_m\gamma_n\ge-\gamma_{m+n}$ です（$\gamma_m+\gamma_n\le\gamma_m+\gamma_n+\gamma_m\gamma_n\le\gamma_{m+n}$）。よって $|\theta|\le\gamma_{m+n}$ です。誤差解析で「$\theta$ の積をまとめる」操作は、この不等式で正当化されます。

演習 6.3標準

$X=[0,1]$ で、同じ関数 $x(1-x)$ を 3 通りの式 $x\times(1-x)$, $x-x\times x$, $\tfrac14-(x-\tfrac12)^2$ で書き、それぞれ区間で評価せよ。ただし最後の式の 2 乗は (a) $Z\times Z$ と (b) 平方 $\{z^2:z\in Z\}$ の 2 通りで計算せよ。真の値域 $[0,\tfrac14]$ と比べよ。

解答

$X\times(1-X)=[0,1]\times[0,1]=[0,1]$。$X-X\times X=[0,1]-[0,1]=[-1,1]$。$Z=X-\tfrac12=[-\tfrac12,\tfrac12]$ について、(a) $Z\times Z=[-\tfrac14,\tfrac14]$ なので $\tfrac14-Z\times Z=[0,\tfrac12]$、(b) 平方は $[0,\tfrac14]$ なので $\tfrac14-[0,\tfrac14]=[0,\tfrac14]$ で真の値域に一致します。どれも 定理 4.4 により値域を含みますが、変数の出現回数が少ない式ほど過大評価が小さく、変数が 1 回しか現れない式 (b) では値域そのものになります（各演算の結果が厳密な像なので、帰納法で等号が保たれるため）。

演習 6.4標準

$f(x)=x^2-2$, $X=[2,3]$, $m=2.5$ で区間 Newton 作用素を計算し、$X$ に零点がないことを示せ。

解答

$f(m)=4.25$, $f'(X)=2X=[4,6]$ です。$4.25\div[4,6]=[4.25/6,4.25/4]=[0.7083\ldots,1.0625]$ なので、$N=2.5-[0.7083\ldots,1.0625]=[1.4375,1.7916\ldots]$ で、$N\cap[2,3]=\emptyset$ です。定理 5.2 の 3 により $[2,3]$ に零点はありません。計算機で行うときは、$4.25/6$ の丸めを外向きにすれば同じ結論が厳密に得られます。

参考文献

• N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002 — 第 2 章（浮動小数）、第 3 章（$\gamma_n$ と内積・和の誤差解析）。
• IEEE Computer Society, IEEE Standard for Floating-Point Arithmetic (IEEE Std 754-2019), 2019 — 第 4 章（丸めの属性）、第 5 章（正しく丸められる演算）、第 9 章（推奨される演算）。
• R. E. Moore, R. B. Kearfott, M. J. Cloud, Introduction to Interval Analysis, SIAM, 2009 — 区間の演算、包含単調性、区間 Newton 法と Krawczyk 法。
• 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
• S. M. Rump, “Verification methods: Rigorous results using floating-point arithmetic”, Acta Numerica 19 (2010), 287–449. doi:10.1017/S096249291000005X