浮動小数点演算 | 数値解析 I | 天游ヶ丘数理学院

1 浮動小数点数系

定義 1.1. $\beta\ge2$ 、 $p\ge1$ 、 $e_{\min}\le e_{\max}$ を整数とし、整数 $e$ に対して $s_e:=\beta^{e+1-p}$ と置く。次の三つの集合の和集合 $F=F(\beta,p,e_{\min},e_{\max})$ を、基数 $\beta$ 、仮数桁数 $p$ 、指数範囲 $[e_{\min},e_{\max}]$ の浮動小数点数系 (floating-point number system) といい、 $F$ の元を浮動小数点数 (floating-point number) という。

$\{0\}$ 。
$\{\sigma m s_e\mid \sigma\in\{1,-1\},\ e,m\in\Z,\ e_{\min}\le e\le e_{\max},\ \beta^{p-1}\le m\le\beta^p-1\}$ 。この集合の元を正規化数 (normalized number) という。
$\{\sigma m s_{e_{\min}}\mid \sigma\in\{1,-1\},\ m\in\Z,\ 1\le m\le\beta^{p-1}-1\}$ 。この集合の元を非正規数 (subnormal number) という。 $p=1$ のとき、この集合は空である。

$N_{\max}:=(\beta^p-1)s_{e_{\max}}=\beta^{e_{\max}+1}(1-\beta^{-p})$ と置く。 $\beta^{e_{\min}}\le|z|\le N_{\max}$ を満たす実数 $z$ の全体を $F$ の正規範囲 (normal range) という。実数 $z$ が $|z|>N_{\max}$ を満たすとき $z$ はオーバーフロー (overflow) の範囲にあるといい、 $0<|z|<\beta^{e_{\min}}$ を満たすとき $z$ はアンダーフロー (underflow) の範囲にあるという。

補題 1.2. $F=F(\beta,p,e_{\min},e_{\max})$ を浮動小数点数系とする。 $0\le t\le N_{\max}$ を満たす実数 $t$ に対し、 $t<\beta^{e_{\min}}$ のとき $e(t):=e_{\min}$ と置き、 $t\ge\beta^{e_{\min}}$ のとき $\beta^{e(t)}\le t<\beta^{e(t)+1}$ を満たす整数を $e(t)$ と置く。

$e_{\min}\le e\le e_{\max}$ を満たす整数 $e$ に対して、 $F\cap[\beta^e,\beta^{e+1})=\{ms_e\mid m\in\Z,\ \beta^{p-1}\le m\le\beta^p-1\}$ が成り立つ。特に $N_{\max}$ は $F$ の最大元である。
$F\cap[0,\beta^{e_{\min}})=\{ms_{e_{\min}}\mid m\in\Z,\ 0\le m\le\beta^{p-1}-1\}$ が成り立つ。
$F$ の各元 $x\ne0$ に対して、 $x=\sigma ms_e$ 、 $\sigma\in\{1,-1\}$ 、 $e_{\min}\le e\le e_{\max}$ を満たし、さらに $\beta^{p-1}\le m\le\beta^p-1$ であるか $e=e_{\min}$ かつ $1\le m\le\beta^{p-1}-1$ であるような整数の組 $(\sigma,e,m)$ はただ一つ存在する。
$x\in F$ が $0\le x<N_{\max}$ を満たすならば、 $x$ より大きい $F$ の元のうち最小のものは $x+s_{e(x)}$ である。
$0\le t\le N_{\max}$ を満たす実数 $t$ に対して、 $a\le t<a+s_{e(t)}$ を満たす $a\in F$ が存在し、 $a+s_{e(t)}\in F$ であるか $t=a=N_{\max}$ である。 $t\ge\beta^{e_{\min}}$ ならば $a\ge\beta^{e(t)}$ である。

証明. $0\le t\le N_{\max}<\beta^{e_{\max}+1}$ ならば $e_{\min}\le e(t)\le e_{\max}$ である。 $\beta^{p-1}\le m\le\beta^p-1$ ならば $\beta^e=\beta^{p-1}s_e\le ms_e\le(\beta^p-1)s_e<\beta^{e+1}$ であるから、正規化数 $\sigma ms_e$ の絶対値は $[\beta^e,\beta^{e+1})$ に属する。これらの区間は $e$ ごとに交わらない。正の非正規数 $ms_{e_{\min}}$ は $ms_{e_{\min}}\le(\beta^{p-1}-1)s_{e_{\min}}<\beta^{e_{\min}}$ を満たす。したがって $F\cap[\beta^e,\beta^{e+1})$ の元は指数 $e$ の正の正規化数であり、(1)の等式が成り立つ。 $F$ の正の元は $[0,\beta^{e_{\min}})$ か、ある $e_{\min}\le e\le e_{\max}$ に対する $[\beta^e,\beta^{e+1})$ に属するので、 $F$ の最大元は $F\cap[\beta^{e_{\max}},\beta^{e_{\max}+1})$ の最大元 $(\beta^p-1)s_{e_{\max}}=N_{\max}$ である。正規化数の絶対値は $\beta^{e_{\min}}$ 以上であるから、 $F\cap[0,\beta^{e_{\min}})$ は $0$ と正の非正規数からなり、(2)が成り立つ。

(3)を示す。 $\sigma$ は $x$ の符号である。 $|x|\ge\beta^{e_{\min}}$ ならば、 $|x|$ は非正規数の絶対値ではなく、 $|x|\in[\beta^e,\beta^{e+1})$ から $e=e(|x|)$ 、 $m=|x|/s_e$ が定まる。 $|x|<\beta^{e_{\min}}$ ならば、 $|x|$ は正規化数の絶対値ではないので $e=e_{\min}$ 、 $m=|x|/s_{e_{\min}}$ である。

