# 浮動小数と区間演算：IEEE 754 の丸めモデルから区間 Newton 法まで
> IEEE 754 の浮動小数の丸めを u = 2⁻⁵³ の誤差モデルにまとめ、Higham の γₙ による内積の誤差評価、区間演算と球算術の包含原理、nextafter による外向き丸め、区間 Newton 法による解の存在と一意性の証明までを扱う。
> https://rikai.mugen-giken.com/mathematics/computer-assisted-proofs/floating-point-interval-arithmetic
## 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 つあることを結論します。計算機の出力が、そのまま存在と一意性の証明になります。
## 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. 浮動小数と丸め
整数 $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$ ずつ動くので、次がわかります。
$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^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}$ の正規化数とつながって同じ間隔で並びます。
$|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)$ が返ります。この要求が、以下のすべての誤差評価の出発点です。
$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}$）。
対称性から $x\gt0$ とし、$2^e\le x\lt2^{e+1}$（$e_{\min}\le e\le e_{\max}$）とします。