今回は整数を使わずに、いわゆる「数値計算」っぽいアルゴリズムで平方根を求めていきます。そういうアルゴリズムであっても correctly-rounded な値を得られるところが面白いですね。
あまくだり
早速ですが、下記のようなアルゴリズムを考えます。
$\gdef\iversonbracket#1{\mathop{[#1]}}$
- 入力: 浮動小数点数 $x\in[1\lldot 4)$
- 出力: $\roundcirc{\sqrt{x}}$
- $y \gets \iversonbracket{x\lt 2}\hat p(a^{[1]}, x) + \iversonbracket{x\ge 2}\hat p(a^{[2]}, x)$ で初期化する。
- $y \gets 0.5 \otimes (y \oplus (x \oslash y) )$ で更新する。
- $y \gets 0.5 \otimes (y \oplus (x \oslash y) )$ で更新する。
- $\roundcirc{y \times (y \ominus 2^{-52}) + (-x)} \ge 0$ であれば下記を行う。
- $y\ominus 2^{-52}$ を出力する。
- $\roundcirc{y \times (y \oplus 2^{-52}) + (-x)} \lt 0$ であれば下記を行う。
- $y\oplus 2^{-52}$ を出力する。
- $y$ を出力する。
ここで、$a^{[1]} = \angled{a^{[1]}_0, a^{[1]}_1, a^{[1]}_2, a^{[1]}_3}$, $a^{[2]} = \angled{a^{[2]}_0, a^{[2]}_1, a^{[2]}_2, a^{[2]}_3}$ とし、各要素は下記の通りです。 $$ \begin{aligned} & \angled{a^{[1]}_0, a^{[1]}_1, a^{[1]}_2, a^{[1]}_3} = \langle \\ & \qquad {\small \phantom- 0.371351660146978}{\footnotesize 290732988625677}{\scriptsize 535310387611389}{\tiny 16015625}, \\ & \qquad {\small \phantom- 0.784942635287931}{\footnotesize 067544775487476}{\scriptsize 726993918418884}{\tiny 27734375}, \\ & \qquad {\small -0.180689144911217}{\footnotesize 069999764817112}{\scriptsize 009041011333465}{\tiny 576171875}, \\ & \qquad {\small \phantom- 0.024476908831259}{\footnotesize 320403761492457}{\scriptsize 306303549557924}{\tiny 2706298828125} \\ & \rangle, \\ & \angled{a^{[2]}_0, a^{[2]}_1, a^{[2]}_2, a^{[2]}_3} = \langle \\ & \qquad {\small \phantom- 0.525170554189621}{\footnotesize 108243102298729}{\scriptsize 354515671730041}{\tiny 50390625}, \\ & \qquad {\small \phantom- 0.555038260254535}{\footnotesize 065207903699047}{\scriptsize 164991497993469}{\tiny 23828125}, \\ & \qquad {\small -0.063883259826760}{\footnotesize 172006565596802}{\scriptsize 829531952738761}{\tiny 90185546875}, \\ & \qquad {\small \phantom- 0.004326947054267}{\footnotesize 089344536945105}{\scriptsize 801351019181311}{\tiny 130523681640625} \\ & \rangle. \end{aligned} $$ また、長さ $2$ 以上の浮動小数点数の列 $a$ に対して $\hat p(a, x)$ を下記で定義します。 $$ \begin{aligned} \hat p(\angled{a_0, a_1}, x) &= \roundcirc{a_1\times x+a_0}, \\ \hat p(\angled{a_0, a_1, \dots, a_n}, x) &= \roundcirc{\hat p(\angled{a_1, \dots, a_n}, x)\times x + a_0}. \end{aligned} $$
加えて、便宜上、$a(x) = \sum_{i=0}^n a_i x^i$ で定義しておきます。 $a(x)$ を浮動小数点型で計算しようとしたのが $\hat p(a, x)$ ということですね。
例
$x_1 = 1.5625 = 1.25^2$ とします。 $$ \begin{aligned} y_1^{(0)} &= \hat p(a^{[1]}, x) \\ &= {\small 1.250060917280526}{\footnotesize 817663144356629}{\scriptsize 345566034317016}{\tiny 6015625}, \\ y_1^{(1)} &= 0.5\otimes(y_1^{(0)}\oplus(x\oslash y_1^{(0)}) ) \\ &= {\small 1.250000001484293}{\footnotesize 576936579484026}{\scriptsize 879072189331054}{\tiny 6875}, \\ y_1^{(2)} &= 0.5\otimes(y_1^{(1)}\oplus(x\oslash y_1^{(1)}) ) \\ &= {\small 1.25}. \end{aligned} $$ $x_2 = 2.25 = 1.5^2$ とします。 $$ \begin{aligned} y_2^{(0)} &= \hat p(a^{[2]}, x) \\ &= {\small 1.499884268179362}{\footnotesize 711848057188035}{\scriptsize 454601049423217}{\tiny 7734375}, \\ y_2^{(1)} &= 0.5\otimes(y_2^{(0)}\oplus(x\oslash y_2^{(0)}) ) \\ &= {\small 1.500000004464962}{\footnotesize 621852919255616}{\scriptsize 143345832824707}{\tiny 03125}, \\ y_2^{(2)} &= 0.5\otimes(y_2^{(1)}\oplus(x\oslash y_2^{(1)}) ) \\ &= {\small 1.5}. \end{aligned} $$
$\hat p(a^{[\ast]}, x)$ の計算で $\sqrt x$ にある程度近い値が得られた後、素早く $\sqrt x$ に収束していることが見て取れます。
グラフ
$\hat p(a^{[1]}, x)$ は $[1\lldot 2)$ の区間で、$\hat p(a^{[2]}, x)$ は $[2\lldot 4)$ の区間で、$\sqrt x$ にぺたーっとくっついていることがわかります。
前提知識たち
上記のアルゴリズムを理解するにあたって必要な知識の紹介や、補題の証明をします。
Equioscillation theorem
等振動定理 (equioscillation theorem) は、実数値関数を多項式で近似するにあたって重要な定理です。
区間 $[a\lldot b]$ で定義される連続な実数値関数全体からなる集合を $\mathscr C[a, b]$、高々 $n$ 次の実数係数多項式全体からなる集合を $\mathscr P_n$ と書きます。また、$f\in\mathscr C[a, b]$ に対し、$\|f\|_{\infty}^{[a\lldot b]} = \sup_{x\in[a\lldot b]} |f(x)|$ で定義します。
Theorem 1 (Chebyshev): $f\in\mathscr C[a, b]$ と $p^\ast\in\mathscr P_n$ に対し、$p^\ast = \argmin_{p\in\mathscr P_n} \|f-p\|_{\infty}^{[a\lldot b]}$ となる必要十分条件は、 $n+2$ 個の点 $a\le x_0\lt x_1\lt \cdots \lt x_{n+1}\le b$ と $\sigma\in\{-1, 1\}$ が存在して、 $$ f(x_i)-p^\ast(x_i) = \sigma\cdot(-1)^i\|f-p^\ast\|_{\infty}^{[a\lldot b]} $$ を満たすことである。
この $p^\ast$ を minimax polynomial と呼びます。これを求めるアルゴリズムを、次で紹介します。
Remez algorithm
Євген Якович Ремез (Evgeny Yakovlevich Remez) によるアルゴリズムです。 実数 $a$, $b$ と関数 $f\in\mathscr C[a, b]$ と非負整数 $n$ に対し、$\argmin_{p\in\mathscr P_n} \|f-p\|_{\infty}^{[a\lldot b]}$（に近いもの）を返します。
Sollya では remez(f, n, [a; b]) のようにして Remez algorithm を使うことができます。
> prec = 20; The precision has been set to 20 bits. > display = hexadecimal; Display mode is hexadecimal numbers. > remez(x^0.5, 3, [1; 2]); 0x1.7c462p-2 + x * (0x1.91e14p-1 + x * (-0x1.7205cp-3 + x * 0x1.90fb2p-6))
今回の記事ではアルゴリズムの実装や定理の証明には触れず、ブラックボックスとして使います。出力される多項式を評価すれば事足りるためです。
精度などの都合で、$p^\ast$ に完全に一致するものが得られるわけではないため、出力される多項式に関する誤差評価自体はどうしても必要になってきます。
Claim 2: 浮動小数点数 $x\in[1\lldot 2)$ に対して、$|a^{[1]}(x)-\sqrt x|\lt 8.20594\times 10^{-5}$ が成り立つ。
Proof
$f(x) = a^{[1]}(x)-\sqrt x$ とし、$f(x)$ の極値を考える。 $$ \begin{aligned} & \angled{m_0, m_1, m_2, m_3} = \langle \\ & \qquad {\phantom- 6689676793045386}, \\ & \qquad {\phantom- 7070134719579883}, \\ & \qquad {-6510012525536406}, \\ & \qquad {\phantom- 7054988639465029} \\ & \rangle, \\ & \angled{e_0, e_1, e_2, e_3} = \angled{-54, -53, -55, -58} \end{aligned} $$ とする。ここで、各 $i\in\{0, 1, 2, 3\}$ に対して $a^{[1]}_i = m_i\times 2^{e_i}$ かつ $2^{52}\le m_i\lt 2^{-53}$ が成り立つ。
このとき、整数 $2^{52}\le m_x\lt 2^{53}$ に対して $f'(2^{-52}\cdot m_x) \gt 0$ は、下記の 2 式が成り立つことと同値である。 $$ \begin{cases} m_1\cdot 2^{2\cdot 52+e_1 + 58} + 2 m_2\cdot 2^{52+e_2 + 58}\cdot m_x + 3 m_3\cdot 2^{e_3 + 58}\cdot m_x^2 \ge 0, \\ 4m_x\cdot ( m_1\cdot 2^{2\cdot 52 + e_1 + 58} + 2m_2\cdot 2^{52+e_2+58}\cdot m_x + 3m_3\cdot 2^{e_3+58}\cdot m_x^2 )^2 - 2^{5\cdot 52+2\cdot 58} \gt 0. \end{cases} $$ 同様に、$f''(2^{-52}\cdot m_x) \gt 0$ は、下記の 2 式が成り立つことと同値である。 $$ \begin{cases} 2m_2\cdot 2^{e_2+58} + 6m_3\cdot 2^{e_3+58}\cdot 2^{52}\cdot m_x \ge 0, \\ 16m_x^3\cdot ( 2m_2 \cdot 2^{52+e_2+58} + 6m_3\cdot 2^{e_3+58}\cdot m_x )^2 - 2^{5\cdot 52+2\cdot 58} \lt 0. \end{cases} $$ さらに $m_3\ge 0$ を踏まえて、$f'''(2^{-52}\cdot m_x)\gt 0$ は下記が成り立つことと同値である。 $$ m_x^5\cdot (16m_3)^2 - 2^{5\cdot 52-2e_3} \gt 0. $$ よってこれらは（多倍長の）整数型を用いて二分探索ができ、その結果を用いて次の増減表を書くことができる。
| $1$ | $\cdots$ | $x_3^-$ | $x_3^+$ | $\cdots$ | $2$ | |
|---|---|---|---|---|---|---|
| $f''(x)$ | $+$ | $\searrow$ | $-$ | $-$ | $\nearrow$ | $+$ |
| $f'''(x)$ | $-$ | $\cdots$ | $-$ | $+$ | $\cdots$ | $+$ |
$$ \begin{aligned} x_3^- &= 6552532042542083 \times 2^{-52} \\ &= {\small 1.454954388644865}{\footnotesize 259649918698414}{\scriptsize 694517850875854}{\tiny 4921875}, \\ x_3^+ &= 6552532042542084 \times 2^{-52} \\ &= {\small 1.454954388644865}{\footnotesize 481694523623446}{\scriptsize 002602577209472}{\tiny 65625}. \end{aligned} $$
| $1$ | $\cdots$ | $x_{2, 1}^-$ | $x_{2, 1}^+$ | $\cdots$ | $x_{2, 2}^-$ | $x_{2, 2}^+$ | $\cdots$ | $2$ | |
|---|---|---|---|---|---|---|---|---|---|
| $f'(x)$ | $-$ | $\nearrow$ | $+$ | $+$ | $\searrow$ | $-$ | $-$ | $\nearrow$ | $+$ |
| $f''(x)$ | $+$ | $\cdots$ | $+$ | $-$ | $\cdots$ | $-$ | $+$ | $\cdots$ | $+$ |
$$ \begin{aligned} x_{2, 1}^- &= 5705497900301872 \times 2^{-52} \\ &= {\small 1.266875027173124}{\footnotesize 834689588169567}{\scriptsize 286968231201171}{\tiny 875}, \\ x_{2, 1}^+ &= 5705497900301873 \times 2^{-52} \\ &= {\small 1.266875027173125}{\footnotesize 056734193094598}{\scriptsize 595052957534790}{\tiny 0390625}, \\ x_{2, 2}^- &= 7549905344458297 \times 2^{-52} \\ &= {\small 1.676415749431624}{\footnotesize 968579512824362}{\scriptsize 609535455703735}{\tiny 3515625}, \\ x_{2, 2}^+ &= 7549905344458298 \times 2^{-52} \\ &= {\small 1.676415749431625}{\footnotesize 190624117749393}{\scriptsize 917620182037353}{\tiny 515625}. \end{aligned} $$
| $1$ | $\cdots$ | $x_{1, 1}^-$ | $x_{1, 1}^+$ | $\cdots$ | $x_{1, 2}^-$ | $x_{1, 2}^+$ | $\cdots$ | $x_{1, 3}^-$ | $x_{1, 3}^+$ | $\cdots$ | $2$ | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| $f(x)$ | $\searrow$ | $\cdots$ | $\searrow$ | $\nearrow$ | $\cdots$ | $\nearrow$ | $\searrow$ | $\cdots$ | $\searrow$ | $\nearrow$ | $\cdots$ | $\nearrow$ |
| $f'(x)$ | $-$ | $\cdots$ | $-$ | $+$ | $\cdots$ | $+$ | $-$ | $\cdots$ | $-$ | $+$ | $\cdots$ | $+$ |
$$ \begin{aligned} x_{1, 1}^- &= 5093190532915904 \times 2^{-52} \\ &= {\small 1.130915479689221}{\footnotesize 569969959091395}{\scriptsize 139694213867187}{\tiny 5}, \\ x_{1, 1}^+ &= 5093190532915905 \times 2^{-52} \\ &= {\small 1.130915479689221}{\footnotesize 792014564016426}{\scriptsize 447778940200805}{\tiny 6640625}, \\ x_{1, 2}^- &= 6620254117772346 \times 2^{-52} \\ &= {\small 1.469991710084072}{\footnotesize 256137233125627}{\scriptsize 972185611724853}{\tiny 515625}, \\ x_{1, 2}^+ &= 6620254117772347 \times 2^{-52} \\ &= {\small 1.469991710084072}{\footnotesize 478181838050659}{\scriptsize 280270338058471}{\tiny 6796875}, \\ x_{1, 3}^- &= 8282230118499356 \times 2^{-52} \\ &= {\small 1.839024514560384}{\footnotesize 737649201269960}{\scriptsize 030913352966308}{\tiny 59375}, \\ x_{1, 3}^+ &= 8282230118499357 \times 2^{-52} \\ &= {\small 1.839024514560384}{\footnotesize 959693806194991}{\scriptsize 338998079299926}{\tiny 7578125}. \end{aligned} $$
よって、$S = \{1, x_{1, 1}^-, x_{1, 1}^+, x_{1, 2}^-, x_{1, 2}^+, x_{1, 3}^-, x_{1, 3}^+, 2\}$ に対して $\max_{x\in S} |a^{[1]}(x)-{\sqrt x}|$ を求めればよい。 数値計算により、$\max_{x\in S} |a^{[1]}(x)-{\sqrt x}| \lt 8.20594\times 10^{-5}$ とわかる。$\qed$
Claim 3: 浮動小数点数 $x\in[2\lldot 4)$ に対して、$|a^{[2]}(x)-\sqrt x|\lt 1.16050\times 10^{-4}$ が成り立つ。
Proof
$f(x) = a^{[2]}(x)-\sqrt x$ とし、$f(x)$ の極値を考える。 $$ \begin{aligned} & \angled{m_0, m_1, m_2, m_3} = \langle \\ & \qquad {\phantom- 4730315824308669}, \\ & \qquad {\phantom- 4999340204117385}, \\ & \qquad {-4603274002416155}, \\ & \qquad {\phantom- 4988630308159777} \\ & \rangle, \\ & \angled{e_0, e_1, e_2, e_3} = \angled{-53, -53, -56, -60} \end{aligned} $$ とする。ここで、各 $i\in\{0, 1, 2, 3\}$ に対して $a^{[2]}_i = m_i\times 2^{e_i}$ かつ $2^{52}\le m_i\lt 2^{-53}$ が成り立つ。
よって、Claim 2 同様にして増減表を書くことができる。
| $2$ | $\cdots$ | $x_3^-$ | $x_3^+$ | $\cdots$ | $4$ | |
|---|---|---|---|---|---|---|
| $f''(x)$ | $+$ | $\searrow$ | $-$ | $-$ | $\nearrow$ | $+$ |
| $f'''(x)$ | $-$ | $\cdots$ | $-$ | $+$ | $\cdots$ | $+$ |
$$ \begin{aligned} x_3^- &= 6552532042542083 \times 2^{-51} \\ &= {\small 2.909908777289730}{\footnotesize 519299837396829}{\scriptsize 389035701751708}{\tiny 984375}, \\ x_3^+ &= 6552532042542084 \times 2^{-51} \\ &= {\small 2.909908777289730}{\footnotesize 963389047246892}{\scriptsize 005205154418945}{\tiny 3125}. \end{aligned} $$
| $2$ | $\cdots$ | $x_{2, 1}^-$ | $x_{2, 1}^+$ | $\cdots$ | $x_{2, 2}^-$ | $x_{2, 2}^+$ | $\cdots$ | $4$ | |
|---|---|---|---|---|---|---|---|---|---|
| $f'(x)$ | $-$ | $\nearrow$ | $+$ | $+$ | $\searrow$ | $-$ | $-$ | $\nearrow$ | $+$ |
| $f''(x)$ | $+$ | $\cdots$ | $+$ | $-$ | $\cdots$ | $-$ | $+$ | $\cdots$ | $+$ |
$$ \begin{aligned} x_{2, 1}^- &= 5705497900301873 \times 2^{-51} \\ &= {\small 2.533750054346250}{\footnotesize 113468386189197}{\scriptsize 190105915069580}{\tiny 078125}, \\ x_{2, 1}^+ &= 5705497900301874 \times 2^{-51} \\ &= {\small 2.533750054346250}{\footnotesize 557557596039259}{\scriptsize 806275367736816}{\tiny 40625}, \\ x_{2, 2}^- &= 7549905344458296 \times 2^{-51} \\ &= {\small 3.352831498863249}{\footnotesize 493069815798662}{\scriptsize 602901458740234}{\tiny 375}, \\ x_{2, 2}^+ &= 7549905344458297 \times 2^{-51} \\ &= {\small 3.352831498863249}{\footnotesize 937159025648725}{\scriptsize 219070911407470}{\tiny 703125}. \end{aligned} $$
| $2$ | $\cdots$ | $x_{1, 1}^-$ | $x_{1, 1}^+$ | $\cdots$ | $x_{1, 2}^-$ | $x_{1, 2}^+$ | $\cdots$ | $x_{1, 3}^-$ | $x_{1, 3}^+$ | $\cdots$ | $4$ | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| $f(x)$ | $\searrow$ | $\cdots$ | $\searrow$ | $\nearrow$ | $\cdots$ | $\nearrow$ | $\searrow$ | $\cdots$ | $\searrow$ | $\nearrow$ | $\cdots$ | $\nearrow$ |
| $f'(x)$ | $-$ | $\cdots$ | $-$ | $+$ | $\cdots$ | $+$ | $-$ | $\cdots$ | $-$ | $+$ | $\cdots$ | $+$ |
$$ \begin{aligned} x_{1, 1}^- &= 5093190532915894 \times 2^{-51} \\ &= {\small 2.261830959378438}{\footnotesize 699047819682164}{\scriptsize 117693901062011}{\tiny 71875}, \\ x_{1, 1}^+ &= 5093190532915895 \times 2^{-51} \\ &= {\small 2.261830959378439}{\footnotesize 143137029532226}{\scriptsize 733863353729248}{\tiny 046875}, \\ x_{1, 2}^- &= 6620254117772374 \times 2^{-51} \\ &= {\small 2.939983420168156}{\footnotesize 946772342053009}{\scriptsize 197115898132324}{\tiny 21875}, \\ x_{1, 2}^+ &= 6620254117772375 \times 2^{-51} \\ &= {\small 2.939983420168157}{\footnotesize 390861551903071}{\scriptsize 813285350799560}{\tiny 546875}, \\ x_{1, 3}^- &= 8282230118499338 \times 2^{-51} \\ &= {\small 3.678049029120761}{\footnotesize 481692625238792}{\scriptsize 970776557922363}{\tiny 28125}, \\ x_{1, 3}^+ &= 8282230118499339 \times 2^{-51} \\ &= {\small 3.678049029120761}{\footnotesize 925781835088855}{\scriptsize 586946010589599}{\tiny 609375}. \end{aligned} $$
よって、$S = \{2, x_{1, 1}^-, x_{1, 1}^+, x_{1, 2}^-, x_{1, 2}^+, x_{1, 3}^-, x_{1, 3}^+, 4\}$ に対して $\max_{x\in S} |a^{[2]}(x)-{\sqrt x}|$ を求めればよい。 数値計算により、$\max_{x\in S} |a^{[2]}(x)-{\sqrt x}| \lt 1.16050\times 10^{-4}$ とわかる。$\qed$
note: Sollya で infnorm なり findzeros なりを使う方法もありますが、証明のコアの部分に魔法の技術を使うのはだめだと思ったので避けました。
Horner’s method
多項式 $\sum_{i=0}^n a_ix^n$ を計算する際に $( (\cdots(a_n\times x+a_{n-1})\times x+\cdots )\times x + a_0)$ と計算する方法を Horner’s method (Horner’s scheme) と呼びます。名前の由来は William George Horner という人ですが、この方法自体はもっと古くからあるようです。$\hat p(a, x)$ は、FMA を用いつつ Horner’s method で計算するアルゴリズムに相当します。
$\hat p(a, x)$ と $a(x)$ の誤差を評価しましょう。
任意の実数 $x$ に対して、$|{\roundcirc{x}}|\lt\infty$ であれば、ある実数 $|\theta|\lt 2^{-53}$ が存在して $\roundcirc{x} = (1+\theta)x$ となります。
Claim 4: ある実数 $|\theta_1|, \dots, |\theta_n|\le 2^{-53}$ が存在して下記が成り立つ。 $$ \hat p(\angled{a_0, \dots, a_n}, x) = \left(\sum_{i=0}^{n-1} \left(\prod_{j=1}^{i+1} (1+\theta_j)\right)a_ix^i\right) + \left(\prod_{j=1}^{n}(1+\theta_j)\right)a_nx^n. $$
Proof
実数 $|\theta_1|, \dots, |\theta_n|\le 2^{-53}$ を用いて $$ \begin{aligned} \hat p(\angled{a_{n-1}, a_n}, x) &= (1+\theta_n)(a_nx+a_{n-1}), \\ \hat p(\angled{a_i, a_{i+1}, \dots, a_n}, x) &= (1+\theta_{i+1})( (\hat p(\angled{a_{i+1}, \dots, a_n}, x)x+a_i) \end{aligned} $$ と表せる。
$0\le k\le n-1$ に対し、ある実数 $|\theta_{k+1}|, \dots, |\theta_n|\le 2^{-53}$ が存在して $$ \hat p(\angled{a_k, \dots, a_n}, x) = \left(\sum_{i=k}^{n-1} \left(\prod_{j=k+1}^{i+1} (1+\theta_j)\right)a_ix^{i-k}\right) + \left(\prod_{j=k+1}^{n}(1+\theta_j)\right)a_nx^{n-k} $$ と表せることを $P(k)$ と書く。$P(0)$ を示す。
To-be-proved 1: $P(n-1)$
$$ \begin{aligned} &\phantom{{}={}} \hat p(\angled{a_{n-1}, a_n}, x) \\ &= \left(\sum_{i=n-1}^{n-1} \left(\prod_{j=n}^{i+1} (1+\theta_j)\right)a_ix^{i-(n-1)}\right) + \left(\prod_{j=1}^{n}(1+\theta_j)\right)a_nx^{n-(n-1)} \\ &= (1+\theta_n)a_{n-1}x^0 + (1+\theta_n)a_nx^1 \\ &= (1+\theta_n)(a_nx+a_{n-1}). \end{aligned} $$
To-be-proved 2: $k\ge 1\wedge P(k) \implies P(k-1)$
$$ \begin{aligned} &\phantom{{}={}} \hat p(\angled{a_{k-1}, a_k, \dots, a_n}, x) \\ &= (\hat p(\angled{a_k, \dots, a_n}, x)x+a_{k-1})(1+\theta_k) \\ &= \hat p(\angled{a_k, \dots, a_n}, x)\cdot (1+\theta_k) x+a_{k-1}(1+\theta_k) \\ &= \left(\sum_{i=k}^{n-1} \left(\prod_{j=k}^{i+1} (1+\theta_j)\right)a_ix^{i-k+1}\right) + \left(\prod_{j=k}^{n}(1+\theta_j)\right)a_nx^{n-k+1} + (1+\theta_k)a_{k-1} \\ &= \left(\sum_{i=k-1}^{n-1} \left(\prod_{j=(k-1)+1}^{i+1} (1+\theta_j)\right)a_ix^{i-(k-1)}\right) \\ &\qquad + \left(\prod_{j=(k-1)+1}^{n}(1+\theta_j)\right)a_nx^{n-(k-1)}.\quad\qed \end{aligned} $$
Corollary 5: $a = \angled{a_0, \dots, a_n}$ に対し、下記が成り立つ。 $$ |\hat p(a, x) - a(x)| \le ( (1+2^{-53})^n-1)\sum_{i=0}^n |a_ix^i|. $$
Proof
各 $1\le i\le n$ に対して $1+\delta_i=\prod_{j=1}^i (1+\theta_j)$ が成り立つように $\delta_i$ を定めると、 $$ \begin{aligned} &\phantom{{}={}} {\left|\hat p(\angled{a_0, \dots, a_n}, x) - \left(\sum_{i=0}^n a_ix^i\right)\right|} \\ &= \left|\left(\sum_{i=0}^{n-1} \delta_i a_ix^i\right) + \delta_n a_nx^n\right| \\ &\le ( (1+2^{-53})^n-1)\sum_{i=0}^n |a_ix^i|. \qquad\qed \end{aligned} $$
Newton’s method
Taylor theorem などを用いた一般論はまたの機会にします。（微分などによって導出される）更新式を反復することで、真の値との誤差を素早く小さくする方法です。
Lemma 6: 実数 $x$, $y$ ($x\ge 0$) に対して、$\varepsilon = y-\sqrt x$ と定義する。このとき、 $$ y' = \frac12\cdot\left(y+\frac xy\right) $$ で定義し、$\varepsilon' = y'-\sqrt x$ とすると、 $$ \varepsilon' = \frac{\varepsilon^2}{2y} $$ が成り立つ。
Proof
$$ \begin{aligned} y' &= \frac12\cdot\left( (\sqrt x+\varepsilon) + \frac x{\sqrt x+\varepsilon}\right) \\ &= \frac12\cdot\left( (\sqrt x+\varepsilon) + \frac{x-\varepsilon^2+\varepsilon^2}{\sqrt x+\varepsilon}\right) \\ &= \frac12\cdot\left( (\sqrt x+\varepsilon) + (\sqrt x-\varepsilon) + \frac{\varepsilon^2}{\sqrt x+\varepsilon}\right) \\ &= \sqrt x + \frac12\cdot\frac{\varepsilon^2}{\sqrt x+\varepsilon} \\ &= \sqrt x + \frac{\varepsilon^2}{2y}. \quad\qed \end{aligned} $$
Lemma 7: 浮動小数点数 $x, y^{(i)} \in[1\lldot 4)$ に対して、$\varepsilon_i = y^{(i)} - \sqrt x$ と定義する。また、$y^{(i+1)} = 0.5\otimes(y^{(i)}\oplus(x\oslash y^{(i)}) )$ および $\varepsilon_{i+1} = y^{(i+1)}-\sqrt x$ とする。$|\varepsilon_i|\le\sqrt{4-3\cdot 2^{-52}}$ ならば、
- $y^{(i+1)}\in[1\lldot 4)$, and
- $\varepsilon_{i+1} \in [-3\cdot 2^{-53}\lldot \frac12\varepsilon_i^2+3\cdot 2^{-53}]$
が成り立つ。なお、$\sqrt{4-3\cdot 2^{-52}}\gt 1.999$ が成り立つ。
Proof
$x$, $y^{(i)}$ の範囲より、$\frac xy\in[\frac14\lldot 4)$ である。ある実数 $|\theta_{\oslash}|\le 2^{-52}$ が存在し、$x\oslash y^{(i)} = \frac x{y^{(i)}}+\theta_{\oslash}$ が成り立つ。
To-be-proved 1: $y^{(i+1)}\ge 1$
$\frac12(y^{(i)}+(x\oslash y^{(i)}) )\ge 1-2^{-54}$ ならば $y^{(i+1)} = 0.5\otimes(y^{(i)}\oplus(x\oslash y^{(i)}) )\ge 1$ なので、$y^{(i)}+(x\oslash y^{(i)})\ge 2-2^{-53}$ を示せばよい。
Case 1: $\frac x{y^{(i)}}\in [2\lldot 4)$
$y^{(i)}\ge 1$ より、 $$ \begin{aligned} y^{(i)}+\frac x{y^{(i)}}+\theta_{\oslash} &\ge 1+2-2^{-52} \\ &\ge 2-2^{-53}. \end{aligned} $$
Case 2: $\frac x{y^{(i)}}\in [\frac14\lldot 2)$
このとき、$|\theta_{\oslash}|\le 2^{-53}$ が成り立つ。$x, y^{(i)}\ge 1$ より、相加・相乗平均の関係から $$ \begin{aligned} y^{(i)}+\frac x{y^{(i)}}+\theta_{\oslash} &\ge 2\cdot\sqrt{y^{(i)}\cdot\frac{x\mathstrut}{y^{(i)}}} + \theta_{\oslash} \\ &= 2\sqrt x + \theta_{\oslash} \\ &\ge 2-2^{-53}. \end{aligned} $$
To-be-proved 2: $y^{(i+1)}\lt 4$
$\frac12(y^{(i)}+(x\oslash y^{(i)}) )\le 4-2^{-52}$ ならば $y^{(i+1)} \lt 4$ なので、$y^{(i)}+(x\oslash y^{(i)})\le 8-2^{-51}$ を示せばよい。 $|\theta_{\oslash}|\le 2^{-52}$ より、$y^{(i)}+\frac x{y^{(i)}}\le 8-3\cdot 2^{-52}$ を示せば十分。
Lemma 6 より $$ \begin{aligned} y^{(i)}+\frac x{y^{(i)}} &= 2\sqrt x + \frac{\varepsilon_i^2}{y^{(i)}} \\ &\le 2\cdot 2 + \frac{\sqrt{4-3\cdot 2^{-52}}^2}1 \\ &= 4 + 4-3\cdot 2^{-52} \end{aligned} $$
To-be-proved 3: $\varepsilon_{i+1}\in [-3\cdot 2^{-53}\lldot \frac12\varepsilon_i^2+3\cdot 2^{-53}]$
$y^{(i)}+(x\oslash y^{(i)})\in[2\lldot 8)$ より、ある実数 $|\theta_{\oplus}|\le 2^{-51}$ が存在して、下記が成り立つ。 $$ \begin{aligned} y^{(i+1)} &= \frac12\cdot\left(y^{(i)}+\frac x{y^{(i)}}+\theta_{\oslash}+\theta_{\oplus}\right) \\ &= \frac12\cdot\left(2\sqrt x+\frac{\varepsilon_i^2}{y^{(i)}} + \theta_{\oslash}+\theta_{\oplus}\right) \\ &= \sqrt x + \frac{\varepsilon_i^2}{2y^{(i)}} + \frac12\cdot\left(\theta_{\oslash}+\theta_{\oplus}\right). \end{aligned} $$ よって、 $$ \begin{aligned} \varepsilon_{i+1} &= y^{(i+1)}-\sqrt x \\ &= \frac{\varepsilon_i^2}{2y^{(i)}} + \frac12\cdot\left(\theta_{\oslash}+\theta_{\oplus}\right) \\ &\le \frac{\varepsilon_i^2}2 + \frac12\cdot(2^{-52} + 2^{-51}) \\ &\le \frac{\varepsilon_i^2}2 + 3\cdot 2^{-53}. \end{aligned} $$ また、 $$ \begin{aligned} \varepsilon_{i+1} &= \frac{\varepsilon_i^2}{2y^{(i)}} + \frac12\cdot\left(\theta_{\oslash}+\theta_{\oplus}\right) \\ &\ge 0 - \frac12\cdot(2^{-52} + 2^{-51}) \\ &\ge - 3\cdot 2^{-53}. \quad\qed \end{aligned} $$
$\delta_- = 3\cdot 2^{-53}$ および $\delta_+ = 2^{-26.5}$ とします。 $\delta_- \gt 3.330\times 10^{-16}$ および $\delta_+ \gt 1.053\times 10^{-8}$ が成り立ちます。
Lemma 8: Lemma 7 の条件に加えて $\varepsilon_i\in[-\delta_-\lldot 0]$ が成り立つとき、$y^{(i+1)}\le 2$ が成り立つ。
Proof
$\frac x{y^{(i)}}\le 4$ であり、ある実数 $|\theta_{\oslash}|\le 2^{-52}$ が存在して $x\oslash y^{(i)} = (x\div y^{(i)})+\theta_{\oslash}$ が成り立つ。
$$ \begin{aligned} &\phantom{{}={}} \tfrac12\cdot (y^{(i)} + (x\oslash y^{(i)}) ) \\ &= \sqrt x + \frac{\varepsilon_i^2}{2y^{(i)}} + \frac{\theta_{\oslash}}2 \\ &\le 2 + 9\cdot 2^{-2\cdot 53-1} + 2^{-53} \\ &\lt 2 + 2^{-52}. \end{aligned} $$ よって、$y^{(i+1)}\le 2$ が成り立つ。$\qed$
Lemma 9: Lemma 7 の条件に加えて $\varepsilon_i\in[0\lldot \delta_+]$ が成り立つとき、$y^{(i+1)}\le 2$ が成り立つ。
Proof
$$ \begin{aligned} \frac x{y^{(i)}} &= \frac x{\sqrt x+\varepsilon_i} \\ &\le \frac x{\sqrt x+0} \\ &= \sqrt x\le 2 \end{aligned} $$ より、ある実数 $|\theta_{\oslash}|\le 2^{-53}$ が存在して $x\oslash y^{(i)} = (x\div y^{(i)})+\theta_{\oslash}$ が成り立つ。
$$ \begin{aligned} &\phantom{{}={}} \tfrac12\cdot (y^{(i)} + (x\oslash y^{(i)}) ) \\ % &= (\sqrt x+\varepsilon_i) + \left(\sqrt x-\varepsilon_i + \frac{\varepsilon_i^2}{\sqrt x+\varepsilon_i} + \theta_{\oslash}\right) \\ &= \sqrt x + \left(\frac{\varepsilon_i^2}{2(\sqrt x+\varepsilon_i)} + \frac{\theta_{\oslash}}2\right) \\ &\le 2 + \frac{\delta_+^2}{2\cdot(1+0)} + 2^{-53-1} \\ &\lt 2 + 2^{-52}. \end{aligned} $$ よって、$y^{(i+1)}\le 2$ が成り立つ。$\qed$
Lemma 10: $x\in(0\lldot \tfrac14)$ ならば $\sqrt{1-3x} \gt 1-2x$ が成り立つ。
Proof
$$ \begin{aligned} \sqrt{1-3x}\gt 1-2x &\iff 1-3x \gt 1-4x+x^2 \\ &\iff x\gt x^2 \\ &\iff (x-\tfrac12)^2\lt \tfrac14. \quad\qed \end{aligned} $$
Lemma 11: Lemma 7 の条件に加えて $\varepsilon_i\in[-\delta_-\lldot 0]$ かつ $x/y^{(i)}\in(2\lldot 4)$ ならば、$(x, y) = (4-2^{-51}, 2-2^{-51})$ が成り立つ。
Proof
$y^{(i)}\in[\sqrt x-3\cdot 2^{-53}\lldot\sqrt x]\cap (x/4\lldot x/2)$。
$y^{(i)}\in[\sqrt x-3\cdot 2^{-53}\lldot x/2)\ne\emptyset$ が必要。$x\ge 0$ に気をつけて $\sqrt x-3\cdot 2^{-53}\le x/2$ を解くと $$ \begin{aligned} &\phantom{{}\iff{}} x-2\sqrt x+3\cdot 2^{-52} \ge 0 \\ &\iff (\sqrt x-1)^2-1+3\cdot 2^{-52} \ge 0 \\ &\iff (\sqrt x-1)^2 \ge 1-3\cdot 2^{-52} \\ &\iff |\sqrt x-1| \ge \sqrt{1-3\cdot 2^{-52}} \\ &\iff \sqrt x \ge 1+\sqrt{1-3\cdot 2^{-52}} \\ &\iff x \ge (1+\sqrt{1-3\cdot 2^{-52}})^2. \end{aligned} $$
$$ \begin{aligned} (1+\sqrt{1-3\cdot 2^{-52}})^2 &= 1 + 2\sqrt{1-3\cdot 2^{-52}} + (1-3\cdot 2^{-52}) \\ &= 2 + 2\sqrt{1-3\cdot 2^{-52}} - 3\cdot 2^{-52} \\ &\ge 2 + 2\cdot(1-2\cdot 2^{-52}) - 3\cdot 2^{-52} \\ &= 4 - 3.5\cdot 2^{-51} \\ \end{aligned} $$ なので、$x\ge 4-3\cdot 2^{-51}$ が必要*1。
Case 1: $x = 4-3\cdot 2^{-51}$
$$ \begin{aligned} \sqrt{4-3\cdot 2^{-51}}-3\cdot 2^{-53} &= 2\sqrt{1-3\cdot 2^{-53}} - 3\cdot 2^{-53} \\ &\gt 2\cdot(1-2\cdot 2^{-53}) - 3\cdot 2^{-53} \\ &= 2-3.5\cdot 2^{-52} \end{aligned} $$ より $y^{(i)}\ge 2-3\cdot 2^{-52}$ が必要。一方 $y^{(i)}\lt x/2 = 2-3\cdot 2^{-52}$ も必要なため、これを満たす $y^{(i)}$ は存在しない。
Case 2: $x = 4-2\cdot 2^{-51}$
$$ \begin{aligned} \sqrt{4-2\cdot 2^{-51}}-3\cdot 2^{-53} &= 2\sqrt{1-2\cdot 2^{-53}} - 3\cdot 2^{-53} \\ &\gt 2\cdot(1-\tfrac43\cdot 2^{-53}) - 3\cdot 2^{-53} \\ &\gt 2-2.9\cdot 2^{-52} \end{aligned} $$ より $y^{(i)}\ge 2-2\cdot 2^{-52}$ が必要。一方 $y^{(i)}\lt x/2 = 2-2\cdot 2^{-52}$ も必要なため、これを満たす $y^{(i)}$ は存在しない。
Case 3: $x = 4-2^{-51}$
$$ \begin{aligned} \sqrt{4-2^{-51}}-3\cdot 2^{-53} &= 2\sqrt{1-2^{-53}} - 3\cdot 2^{-53} \\ &\gt 2\cdot(1-\tfrac23\cdot 2^{-53}) - 3\cdot 2^{-53} \\ &\gt 2-2.2\cdot 2^{-52} \end{aligned} $$ より $y^{(i)}\ge 2-2\cdot 2^{-52}$ が必要。一方 $y^{(i)}\lt x/2 = 2-2^{-52}$ も必要なため、これを満たす $y^{(i)}$ は $2-2\cdot 2^{-52} = 2-2^{-51}$ しか存在しない。逆に、下記がすべて成り立つことは容易に確かめられる。 $$ \begin{aligned} 2-2^{-51} &\ge \sqrt{4-2^{-51}}-3\cdot 2^{-53}; \\ 2-2^{-51} &\le \sqrt{4-2^{-51}}; \\ 2-2^{-51} &\gt \tfrac14\cdot(4-2^{-51}); \\ 2-2^{-51} &\lt \tfrac12\cdot(4-2^{-51}). \quad\qed \end{aligned} $$
Lemma 12: Lemma 7 の条件に加えて $\varepsilon_i\in[-\delta_-\lldot 0]$ が成り立つとき、$|\varepsilon_{i+1}|\le 2^{-52}$ が成り立つ。
Proof
To-be-proved 1: $\varepsilon_{i+1}\ge -2^{-52}$
$\frac x{y^{(i)}}\in[\frac14\lldot 4)$ より、実数 $|\theta_{\oslash}|\le 2^{-52}$ が存在して $x\oslash y^{(i)} = (x\div y^{(i)})+\theta_{\oslash}$ が成り立つ。 また、Lemma 8 より、実数 $|\theta_{\oplus}|\le 2^{-52}$ が存在して $y^{(i)}\oplus(x\oslash y^{(i)}) = y^{(i)}+(x\oslash y^{(i)})+\theta_{\oplus}$ が成り立つ。 $$ \begin{aligned} y^{(i+1)} &= \frac12\cdot \left(y^{(i)} + \frac x{y^{(i)}} + \theta_{\oslash} + \theta_{\oplus}\right) \\ &= \sqrt x + \left(\frac{\varepsilon_i^2}{2(\sqrt x+\varepsilon_i)} + \frac{\theta_{\oslash} + \theta_{\oplus}}2\right) \end{aligned} $$ より、 $$ \begin{aligned} \varepsilon_{i+1} &= y^{(i+1)} - \sqrt x \\ &= \frac{\varepsilon_i^2}{2(\sqrt x+\varepsilon_i)} + \frac12\left(\theta_{\oslash} + \theta_{\oplus}\right) \\ &\ge \frac{0^2}{2(\sqrt x+\varepsilon_i)} - 2^{-52}. \end{aligned} $$
To-be-proved 2: $\varepsilon_{i+1}\le 2^{-52}$
Case 1: $\frac x{y^{(i)}}\in [\tfrac14\lldot 2]$
$|\theta_{\oslash}|\le 2^{-53}$ が成り立つため、 $$ \begin{aligned} \varepsilon_{i+1} &= \frac{\varepsilon_i^2}{2(\sqrt x+\varepsilon_i)} + \frac12\left(\theta_{\oslash} + \theta_{\oplus}\right) \\ &\le \frac{\delta_-^2}{2(1-\delta_-)} + 2^{-54} + 2^{-53} \\ &= \frac{9\cdot 2^{-106}}{2(1-3\cdot 2^{-53})} + 3\cdot 2^{-54} \\ &\lt 2^{-52}. \end{aligned} $$
Case 2: $\frac x{y^{(i)}}\in (2\lldot 4)$
Lemma 11 より、条件を満たす $(x, y^{(i)})$ は $(4-2^{-51}, 2-2^{-51})$ のみである。
$$ \begin{aligned} \frac x{y^{(i)}} &= \frac{4-2^{-51}}{2-2^{-51}} \\ &= \frac{4-2^{-50}+2^{-51}}{2-2^{-51}} \\ &= 2 + \frac{2^{-51}}{2-2^{-51}} \end{aligned} $$ $2^{-52} \lt \frac{2^{-51}}{2-2^{-51}} \lt 3\cdot 2^{-52}$ より、$x\oslash y^{(i)} = 2+2^{-51}$ となる。よって、 $$ y^{(i)} \oplus (x\oslash y^{(i)}) = \roundcirc{(2-2^{-51}) + (2+2^{-51})} = 4 $$ であり、$y^{(i+1)}=2$ となる。
$$ \begin{aligned} \varepsilon_{i+1} &= y^{(i+1)} - \sqrt x \\ &= 2 - \sqrt{4-2^{-51}} \\ &= 2 - 2\sqrt{1-2^{-53}} \\ &\lt 2 - 2\cdot(1-\tfrac23\cdot 2^{-53}) \\ &\lt 2^{-52}. \quad\qed \end{aligned} $$
Lemma 13: Lemma 7 の条件に加えて $\varepsilon_i\in[0\lldot \delta_+]$ が成り立つとき、$|\varepsilon_{i+1}|\le 2^{-52}$ が成り立つ。
Proof
$$ \begin{aligned} |(x\div y^{(i)}) - (x\oslash y^{(i)})| &\le 2^{-53}, \\ |(y^{(i)}+(x\oslash y^{(i)}) )-(y^{(i)}\oplus(x\oslash y^{i}) )| &\le 2^{-52} \end{aligned} $$ が成り立つので、三角不等式から $$ \left|\left(\sqrt x+\frac{\varepsilon_i^2}{2y^{(i)}}\right) - \left(0.5\otimes(y^{(i)}\oplus (x\oslash y^{(i)}) ) \right)\right| \le 2^{-53} + 2^{-54} $$ であり、 $$ \begin{aligned} |\sqrt x-y^{(i+1)}| &\le \frac{\varepsilon_i^2}{2y^{(i)}} + 3\cdot 2^{-54} \\ &\lt 2^{2\cdot (-26.5)-1} + 3\cdot 2^{-54} \\ &= 2^{-52}. \quad\qed \end{aligned} $$
Lemmata 8–13 をまとめると、$\varepsilon_i\in[-\delta_-\lldot \delta_+]$ ならば $|\varepsilon_{i+1}|\le 2^{-52}$ かつ $y^{(i+1)}\in [1\lldot 2]$ が成り立つということです。
正当性の証明
いよいよやっていきましょう。
Claim 14: $y^{(0)} = \iversonbracket{x\lt 2}\hat p(a^{[1]}, x) + \iversonbracket{x\ge 2}\hat p(a^{[2]}, x)$ とし、$\varepsilon_0 = y^{(0)}-\sqrt x$ とする。 このとき、$|\varepsilon_0|\lt 1.161\times 10^{-4}$ が成り立つ。
Proof
$(1+2^{-53})^3-1 \lt 3.331\times 10^{-16}$ が成り立つ。
Corollary 5 より、$x\in[1\lldot 2)$ の範囲で $$ \begin{aligned} |\hat p(a^{[1]}, x)-a^{[1]}(x)| &\le ( (1+2^{-53})^3-1)\cdot \sum_{i=0}^3 {|a^{[1]}_i\cdot x^i|} \\ &\lt (3.331\times 10^{-16}) \cdot \sum_{i=0}^3 {|a^{[1]}_i\cdot 2^i|} \\ &\lt (3.331\times 10^{-16}) \cdot 2.860 \\ &\lt 9.527 \times 10^{-16} \end{aligned} $$ が成り立つ。Claim 2 と三角不等式から $$ \begin{aligned} |\hat p(a^{[1]}, x)-\sqrt x| &\le |\hat p(a^{[1]}, x)-a^{[1]}(x)| + |a^{[1]}(x)-\sqrt x| \\ &\lt 9.527 \times 10^{-16} + 8.20594\times 10^{-5} \\ &\lt 8.206\times 10^{-5}. \end{aligned} $$
同様に、$x\in[2\lldot 4)$ の範囲で $$ \begin{aligned} |\hat p(a^{[2]}, x)-a^{[2]}(x)| &\le ( (1+2^{-53})^3-1)\cdot \sum_{i=0}^3 {|a^{[2]}_i\cdot x^i|} \\ &\lt (3.331\times 10^{-16}) \cdot \sum_{i=0}^3 {|a^{[2]}_i\cdot 4^i|} \\ &\lt (3.331\times 10^{-16}) \cdot 4.045 \\ &\lt 1.348 \times 10^{-15} \end{aligned} $$ が成り立つ。三角不等式から $$ \begin{aligned} |\hat p(a^{[2]}, x)-\sqrt x| &\le |\hat p(a^{[2]}, x)-a^{[2]}(x)| + |a^{[2]}(x)-\sqrt x| \\ &\lt 1.348\times 10^{-15} + 1.16050\times 10^{-4} \\ &\lt 1.161\times 10^{-4}. \quad\qed \end{aligned} $$
Claim 15: $y^{(1)} = 0.5\otimes(y^{(0)} \oplus (x\oslash y^{(0)}) )$ とし、$\varepsilon_1 = y^{(1)}-\sqrt x$ とする。 このとき、$\varepsilon_1\in[-\delta_-\lldot 6.740\times 10^{-9})$ が成り立つ。
Proof
Claim 14 から $|\varepsilon_0| \lt 1.161\times 10^{-4}$ なので、Lemma 7 より $\varepsilon_1 \in [-\delta_- \lldot \frac12\varepsilon_0^2+3\cdot 2^{-53}]$ が成り立つ。 $$ \begin{aligned} \tfrac12\varepsilon_0^2+3\cdot 2^{-53} &\lt 0.5\cdot (1.161\times 10^{-4})^2 + 3.331 \times 10^{-16} \\ &\lt 6.740\times 10^{-9}. \quad\qed \end{aligned} $$
Claim 16: $y^{(2)} = 0.5\otimes(y^{(1)} \oplus (x\oslash y^{(1)}) )$ とし、$\varepsilon_2 = y^{(2)}-\sqrt x$ とする。 このとき、$y^{(2)}\in[1\lldot 2]$ かつ $|\varepsilon_2|\lt 2^{-52}$ が成り立つ。
Proof
Claim 15 から、$\varepsilon_1 \in [-\delta_- \lldot 6.740\times 10^{-9})\subset[-\delta_-\lldot\delta_+]$ が成り立つ。よって、Lemmata 7–13 より従う。$\qed$
Lemma 17: $\roundcirc{\sqrt x} \in \{y^{(2)}\ominus 2^{-52}, y^{(2)}, y^{(2)}\oplus 2^{-52}\}$ が成り立つ。
Proof
Claim 16 より $y^{(2)}-2^{-52}\le \sqrt x\le y^{(2)}+2^{-52}$ である。
Case 1: $y^{(2)} = 1$
$x\in[1\lldot 4)$ であるから、$\sqrt x\in [1\lldot 2)\cap[1-2^{-52}\lldot 1+2^{-52}] = [1\lldot 1+2^{-52}]$ である。 よって、$\roundcirc{\sqrt x} \in\{1, 1+2^{-52}\}\subset\{y^{(2)}\ominus 2^{-52}, y^{(2)}, y^{(2)}\oplus 2^{-52}\}$ が成り立つ。
Case 2: $y^{(2)} = 2$
$2\oplus 2^{-52} = 2$ であるから、$\sqrt x\in\{2-2^{-52}, 2\}\subset\{y^{(2)}\ominus 2^{-52}, y^{(2)}, y^{(2)}\oplus 2^{-52}\}$ が成り立つ。
Case 3: $y^{(2)}\in(1\lldot 2)$
$y^{(2)}\ominus 2^{-52} = y^{(2)}-2^{-52}$ および $y^{(2)}\oplus 2^{-52} = y^{(2)}+2^{-52}$ であり、 $[1\lldot 2]$ の範囲の浮動小数点数は $2^{-52}$ の倍数であることから従う。$\qed$
さて、 $$ \begin{aligned} I_1 &= [y^{(2)}-2^{-52}\lldot y^{(2)}-2^{-53}), \\ I_2 &= (y^{(2)}-2^{-53}\lldot y^{(2)}+2^{-53}), \\ I_3 &= (y^{(2)}-2^{-53}\lldot y^{(2)}+2^{-52}] \end{aligned} $$ とすると、$\sqrt x\in I_1\sqcup I_2\sqcup I_3$ が成り立ちます。$\sqrt x \notin \{y^{(2)}-2^{-53}, y^{(2)}+2^{-53}\}$ は (1) の記事 の Claim 2 から従います。よって、下記が成り立ちます。 $$ \begin{aligned} \roundcirc{\sqrt x} &= \begin{cases} y^{(2)}\ominus 2^{-53}, & \text{if}\: y^{(2)}\in I_1; \\ y^{(2)}, & \text{if}\: y^{(2)}\in I_2; \\ y^{(2)}\oplus 2^{-53}, & \text{if}\: y^{(2)}\in I_3 \end{cases} \\ &= \begin{cases} y^{(2)}\ominus 2^{-53}, & \text{if}\: \sqrt x \lt y^{(2)}-2^{-53}; \\ y^{(2)}\oplus 2^{-53}, & \text{if}\: \sqrt x \gt y^{(2)}+2^{-53}; \\ y^{(2)}, & \text{otherwise} \end{cases} \\ &= \begin{cases} y^{(2)}\ominus 2^{-53}, & \text{if}\: x \lt (y^{(2)}-2^{-53})^2; \\ y^{(2)}\oplus 2^{-53}, & \text{if}\: x \gt (y^{(2)}+2^{-53})^2; \\ y^{(2)}, & \text{otherwise}. \end{cases} \\ \end{aligned} $$
Claim 18: 浮動小数点数 $x\in[1\lldot 4)$ および $y\in [1\lldot 2)$ に対して、$$(y+2^{-53})^2\lt x \iff \roundcirc{y\times(y\oplus 2^{-52}) + (-x)} \lt 0$$ が成り立つ。
Proof: (1) の記事 の Claims 2, 3 より従う。$\qed$
Claim 19: 浮動小数点数 $x\in[1\lldot 4)$ および $y\in [1\lldot 2)$ に対して、$$(y-2^{-53})^2\gt x \iff \roundcirc{y\times(y\ominus 2^{-52}) + (-x)} \ge 0$$ が成り立つ。
Proof
$y\in(1\lldot 2)$ のとき、Claim 18 の $y$ を $y-2^{-52}$ で置き換えればよい。
以下、$y=1$ について考える。 左辺について、$(1-2^{-53})^2\ge x$ かつ $x\ge 1$ なる $x$ は存在しない。 右辺について、 $$ \begin{aligned} 1\times(1\ominus 2^{-52}) + (-x) &= (1-2^{-52})-x \\ &\le -2^{-52} \end{aligned} $$ より、$\roundcirc{1\times(1\ominus 2^{-52}) + (-x)}\lt 0$ であり、こちらも条件を満たす $x$ は存在しない。$\qed$
Claim 20: 上記の手続きによって返される値は $\roundcirc{\sqrt x}$ に等しい。
Proof
Case 1: $\sqrt x\in I_1$
Claim 19 より従う。
Case 2: $\sqrt x\in I_2$
$y^{(2)}\lt 2$ の場合は Claims 18–19 より従う。$y^{(2)}=2$ の場合は $y^{(2)}=y^{(2)}\oplus 2^{-52}=2$ より従う。
Case 3: $\sqrt x\in I_3$
$y^{(2)}\lt 2$ の場合は Claim 18 より従う。$y^{(2)}=2$ の場合は $y^{(2)}=y^{(2)}\oplus 2^{-52}=2$ より従う。$\qed$
実装
実装
const A0: f64 = 0.37135166014697829073298862567753531038761138916015625; const A1: f64 = 0.78494263528793106754477548747672699391841888427734375; const A2: f64 = -0.180689144911217069999764817112009041011333465576171875; const A3: f64 = 0.0244769088312593204037614924573063035495579242706298828125; const B0: f64 = 0.52517055418962110824310229872935451567173004150390625; const B1: f64 = 0.55503826025453506520790369904716499149799346923828125; const B2: f64 = -0.06388325982676017200656559680282953195273876190185546875; const B3: f64 = 0.004326947054267089344536945105801351019181311130523681640625; const TWO_PM52: f64 = 2.220446049250313080847263336181640625e-16; fn sqrt(x: f64) -> f64 { if x < 0.0 { return f64::NAN; } if x == 0.0 || x.is_infinite() || x.is_nan() { return x; } let (x_m, x_e) = match frexp(x) { (m, e) if e % 2 == 0 => (4.0 * m, e - 2), (m, e) => (2.0 * m, e - 1), }; assert_eq!(x_e % 2, 0); assert!(1.0 <= x_m && x_m < 4.0); let mut y = if x_m < 2.0 { A3.mul_add(x_m, A2).mul_add(x_m, A1).mul_add(x_m, A0) } else { B3.mul_add(x_m, B2).mul_add(x_m, B1).mul_add(x_m, B0) }; y = 0.5 * (y + x_m / y); y = 0.5 * (y + x_m / y); y = if y.mul_add(y - TWO_PM52, -x_m) >= 0.0 { y - TWO_PM52 } else if y.mul_add(y + TWO_PM52, -x_m) < 0.0 { y + TWO_PM52 } else { y }; ldexp(y, x_e / 2) }
所感
疲れました。
不等式評価が下手すぎて一生かかったり、偽の命題を示そうとしていたり、示したと思ったら off-by-one を間違っていて破綻していたり、なんかもういろいろです。 とはいえ Newton 法でささっとよい感じにできるのを示せたので満足ではあります。
今回は、（証明のパートで使ったの以外は）整数型なしで計算を完結させられたので、ちょっとした達成感がありました。FMA を使うのは甘えかもしれません？（？？）
定義域が離散値なので整数型に逃げましたが、関数の最大値を求めるのってなんかよい方法がありますか？ 次数が小さければどうにでもなる気もしないでもないかもしれないですが、たとえば $7$ 次とかの多項式の特定の区間での最大値とかだとどうするんでしょう。二分法とかで「この（微小な）区間の中に極値があることはわかっているんだが」という状況までいったとして、そこからどうにかなるものですか？
（追記）該当の区間で凸であれば、端点での接線の交点を考えればよさそうな気がしました。
ところで rust-lang/libm/src/math/generic/sqrt.rs を見てみたところ、Newton–Raphson iterations ではなく Goldschmidt iterations というのが使われているらしいなぁとなりました。$1/\sqrt x$ も同時に計算できるらしいです？
さて、(1) から (3) まで平方根の話をしてさすがに飽きたので、そろそろ超越関数と仲よくなっていきたいです。 今回の平方根みたいなのだと「とりあえず Newton 法でばーっと詰めて細かいところはいい感じの式で合わせる」ということができましたが、超越関数だとそうもいかないような気がしています。「ばーっと詰める」で実際に近づいていることを示せたの自体は楽しかったですが、そういう手法が使える箱庭の中でしか生きていけない解法に過ぎないんだという感覚もあります。
今のところのそれっぽい知識が Lindemann–Weierstrass theorem くらいしかありません。 ちらっと本を読んだところだと Cody and Waite’s method というのがあるらしかったり、worst case の解析に連分数展開を使うらしかったり、わくわく感はあります。
おわり
おわりです。ごはんをたべましょう。
*1:$x\in[2\lldot 4)$ のとき $x$ は $2^{-51}$ の倍数であるため。