(4)を示す。 $x\ge\beta^{e_{\min}}$ とし、 $e:=e(x)$ と置く。(1)により $x=ms_e$ 、 $\beta^{p-1}\le m\le\beta^p-1$ であり、 $(x,\beta^{e+1})$ に属する $F$ の元は $m<m'\le\beta^p-1$ を満たす $m's_e$ である。 $m\le\beta^p-2$ ならば、その最小元は $(m+1)s_e=x+s_e$ である。 $m=\beta^p-1$ ならば $(x,\beta^{e+1})$ に $F$ の元は無く、 $x<N_{\max}$ であるから $e<e_{\max}$ であり、 $x+s_e=\beta^{e+1}=\beta^{p-1}s_{e+1}\in F$ である。 $\beta^{e+1}$ 以上の元は $x+s_e$ 以上であるから、どちらの場合も最小元は $x+s_e$ である。 $0\le x<\beta^{e_{\min}}$ の場合は、(2)により $x=ms_{e_{\min}}$ 、 $0\le m\le\beta^{p-1}-1$ であり、 $(m+1)s_{e_{\min}}$ は $m+1\le\beta^{p-1}-1$ ならば非正規数、 $m+1=\beta^{p-1}$ ならば正規化数 $\beta^{e_{\min}}$ である。(2)により $[0,\beta^{e_{\min}})$ に属する $F$ の元は $s_{e_{\min}}$ の整数倍であり、正規化数は $\beta^{e_{\min}}$ 以上であるから、 $(x,(m+1)s_{e_{\min}})$ に $F$ の元は無く、最小元は $x+s_{e_{\min}}$ である。

(5)を示す。 $a:=\max\{y\in F\mid y\le t\}$ と置く。 $0\in F$ であるから $a$ は存在する。 $t\ge\beta^{e_{\min}}$ ならば $\beta^{e(t)}=\beta^{p-1}s_{e(t)}\in F$ かつ $\beta^{e(t)}\le t$ であるから $\beta^{e(t)}\le a\le t<\beta^{e(t)+1}$ であり、 $e(a)=e(t)$ である。 $t<\beta^{e_{\min}}$ ならば $0\le a<\beta^{e_{\min}}$ であり、 $e(a)=e_{\min}=e(t)$ である。 $a<N_{\max}$ ならば、(4)により $a+s_{e(t)}$ は $a$ より大きい最小の $F$ の元であり、 $a$ の最大性から $a+s_{e(t)}>t$ である。 $a=N_{\max}$ ならば $N_{\max}=a\le t\le N_{\max}$ より $t=a$ である。▨

補題 1.3. $F=F(\beta,p,e_{\min},e_{\max})$ を浮動小数点数系とし、 $0\le t\le N_{\max}$ を満たす実数 $t$ に対して $e(t)$ を補題 1.2の整数とする。

整数 $e$ が $e_{\min}\le e\le e_{\max}$ を満たし、整数 $k$ が $|k|\le\beta^p$ と $|ks_e|\le N_{\max}$ を満たすならば、 $ks_e\in F$ である。
$x\in F$ が $x\ne0$ を満たすならば、 $1\le|m|\le\beta^p-1$ を満たす整数 $m$ が存在して $x=ms_{e(|x|)}$ が成り立つ。さらに、 $e_{\min}\le e\le e(|x|)$ を満たす任意の整数 $e$ について、 $x$ は $s_e$ の整数倍である。
$0\le t\le t'\le N_{\max}$ ならば $e(t)\le e(t')$ である。

証明. 演習とする（問題 6.2）。▨

定義 1.4. $F=F(\beta,p,e_{\min},e_{\max})$ を浮動小数点数系とする。

$u:=\beta^{1-p}/2$ を $F$ の単位丸め誤差 (unit roundoff) という。
$e_{\min}\le0<e_{\max}$ のとき、 $1=\beta^{p-1}s_0$ と $N_{\max}\ge\beta^{e_{\max}}>1$ は $F$ の元である。 $1$ より大きい $F$ の元のうち最小のものから $1$ を引いた値 $\epsilon_{\mathrm{mach}}$ を $F$ の機械イプシロン (machine epsilon) という。

系 1.5. $F=F(\beta,p,e_{\min},e_{\max})$ を $e_{\min}\le0<e_{\max}$ を満たす浮動小数点数系とする。 $F$ の機械イプシロン $\epsilon_{\mathrm{mach}}$ と単位丸め誤差 $u$ は $\epsilon_{\mathrm{mach}}=\beta^{1-p}=2u$ を満たす。

証明. $e(1)=0$ かつ $1<N_{\max}$ であるから、補題 1.2 (4)により $1$ より大きい最小の $F$ の元は $1+s_0=1+\beta^{1-p}$ である。▨

例 1.6. $F=F(2,3,-2,2)$ とする。 $s_e=2^{e-2}$ であり、 $s_{-2}=1/16$ 、 $s_{-1}=1/8$ 、 $s_0=1/4$ 、 $s_1=1/2$ 、 $s_2=1$ である。補題 1.2 (1)と補題 1.2 (2)により、 $F$ の正の元は次の $23$ 個であり、 $|F|=2\cdot23+1=47$ である。

正の非正規数は $1/16,\ 1/8,\ 3/16$ である。
$[1/4,1/2)$ に属する元は $1/4,\ 5/16,\ 3/8,\ 7/16$ である。
$[1/2,1)$ に属する元は $1/2,\ 5/8,\ 3/4,\ 7/8$ である。
$[1,2)$ に属する元は $1,\ 5/4,\ 3/2,\ 7/4$ である。
$[2,4)$ に属する元は $2,\ 5/2,\ 3,\ 7/2$ である。
$[4,8)$ に属する元は $4,\ 5,\ 6,\ 7$ である。

$N_{\max}=7$ 、正規範囲の下端は $\beta^{e_{\min}}=1/4$ 、単位丸め誤差は $u=1/8$ 、機械イプシロンは $1/4$ である。正規化数 $ms_e$ は二進 $3$ 桁の仮数による指数表示 $(m/4)\cdot2^e$ に対応し、たとえば $5/4=5\cdot2^{-2}=(1.01)_2\times2^0$ である。隣接する元の間隔は $[0,1/2)$ で $1/16$ であり、 $e\ge-1$ の帯 $[2^e,2^{e+1})$ で $s_e$ である。帯 $[2^e,2^{e+1})$ の元 $x$ に対する間隔の比は $1/8<s_e/x\le1/4$ を満たし、非正規数 $x=1/16$ では $s_{-2}/x=1$ である。

例 1.7. IEEE 754 の binary64 形式の有限数の全体は $F(2,53,-1022,1023)$ である。 $N_{\max}=(2^{53}-1)2^{971}=(2-2^{-52})2^{1023}\approx1.7976931348623157\times10^{308}$ 、正規範囲の下端は $2^{-1022}\approx2.2250738585072014\times10^{-308}$ 、最小の正の非正規数は $s_{-1022}=2^{-1074}\approx4.94\times10^{-324}$ である。系 1.5により $\epsilon_{\mathrm{mach}}=2^{-52}\approx2.220446049250313\times10^{-16}$ 、 $u=2^{-53}\approx1.1102230246251565\times10^{-16}$ である。 CPython の sys.float_info は max・min・epsilon にそれぞれ $N_{\max}$ ・ $2^{-1022}$ ・ $2^{-52}$ を与える。同じ sys.float_info の min_exp=-1021 と max_exp=1024 は仮数を $[1/2,1)$ に置く C の規約による指数であり、 $e_{\min}=-1022$ 、 $e_{\max}=1023$ とはそれぞれ $1$ だけ異なる。

2 最近接丸め

定義 2.1. $F=F(\beta,p,e_{\min},e_{\max})$ を浮動小数点数系とし、 $z$ を $|z|\le N_{\max}$ を満たす実数とする。 $x\in F$ が任意の $y\in F$ に対して $|x-z|\le|y-z|$ を満たすとき、 $x$ を $z$ の最近接点 (nearest point) という。 $F$ は空でない有限集合であるから、 $z$ の最近接点は存在する。写像 $\operatorname{fl}\colon[-N_{\max},N_{\max}]\to F$ であって、各 $z$ に対して $\operatorname{fl}(z)$ が $z$ の最近接点であるものを、 $F$ の最近接丸め (rounding to nearest) という。

定理 2.2. $F=F(\beta,p,e_{\min},e_{\max})$ を浮動小数点数系、 $u$ をその単位丸め誤差、 $\operatorname{fl}$ を $F$ の最近接丸めとし、 $z$ を $|z|\le N_{\max}$ を満たす実数とする。

$z\in F$ ならば $\operatorname{fl}(z)=z$ である。特に $\operatorname{fl}(0)=0$ である。
$z$ が正規範囲にあるならば、 $|\delta|\le u$ を満たす実数 $\delta$ が存在して $\operatorname{fl}(z)=z(1+\delta)$ が成り立つ。
$|z|<\beta^{e_{\min}}$ ならば $|\operatorname{fl}(z)-z|\le u\beta^{e_{\min}}=s_{e_{\min}}/2$ が成り立つ。

証明. $z\in F$ ならば、 $z$ との距離が $0$ である $F$ の元は $z$ だけであるから、(1)が成り立つ。

$t:=|z|$ と置く。 $F=-F$ であるから $|\operatorname{fl}(z)-z|=\min_{y\in F}|y-z|=\min_{y\in F}|y-t|$ である。補題 1.2 (5)により、 $a\le t<a+s_{e(t)}$ を満たす $a\in F$ があり、 $a+s_{e(t)}\in F$ であるか $t=a$ である。 $a+s_{e(t)}\in F$ ならば $\min_{y\in F}|y-t|\le\min\{t-a,\ a+s_{e(t)}-t\}\le s_{e(t)}/2$ であり、 $t=a$ ならば $\min_{y\in F}|y-t|=0$ である。いずれの場合も

|\operatorname{fl}(z)-z|\le\frac{s_{e(t)}}{2}=u\beta^{e(t)}

が成り立つ。 $z$ が正規範囲にあるならば $\beta^{e(t)}\le t=|z|$ であるから $|\operatorname{fl}(z)-z|\le u|z|$ であり、 $\delta:=(\operatorname{fl}(z)-z)/z$ と置けば(2)が成り立つ。 $|z|<\beta^{e_{\min}}$ ならば $e(t)=e_{\min}$ であり、(3)が成り立つ。▨

補題 2.3. $F=F(\beta,p,e_{\min},e_{\max})$ を浮動小数点数系とする。 $x\in F$ に対し、 $x=0$ のとき $m(x):=0$ と置き、 $x\ne0$ のとき補題 1.2 (3)の整数 $m$ を $m(x)$ と置く。

$\beta$ が奇数であるか $p\ge2$ であるとする。 $|z|\le N_{\max}$ を満たす実数 $z$ が二つの最近接点を持つならば、そのうち $m(x)$ が偶数であるもの $x$ はちょうど一つである。
$\beta$ が偶数で $p=1$ であり、整数 $e$ が $e_{\min}\le e<e_{\max}$ を満たすとする。 $z:=(\beta-1/2)\beta^e$ の最近接点は $(\beta-1)\beta^e$ と $\beta^{e+1}$ の二つであり、どちらも $m(x)$ は奇数である。

証明.(1)を示す。 $z$ の二つの最近接点を $x<x'$ とする。 $z=(x+x')/2$ であり、 $(x,x')$ に属する $y\in F$ は $|y-z|<|x-z|$ を満たすので、 $(x,x')$ に $F$ の元は無い。 $F=-F$ かつ $m(-y)=m(y)$ であるから、 $z$ を $-z$ に替えて $z\ge0$ としてよい。このとき $x'>0$ であり、 $0\in F$ が $(x,x')$ に属さないので $x\ge0$ である。 $x<x'\le N_{\max}$ であるから、補題 1.2 (4)により $x'=x+s_{e(x)}$ である。 $x<\beta^{e_{\min}}$ ならば、補題 1.2 (2)により $x=ms_{e_{\min}}$ 、 $m=m(x)\le\beta^{p-1}-1$ であり、 $x'=(m+1)s_{e_{\min}}$ は $m(x')=m+1$ を満たす。 $x\ge\beta^{e_{\min}}$ かつ $m(x)\le\beta^p-2$ ならば、補題 1.2 (1)により $m(x')=m(x)+1$ である。これらの場合、 $m(x)$ と $m(x')$ は連続する整数であり、ちょうど一方が偶数である。残る場合は $x\ge\beta^{e_{\min}}$ かつ $m(x)=\beta^p-1$ であり、このとき $x'=\beta^{e(x)+1}=\beta^{p-1}s_{e(x)+1}$ から $m(x')=\beta^{p-1}$ である。 $\beta$ が奇数ならば $\beta^p-1$ は偶数で $\beta^{p-1}$ は奇数である。 $\beta$ が偶数で $p\ge2$ ならば $\beta^p-1$ は奇数で $\beta^{p-1}$ は偶数である。

(2)を示す。 $p=1$ であるから $s_e=\beta^e$ であり、補題 1.2 (1)により $(\beta-1)\beta^e\in F$ で $m((\beta-1)\beta^e)=\beta-1$ である。 $e+1\le e_{\max}$ であるから、補題 1.2 (4)により $(\beta-1)\beta^e$ より大きい最小の $F$ の元は $\beta^{e+1}=1\cdot s_{e+1}$ であり、 $m(\beta^{e+1})=1$ である。 $z$ はこの隣接する二元の中点であるから、最近接点はこの二元である。 $\beta$ が偶数であるから $\beta-1$ は奇数である。▨

定義 2.4. $F=F(\beta,p,e_{\min},e_{\max})$ を、 $\beta$ が奇数であるか $p\ge2$ である浮動小数点数系とし、 $m(x)$ を補題 2.3の整数とする。 $F$ の最近接丸め $\operatorname{fl}$ であって、二つの最近接点を持つ各 $z$ に対して $m(\operatorname{fl}(z))$ が偶数であるものを、 $F$ の最近接偶数丸め (round to nearest, ties to even) という。補題 2.3 (1)により、最近接偶数丸めはただ一つ存在する。

例 2.5. $F=F(2,3,-2,2)$ とし、 $\operatorname{fl}$ を $F$ の最近接偶数丸めとする。 $u=1/8$ 、 $\beta^{e_{\min}}=1/4$ 、 $s_{e_{\min}}=1/16$ である。 $9/8$ の最近接点は $1=4s_0$ と $5/4=5s_0$ であり、 $\operatorname{fl}(9/8)=1$ 、 $\delta=-1/9$ である。 $11/8$ の最近接点は $5/4=5s_0$ と $3/2=6s_0$ であり、 $\operatorname{fl}(11/8)=3/2$ 、 $\delta=1/11$ である。どちらも $|\delta|\le u$ を満たす。 $15/64$ はアンダーフローの範囲にあり、最近接点は $1/4$ だけであって、 $|\operatorname{fl}(15/64)-15/64|=1/64$ である。 $3/32$ の最近接点は $1/16=1\cdot s_{-2}$ と $1/8=2s_{-2}$ であり、 $\operatorname{fl}(3/32)=1/8$ である。絶対誤差は $1/32=u\beta^{e_{\min}}$ であり、定理 2.2 (3)の評価は等号で成り立つ。相対誤差は $1/3>u$ である。 $1/32$ の最近接点は $0$ と $1/16$ であり、 $m(0)=0$ であるから $\operatorname{fl}(1/32)=0$ である。絶対誤差は $1/32=u\beta^{e_{\min}}$ 、相対誤差は $1$ である。

注意 2.6. 最近接丸めは $|z|\le N_{\max}$ の範囲だけで定めた。 $z>N_{\max}$ ならば $F$ の元はすべて $z$ より小さく、 $F$ の元のうち $z$ に最も近いものは $N_{\max}$ であるが、相対誤差 $(z-N_{\max})/z$ は $z\to\infty$ で $1$ に近づき、 $u$ で押さえられない。 binary64 の最近接偶数丸めについて、IEEE 754 は $N_{\max}<|z|<(2-2^{-53})2^{1023}$ の $z$ を $\pm N_{\max}$ へ、 $|z|\ge(2-2^{-53})2^{1023}$ の $z$ を符号付きの無限大へ写す。 $(2-2^{-53})2^{1023}=N_{\max}+2^{970}$ は $N_{\max}$ と $2^{1024}$ の中点である。 CPython の浮動小数点加算で、 $N_{\max}+2^{969}$ の計算値は $N_{\max}$ 、 $N_{\max}+2^{970}$ の計算値は無限大である。

3 四則演算のモデル

定義 3.1. $F=F(\beta,p,e_{\min},e_{\max})$ を浮動小数点数系、 $\operatorname{fl}$ を $F$ の最近接丸めとし、 $\circ\in\{+,-,\times,/\}$ とする。 $x,y\in F$ が、 $\circ=/$ のとき $y\ne0$ を満たし、 $|x\circ y|\le N_{\max}$ を満たすとする。実数 $x\circ y$ を演算の厳密な結果といい、 $\operatorname{fl}(x\circ y)$ を正しく丸めた演算 (correctly rounded operation) の結果という。

系 3.2. $F=F(\beta,p,e_{\min},e_{\max})$ を浮動小数点数系、 $u$ をその単位丸め誤差、 $\operatorname{fl}$ を $F$ の最近接丸めとし、 $\circ\in\{+,-,\times,/\}$ とする。 $x,y\in F$ が、 $\circ=/$ のとき $y\ne0$ を満たし、 $|x\circ y|\le N_{\max}$ を満たすとする。

$x\circ y\in F$ ならば $\operatorname{fl}(x\circ y)=x\circ y$ である。特に $x\circ y=0$ ならば $\operatorname{fl}(x\circ y)=0$ である。
$x\circ y$ が正規範囲にあるならば、 $|\delta|\le u$ を満たす実数 $\delta$ が存在して $\operatorname{fl}(x\circ y)=(x\circ y)(1+\delta)$ が成り立つ。
$0<|x\circ y|<\beta^{e_{\min}}$ ならば $|\operatorname{fl}(x\circ y)-x\circ y|\le u\beta^{e_{\min}}$ が成り立つ。
$\circ\in\{+,-\}$ かつ $|x\circ y|<\beta^{e_{\min}}$ ならば $\operatorname{fl}(x\circ y)=x\circ y$ である。

証明.(1)、(2)、(3)は、定理 2.2を $z=x\circ y$ に適用したものである。(4)を示す。 $x\pm y=0$ ならば、(1)により $\operatorname{fl}(x\pm y)=x\pm y$ である。 $x\pm y\ne0$ とする。 $x$ と $y$ のそれぞれは、 $0$ であるか、補題 1.3 (2)を $e=e_{\min}$ として適用して $s_{e_{\min}}$ の整数倍である。したがって $x\pm y=ks_{e_{\min}}$ を満たす整数 $k$ があり、 $|x\pm y|<\beta^{e_{\min}}=\beta^{p-1}s_{e_{\min}}$ から $|k|<\beta^{p-1}\le\beta^p$ である。 $|ks_{e_{\min}}|=|x\pm y|\le N_{\max}$ であるから、補題 1.3 (1)により $x\pm y\in F$ であり、(1)により $\operatorname{fl}(x\pm y)=x\pm y$ である。▨

系 3.3. $F=F(\beta,p,e_{\min},e_{\max})$ を浮動小数点数系、 $u$ をその単位丸め誤差、 $\operatorname{fl}$ を $F$ の最近接丸めとする。 $y,y'\in F$ が $|y+y'|\le N_{\max}$ を満たすならば、 $|\delta|\le u$ を満たす実数 $\delta$ が存在して $\operatorname{fl}(y+y')=(y+y')(1+\delta)$ が成り立つ。

証明. $y+y'$ が正規範囲にあるならば、系 3.2 (2)により主張が成り立つ。 $|y+y'|<\beta^{e_{\min}}$ ならば、系 3.2 (4)により $\operatorname{fl}(y+y')=y+y'$ であり、 $\delta:=0$ とすればよい。▨

注意 3.4. IEEE 754 は、四則演算の結果を厳密な結果の丸めとして返すことを定め、既定の丸め方向を最近接偶数丸めとする。 binary64 の既定の丸め方向の下では、厳密な結果が $[-N_{\max},N_{\max}]$ にある演算の結果は $\operatorname{fl}(x\circ y)$ であって、系 3.2が $u=2^{-53}$ で適用される。 IEEE 754 は、 $0$ への丸め、 $+\infty$ への丸め、 $-\infty$ への丸めも定める。これらの丸め方向では計算値は補題 1.2 (5)の $\pm a$ または $\pm(a+s_{e(|z|)})$ であり、正規範囲の $z$ について相対誤差は $s_{e(|z|)}/\beta^{e(|z|)}=\beta^{1-p}=2u$ 以下である。中間結果を binary64 より広い形式で保持してから binary64 へ丸め直す実行環境や、非正規数の結果を $0$ へ置き換える実行環境では、計算値が $\operatorname{fl}(x\circ y)$ と一致するとは限らず、系 3.2はそのままでは適用されない。既定の丸め方向（最近接偶数丸め）で実行した CPython の float では、本記事の binary64 の例に現れる演算のうち、厳密な結果が $[-N_{\max},N_{\max}]$ にあるすべての演算について、計算値が $\operatorname{fl}(x\circ y)$ と一致した。

例 3.5. $F=F(2,53,-1022,1023)$ とし、 $\operatorname{fl}$ を最近接偶数丸めとする。 $F$ の元は分母が $2$ の冪である有理数であるから、 $1/10$ と $1/5$ は $F$ に属さない。十進の入力 $0.1$ と $0.2$ の変換値を $\hat x:=\operatorname{fl}(1/10)$ 、 $\hat y:=\operatorname{fl}(1/5)$ とすると、

\hat x=7205759403792794\cdot2^{-56},\qquad \hat y=7205759403792794\cdot2^{-55},

\hat x-\frac1{10}=\frac{2^{-55}}{5},\qquad \hat y-\frac15=\frac{2^{-54}}{5}

であり、入力の相対誤差はどちらも $2^{-54}=u/2$ である。 $\hat x+\hat y=5404319552844595.5\cdot2^{-54}$ は $[1/4,1/2)$ に属し、 $s_{-2}=2^{-54}$ であるから、最近接点は $5404319552844595\cdot2^{-54}$ と $5404319552844596\cdot2^{-54}$ の二つである。最近接偶数丸めにより $\operatorname{fl}(\hat x+\hat y)=5404319552844596\cdot2^{-54}$ であり、演算の誤差は $\operatorname{fl}(\hat x+\hat y)-(\hat x+\hat y)=2^{-55}$ 、その相対誤差は $2^{-55}/(\hat x+\hat y)\approx9.25\times10^{-17}\le u$ である。入力の誤差は $(\hat x+\hat y)-3/10=3\cdot2^{-55}/5$ であり、全体の誤差は

\operatorname{fl}(\hat x+\hat y)-\frac3{10}=\frac{3\cdot2^{-55}}{5}+2^{-55}=\frac{2^{-52}}{5}

である。全体の相対誤差 $2^{-52}\cdot2/3\approx1.48\times10^{-16}$ は $u$ を超え、系 3.2が評価するのは演算の誤差 $2^{-55}$ だけである。一方 $\operatorname{fl}(3/10)=5404319552844595\cdot2^{-54}$ であるから、 $\operatorname{fl}(\hat x+\hat y)=\operatorname{fl}(3/10)+2^{-54}$ であり、二つの計算値は $F$ の隣接する異なる元である。

4 丸め誤差の積

補題 4.1. $u\ge0$ を実数、 $k\ge1$ を $ku<1$ を満たす整数とし、 $\gamma_k:=ku/(1-ku)$ と置く。実数 $\delta_1,\dots,\delta_k$ が $|\delta_i|\le u$ を満たし、 $\varepsilon_1,\dots,\varepsilon_k\in\{1,-1\}$ とする。このとき

\prod_{i=1}^k(1+\delta_i)^{\varepsilon_i}=1+\theta_k,\qquad |\theta_k|\le\gamma_k

を満たす実数 $\theta_k$ が存在する。

証明. $k\ge1$ かつ $ku<1$ であるから $0\le u<1$ である。各 $i$ について $1-u\le1+\delta_i\le1+u$ であり、 $(1+u)(1-u)=1-u^2\le1$ から $1+u\le(1-u)^{-1}$ であるので、 $1-u\le1+\delta_i\le(1-u)^{-1}$ が成り立つ。 $1+\delta_i\ge1-u>0$ であるから、逆数をとって $1-u\le(1+\delta_i)^{-1}\le(1-u)^{-1}$ も成り立つ。したがって $P:=\prod_{i=1}^k(1+\delta_i)^{\varepsilon_i}$ は $(1-u)^k\le P\le(1-u)^{-k}$ を満たす。整数 $j\ge0$ に対して $(1-u)^j\ge1-ju$ が成り立つ。実際、 $j=0$ では等号であり、 $(1-u)^j\ge1-ju$ の両辺に $1-u\ge0$ を掛けると

(1-u)^{j+1}\ge(1-ju)(1-u)=1-(j+1)u+ju^2\ge1-(j+1)u

が得られる。 $j=k$ として $(1-u)^k\ge1-ku>0$ であるから、 $\theta_k:=P-1$ は

-ku\le\theta_k\le\frac{1}{1-ku}-1=\gamma_k

を満たす。 $0<1-ku\le1$ から $ku\le\gamma_k$ であり、 $|\theta_k|\le\gamma_k$ が成り立つ。▨

例 4.2. $F$ を浮動小数点数系、 $u$ をその単位丸め誤差、 $\operatorname{fl}$ を $F$ の最近接丸めとし、 $n\ge2$ を $(n-1)u<1$ を満たす整数、 $x_1,\dots,x_n\in F$ とする。 $\hat p_1:=x_1$ 、 $\hat p_j:=\operatorname{fl}(\hat p_{j-1}x_j)$ （ $j=2,\dots,n$ ）と置き、各 $j$ について厳密な積 $\hat p_{j-1}x_j$ が正規範囲にあると仮定する。系 3.2 (2)により $\hat p_j=\hat p_{j-1}x_j(1+\delta_j)$ 、 $|\delta_j|\le u$ であり、 $j$ に関する帰納法で $\hat p_n=x_1\cdots x_n\prod_{j=2}^n(1+\delta_j)$ が成り立つ。補題 4.1を $k=n-1$ 、 $\varepsilon_j=1$ に適用して $\hat p_n=x_1\cdots x_n(1+\theta_{n-1})$ 、 $|\theta_{n-1}|\le\gamma_{n-1}$ を得る。 $\varepsilon_j=-1$ に適用すると $x_1\cdots x_n=\hat p_n\prod_{j=2}^n(1+\delta_j)^{-1}=\hat p_n(1+\theta'_{n-1})$ 、 $|\theta'_{n-1}|\le\gamma_{n-1}$ を得る。 $F=F(2,3,-2,2)$ 、 $\operatorname{fl}$ を最近接偶数丸め、 $n=4$ 、 $x_1=\dots=x_4=5/4$ とすると、 $3u=3/8<1$ 、 $\gamma_3=3/5$ である。 $\hat p_2=\operatorname{fl}(25/16)=3/2$ 、 $\hat p_3=\operatorname{fl}(15/8)=2$ （最近接点 $7/4=7s_0$ と $2=4s_1$ のうち $m$ が偶数である方）、 $\hat p_4=\operatorname{fl}(5/2)=5/2$ であり、厳密な積 $25/16$ 、 $15/8$ 、 $5/2$ はいずれも正規範囲 $[1/4,7]$ に属する。厳密な値は $(5/4)^4=625/256$ であり、 $\theta_3=(5/2)/(625/256)-1=3/125$ 、 $\theta'_3=(625/256)/(5/2)-1=-3/128$ である。

注意 4.3. 「安定な和・内積・多項式評価」は、逐次和・内積・Horner 法の計算値を各演算の誤差因子の積で表し、補題 4.1によって後退誤差を $\gamma_k$ で評価する。

5 桁落ち・情報落ち・オーバーフロー・アンダーフロー

定義 5.1. $\beta\ge2$ と $k\ge1$ を整数とし、 $x,y$ を実数とする。 $xy>0$ 、 $x\ne y$ 、 $|x-y|\le\beta^{-k}\max\{|x|,|y|\}$ が成り立つとき、差 $x-y$ で基数 $\beta$ の $k$ 桁の桁落ち (cancellation) が起きるという。

命題 5.2. $x,y$ を $x\ne0$ 、 $y\ne0$ 、 $x\ne y$ を満たす実数とし、 $\rho\ge0$ とする。

$|\hat x-x|\le\rho|x|$ と $|\hat y-y|\le\rho|y|$ を満たす任意の実数 $\hat x,\hat y$ に対して $|(\hat x-\hat y)-(x-y)|\le\rho(|x|+|y|)$ が成り立つ。
$xy<0$ ならば $|x-y|=|x|+|y|$ であり、 $|\hat x-x|\le\rho|x|$ と $|\hat y-y|\le\rho|y|$ を満たす任意の実数 $\hat x,\hat y$ に対して $|(\hat x-\hat y)-(x-y)|/|x-y|\le\rho$ が成り立つ。
整数 $\beta\ge2$ 、 $k\ge1$ について差 $x-y$ で基数 $\beta$ の $k$ 桁の桁落ちが起き、 $\rho>0$ であるとする。このとき $\hat x:=x(1+\rho)$ 、 $\hat y:=y(1-\rho)$ は $|\hat x-x|=\rho|x|$ 、 $|\hat y-y|=\rho|y|$ と $|(\hat x-\hat y)-(x-y)|/|x-y|=\rho(|x|+|y|)/|x-y|\ge\rho\beta^k$ を満たす。

証明. $(\hat x-\hat y)-(x-y)=(\hat x-x)-(\hat y-y)$ に三角不等式を適用して(1)を得る。 $xy<0$ ならば $x$ と $-y$ は同符号であるから $|x-y|=|x|+|-y|=|x|+|y|$ である。(1)により $|(\hat x-\hat y)-(x-y)|\le\rho(|x|+|y|)=\rho|x-y|$ であり、(2)が成り立つ。差 $x-y$ で基数 $\beta$ の $k$ 桁の桁落ちが起きるならば $xy>0$ である。 $\hat x=x(1+\rho)$ 、 $\hat y=y(1-\rho)$ に対して $(\hat x-\hat y)-(x-y)=\rho(x+y)$ であり、 $x,y$ が同符号であるから $|x+y|=|x|+|y|$ である。桁落ちの定義により $|x|+|y|\ge\max\{|x|,|y|\}\ge\beta^k|x-y|$ であるから、(3)が成り立つ。▨

例 5.3. $F=F(2,53,-1022,1023)$ 、 $\operatorname{fl}$ を最近接偶数丸めとし、 $x=11/10$ 、 $y=1$ とする。 $\hat x:=\operatorname{fl}(11/10)=4953959590107546\cdot2^{-52}$ であり、 $\hat x-11/10=2^{-51}/5$ 、その相対誤差は $2^{-50}/11\approx8.07\times10^{-17}$ である。 $\hat y:=\operatorname{fl}(1)=1$ である。 $\hat x-\hat y=450359962737050\cdot2^{-52}=7205759403792800\cdot2^{-56}$ であり、 $2^{52}\le7205759403792800\le2^{53}-1$ であるから補題 1.2 (1)により $\hat x-\hat y\in F$ である。系 3.2 (1)により $\operatorname{fl}(\hat x-\hat y)=\hat x-\hat y$ である。 $(\hat x-\hat y)-1/10=2^{-51}/5$ であるから、差の相対誤差は $2^{-50}\approx8.88\times10^{-16}=8u$ であり、入力の相対誤差の $11$ 倍である。 $\max\{|x|,|y|\}/|x-y|=(11/10)/(1/10)=11$ であり、 $2^3\le11<2^4$ であるから、差 $x-y$ で基数 $2$ の $3$ 桁の桁落ちが起き、 $4$ 桁の桁落ちは起きない。この値は命題 5.2 (1)の上界 $\rho(|x|+|y|)/|x-y|=21\rho$ （ $\rho=2^{-50}/11$ ）以下である。CPython で 1.1 - 1 を計算すると 0.10000000000000009 が表示される。

注意 5.4. $a,b,c$ を $a\ne0$ 、 $b>0$ 、 $ac>0$ 、 $b^2-4ac>0$ を満たす実数とする。二次方程式 $ax^2+bx+c=0$ の根 $(-b+\sqrt{b^2-4ac})/(2a)$ の分子は正の二数 $\sqrt{b^2-4ac}$ と $b$ の差であり、命題 5.2 (1)の比は

\frac{b+\sqrt{b^2-4ac}}{b-\sqrt{b^2-4ac}}=\frac{\bigl(b+\sqrt{b^2-4ac}\bigr)^2}{4ac}

である。この根の計算式の変形と、近い二つの浮動小数点数の差が $F$ に属することの一般形は、「安定な和・内積・多項式評価」で扱う。

定義 5.5. $F$ を浮動小数点数系、 $\operatorname{fl}$ を $F$ の最近接丸めとし、 $x,y\in F$ が $y\ne0$ と $|x+y|\le N_{\max}$ を満たすとする。 $\operatorname{fl}(x+y)=x$ が成り立つとき、加算 $x+y$ で情報落ち (absorption) が起きるという。

命題 5.6. $F=F(\beta,p,e_{\min},e_{\max})$ を浮動小数点数系、 $\operatorname{fl}$ を $F$ の最近接丸めとする。 $x\in F$ が $|x|\ge\beta^{e_{\min}}$ を満たし、 $e:=e(|x|)$ を補題 1.2の整数とする。実数 $y$ が $|x+y|\le N_{\max}$ 、 $|y|<s_e/2$ を満たし、さらに $|x|>\beta^e$ であるか $xy\ge0$ であるならば、 $\operatorname{fl}(x+y)=x$ が成り立つ。

証明. $t:=x+y$ 、 $s:=s_e$ と置く。 $F=-F$ であり、 $x$ が $t$ のただ一つの最近接点であるという性質は $(x,y)$ を $(-x,-y)$ に替えても保たれるので、 $x>0$ としてよい。 $|t-x|=|y|<s/2$ である。 $x'\in F$ 、 $x'\ne x$ とする。 $x'>x$ ならば $x<N_{\max}$ であり、補題 1.2 (4)により $x'\ge x+s$ である。 $y\ge0$ ならば $|x'-t|\ge x+s-t=s-|y|>s/2>|y|$ であり、 $y<0$ ならば $|x'-t|>x-t=|y|$ である。 $x'<x$ かつ $y\ge0$ ならば $|x'-t|>t-x=|y|$ である。 $x'<x$ かつ $y<0$ ならば、仮定から $x>\beta^e$ であり、補題 1.2 (1)により $x-s\in F$ かつ $(x-s,x)$ に $F$ の元は無いので $x'\le x-s$ である。したがって $|x'-t|\ge t-(x-s)=s-|y|>|y|$ である。いずれの場合も $|x'-t|>|x-t|$ であるから、 $x$ は $t$ のただ一つの最近接点であり、 $\operatorname{fl}(x+y)=x$ である。▨

例 5.7. $F=F(2,53,-1022,1023)$ 、 $x=1$ とする。 $e(1)=0$ 、 $s_0=2^{-52}$ である。 $y=2^{-54}$ は $|y|<2^{-53}$ と $xy>0$ を満たすので、命題 5.6により任意の最近接丸めで $\operatorname{fl}(1+2^{-54})=1$ であり、加算 $1+y$ で情報落ちが起きる。このとき $\operatorname{fl}(\operatorname{fl}(1+y)-1)=0\ne y$ である。一方 $1=(1+2^{-54})(1+\delta)$ の $\delta=-2^{-54}/(1+2^{-54})$ は $|\delta|<u$ を満たし、系 3.2 (2)の評価は和 $1+y$ に対して成り立つ。 $y=2^{-53}=s_0/2$ では $1+y$ の最近接点は $1$ と $1+2^{-52}$ の二つである。 $m(1)=2^{52}$ 、 $m(1+2^{-52})=2^{52}+1$ であるから、最近接偶数丸めでは $\operatorname{fl}(1+2^{-53})=1$ となり、他方の最近接点を選ぶ最近接丸めでは情報落ちは起きない。 $y=2^{-53}+2^{-105}=(2^{52}+1)s_{-53}\in F$ は $|y|>s_0/2$ を満たし、 $1+y$ のただ一つの最近接点は $1+2^{-52}$ である。 $y=-3\cdot2^{-55}\in F$ は $|y|<s_0/2$ を満たすが、 $x=\beta^0$ かつ $xy<0$ である。 $1$ の前の $F$ の元は $1-2^{-53}=(2^{53}-1)s_{-1}$ であり、 $1-3\cdot2^{-55}$ との距離は $2^{-55}$ 、 $1$ との距離は $3\cdot2^{-55}$ であるから、 $\operatorname{fl}(1+y)=1-2^{-53}\ne1$ である。

例 5.8. $F=F(2,53,-1022,1023)$ 、 $\operatorname{fl}$ を最近接偶数丸めとし、 $c:=\operatorname{fl}(10^{308})$ とする。 $c\ge10^{308}(1-u)>N_{\max}/2$ であるから、 $c+c$ はオーバーフローの範囲にあり、系 3.2の仮定は $(c+c)/2$ の最初の演算で満たされない。 binary64 の演算はこの加算の結果を $+\infty$ とし、 $(c+c)/2$ の計算値も $+\infty$ である。厳密な値 $(c+c)/2=c$ は正規範囲にある。 $c=ms_{1023}$ と書くと $c/2=ms_{1022}\in F$ であるから、系 3.2 (1)により $\operatorname{fl}(c/2)=c/2$ であり、 $c/2+c/2=c\in F$ から $\operatorname{fl}(c/2+c/2)=c$ である。計算式 $c/2+c/2$ のすべての演算の厳密な結果は正規範囲にある。

例 5.9. $F=F(2,53,-1022,1023)$ 、 $\operatorname{fl}$ を最近接偶数丸めとし、 $d:=\operatorname{fl}(10^{-200})$ とする。 $d$ は正規範囲にあり、 $d^2<2^{-1075}=s_{-1022}/2$ である。 $F$ の正の最小元は $s_{-1022}$ であるから、 $d^2$ のただ一つの最近接点は $0$ であり、 $\operatorname{fl}(d\cdot d)=0$ である。絶対誤差 $d^2$ は系 3.2 (3)の上界 $u\beta^{e_{\min}}=2^{-1075}$ 以下であり、相対誤差は $1$ である。続く除算の計算値は $\operatorname{fl}(0/d)=0$ であり、厳密な値 $(d\cdot d)/d=d$ とは異なる。 $x=3\cdot2^{-1023}$ 、 $y=2^{-1022}$ は正規化数であり、 $x-y=2^{-1023}$ はアンダーフローの範囲にある。系 3.2 (4)により $\operatorname{fl}(x-y)=2^{-1023}$ である。

6 演習

問題 6.1. $F=F(\beta,p,e_{\min},e_{\max})$ を浮動小数点数系、 $u$ をその単位丸め誤差、 $\operatorname{fl}$ を $F$ の最近接丸めとする。正規範囲の任意の実数 $z$ に対して $|\operatorname{fl}(z)-z|\le u|z|/(1+u)$ が成り立つことを示せ。また、整数 $e$ が $e_{\min}\le e\le e_{\max}$ と $\beta^e(1+u)\le N_{\max}$ を満たすならば、 $z=\beta^e(1+u)$ で等号が成り立つことを示せ。

解答.

$t:=|z|$ 、 $e:=e(t)$ 、 $s:=s_e$ と置く。 $s/2=u\beta^e$ である。補題 1.2 (5)により $\beta^e\le a\le t<a+s$ を満たす $a\in F$ があり、 $a+s\in F$ であるか $t=a$ である。 $d:=|\operatorname{fl}(z)-z|=\min_{y\in F}|y-t|$ と置く。 $t=a$ ならば $d=0$ である。 $a<t\le a+s/2$ ならば、 $d\le t-a$ であり、 $1-a/t$ は $t$ について増加するので

\frac dt\le1-\frac at\le1-\frac{a}{a+s/2}=\frac{s/2}{a+s/2}\le\frac{s/2}{\beta^e+s/2}=\frac{u}{1+u}

である。 $t>a+s/2$ ならば、 $t\ne a$ であるから $a+s\in F$ であり、 $d\le a+s-t<s/2$ かつ $t>a+s/2\ge\beta^e+s/2$ であるから、 $d/t<(s/2)/(\beta^e+s/2)=u/(1+u)$ である。

$z=\beta^e(1+u)=\beta^e+s/2$ とする。 $\beta^e\in F$ との距離は $s/2$ である。 $\beta^e<N_{\max}$ であるから、補題 1.2 (4)により $\beta^e$ より大きい $F$ の元は $\beta^e+s$ 以上であり、 $z$ との距離は $s/2$ 以上である。 $\beta^e$ より小さい $F$ の元と $z$ との距離は $s/2$ より大きい。したがって $d=s/2$ であり、 $d/z=(s/2)/(\beta^e+s/2)=u/(1+u)$ である。▨

問題 6.2.補題 1.3の証明を完成させよ。

解答.

補題 1.3 (1)を示す。 $F=-F$ であるから $k\ge0$ としてよい。 $k=0$ ならば $ks_e=0\in F$ である。 $\beta^{p-1}\le k\le\beta^p-1$ ならば $ks_e$ は正規化数である。 $k=\beta^p$ ならば $ks_e=\beta^{e+1}$ であり、 $\beta^{e+1}\le N_{\max}<\beta^{e_{\max}+1}$ から $e+1\le e_{\max}$ であるので、 $\beta^{e+1}=\beta^{p-1}s_{e+1}$ は正規化数である。 $1\le k\le\beta^{p-1}-1$ とする。 $k\beta^j\ge\beta^{p-1}$ を満たす最小の整数 $j\ge0$ をとると $j\ge1$ であり、 $k\beta^{j-1}\le\beta^{p-1}-1$ から $k\beta^j\le\beta^p-\beta\le\beta^p-1$ である。 $e-j\ge e_{\min}$ ならば $ks_e=(k\beta^j)s_{e-j}$ は指数 $e-j$ の正規化数である。 $e-j<e_{\min}$ ならば $i:=e-e_{\min}$ は $0\le i<j$ を満たし、 $j$ の最小性から $k\beta^i\le\beta^{p-1}-1$ であるので、 $ks_e=(k\beta^i)s_{e_{\min}}$ は非正規数である。

補題 1.3 (2)を示す。補題 1.2 (3)により $x=\sigma ms_e$ と書く。 $\beta^{p-1}\le m\le\beta^p-1$ ならば $|x|\in[\beta^e,\beta^{e+1})$ であるから $e=e(|x|)$ である。 $e=e_{\min}$ かつ $1\le m\le\beta^{p-1}-1$ ならば $|x|<\beta^{e_{\min}}$ であるから $e(|x|)=e_{\min}=e$ である。いずれの場合も $x=(\sigma m)s_{e(|x|)}$ 、 $1\le m\le\beta^p-1$ である。 $e_{\min}\le e'\le e(|x|)$ ならば $s_{e(|x|)}=\beta^{e(|x|)-e'}s_{e'}$ であるから、 $x$ は $s_{e'}$ の整数倍である。

補題 1.3 (3)を示す。 $t'<\beta^{e_{\min}}$ ならば $e(t)=e(t')=e_{\min}$ である。 $t<\beta^{e_{\min}}\le t'$ ならば $e(t)=e_{\min}\le e(t')$ である。 $\beta^{e_{\min}}\le t$ ならば $\beta^{e(t)}\le t\le t'<\beta^{e(t')+1}$ であるから $e(t)<e(t')+1$ であり、 $e(t)\le e(t')$ である。▨

§E20.1浮動小数点演算