倍精度浮動小数点型であるところの double ですが、現代においてはもはや double が「ふつう」になり、単精度と呼ばれている float は実質的に「半精度」のような感覚になっている気がします。
なにかしらを計算したいときに、計算結果は double だとしても計算途中は double よりも高い精度で行う必要がある場合はしばしばあります*1。
そうした場合に 64-bit や 113-bit の仮数部を持つ long double や __float128 などの型を扱える環境ならそれを使うのも手ですが、そうでない環境もしばしばありますし、将来 __float128 が「ふつう」になったときに「__float256 なしにはなにも打つ手がない...」となってしまいます。
内部実装として多倍長整数を使うような任意精度の型を使う方法もありますが、これも標準で用意されていない言語もしばしばあります*2し、速度面も気になってきます。
ともかく、今回からの数回は「ふつう」の浮動小数点数だけを使って、より高い精度で計算する方法について触れようと思います。
導入
一旦、円周率 $\pi$ について考えてみましょう。
$$
\pi = {\small 3.141592653589793}{\footnotesize 238462643383279}{\scriptsize 502884197169399}{\tiny 375105820974944592307{\dots}}.
$$
これを普通に double で表そうとすると以下のようになります。
$$
\begin{aligned}
\roundcirc{\pi} &= {\small 3.141592653589793}{\footnotesize 115997963468544}{\scriptsize 185161590576171}{\tiny 875} \\
&= 7074237752028440\times 2^{-51}.
\end{aligned}
$$
$\hat\pi_1 = \roundcirc{\pi}$ として、下記が成り立ちます。
$$
\pi - \hat\pi_1 = {\small 1.224646799147353}{\footnotesize 177226065932275}{\scriptsize 001058209749445}{\tiny 923078164062862{\dots}}\times 10^{-16}
$$
よって、$\hat\pi_2 = \roundcirc{\pi - \hat\pi_1}$ とすると下記のようになります。
$$
\begin{aligned}
\hat\pi_2 &= ({\small 1.224646799147353}{\footnotesize 207173764029458}{\scriptsize 396604625692124}{\tiny 677580063796256} \\
&\qquad\quad + {\tiny 12680683843791484832763671875}\times 10^{-60})\times 10^{-16} \\
&= 4967757600021511\times 2^{-105}.
\end{aligned}
$$
同様にして
$$
\begin{aligned}
\hat\pi_3 &= \roundcirc{\pi - (\hat\pi_1 + \hat\pi_2)} \\
&= -({\small 2.994769809718339}{\footnotesize 665887016354211}{\scriptsize 984016705500952}{\tiny 339037967798563} \\
& \qquad \quad + {\tiny 706161978508415962547495152434873233460166375152766704559326171875}\times 10^{-60}) \times 10^{-33} \\
&= 8753721960665020 \times 2^{-161}
\end{aligned}
$$
とすると、
$$
\pi - (\hat\pi_1 + \hat\pi_2 + \hat\pi_3) \approx 1.112\times 10^{-49}
$$
となり、$\hat\pi_1$ と比べて 3 倍程度の桁数で一致していることがわかります（3 倍程度の桁数を使っているのでそれはそう）。
ただし、もちろん $\hat\pi_1 \oplus \hat\pi_2 \oplus \hat\pi_3$ のように計算してしまうと $\hat\pi_1$ と同程度の精度に落ちてしまうため、内部表現としては $(\hat\pi_1, \hat\pi_2, \hat\pi_3)$ のように別個で持っておくものとします。 この $(\hat\pi_1, \hat\pi_2)$ や $(\hat\pi_1, \hat\pi_2, \hat\pi_3)$ のような表現を、それぞれ double-double, triple-double と呼んでいます。
上記の例では、「何らかの手段で別途与えた定数を 3 つの double に分解して表す」ということだけをしましたが、もちろんこれを計算に使う方法も必要になってきます。 すなわち、たとえば double と double-double の和や積を triple-double として得る（内部的には double の演算のみを用いる）方法などについても話していきます。
本題
今回は、double と double の和と積を double-double として表す方法について述べ、それ以外については次回以降にしたいと思います*3。
これを明らかと感じるかはわからないのですが、浮動小数点数 $a$, $b$ に対して $(a\oplus b)-(a+b)$ や $(a\otimes b)-(a\times b)$ は（underflow などが絡む議論は別途するとして、基本的に）同じフォーマットの浮動小数点数として表すことができます。 すなわち、double 同士の和・積の誤差を double で正確に表せるということです。よって、double 同士の和・積を誤差なく double-double として表せるということになります。
さて、「値としては表せる」という事実だけでは実際に計算できることにならないので、手続きについて説明していきます。もちろん、表せるという事実についてもちゃんと証明します。
記号に関して
任意の実数 $x \gt 0$ に対して、$2^k\le x\lt 2^{k+1}$ なる整数 $k$ が一意に存在するので、その $k$ を用いて $\hfloor{x} = 2^k$ と定義します。 また、便宜上 $\hfloor0 = 0$ としておきます。
任意の実数 $x\ge 0$ と整数 $k$ に対して明らかに $\hfloor{x\cdot 2^k} = \hfloor x\cdot 2^k$ が成り立ちます。
Property 1: 任意の実数 $a$, $b$ に対して $\hfloor{|a|}\ge \hfloor{|b|}$ ならば $\hfloor{|a+b|}\le 2\hfloor{|a|}$ が成り立つ。
Proof
Case 1: $\hfloor{|a|} = \hfloor{|b|} = 0$
このとき、$a = b = 0$ であり、$\hfloor{|a+b|} = 2\hfloor{|a|} = 0$ である。
Case 2: $\hfloor{|a|} = \hfloor{|b|} \gt 0$
整数 $k$ と実数 $a', b' \in (-2\lldot -1]\sqcup[1\lldot 2)$ が一意に存在して $$ \begin{aligned} 2^k &= \hfloor{|a|} = \hfloor{|b|}, \\ a &= a'\cdot 2^k, \\ b &= b'\cdot 2^k \end{aligned} $$ が成り立つ。ここで $|a'+b'|\lt 4$ であり、$\hfloor{|a'+b'|}\le 2$ となる。 $$ \begin{aligned} \hfloor{|a+b|} &= \hfloor{|(a'+b')\cdot 2^k|} \\ &= \hfloor{|a'+b'|}\cdot 2^k \\ &\le 2\cdot 2^k \\ &= 2\hfloor{|a|}. \end{aligned} $$
Case 3: $\hfloor{|a|} \gt \hfloor{|b|}$
このとき、$|a| \gt |b|$ が成り立つ。 $$ \begin{aligned} \hfloor{|a+b|} &\le \hfloor{|a| + |b|} \\ &\le \hfloor{|a| + |a|} \\ &= \hfloor{2|a|} \\ &= 2\hfloor{|a|}. \quad\qed \end{aligned} $$
Property 2: 任意の実数 $a, b\ge 0$ に対して、$a\le b$ ならば $\hfloor a\le \hfloor b$ が成り立つ。
Proof
このとき、$a\le b$ より $\hfloor a\lt 2\hfloor b$ である。すなわち、$\hfloor{a}\le \hfloor b$ が成り立つ。$\qed$
正負のゼロについて、0.0 を $0_+$、-0.0 を $0_-$ と表記します。以下では単に $0$ と書いた場合は $0_+$ を意味するものとします。
特に $0_-$ について考えなくていいかなと思っている主な理由は下記の通りです。
- $\roundcirc0 = 0_+$ であるため
- 最終的に実装したい関数において、入力が $0_-$ だった場合は別途処理すればよいため
- double-double, triple-double の計算途中で $0_-$ が出てこないのが期待されるため
- 計算途中で $0_-$ が出てくる場合でも、$\pm\infty$ や $0_{\pm}$ 以外との演算では符号は影響しないため
さて、double として考えている値たちを改めて定式化しておきます。
$-1022\le k\le 1023$ に対して、下記を定義します。それぞれ符号ごとの正規化数の集合です。 $$ \begin{aligned} F_k^- &= \{-m\cdot 2^{k-52} \mid m\in [2^{52}\lldot 2^{53})\cap\Z\}, \\ F_k^+ &= \{+m\cdot 2^{k-52} \mid m\in [2^{52}\lldot 2^{53})\cap\Z\}. \end{aligned} $$ また、非正規化数の集合を下記で定義します*4。 $$ \begin{aligned} F_{-1023}^- &= \{-m\cdot 2^{-1022-52} \mid m\in[1\lldot 2^{52})\cap\Z\}, \\ F_{-1023}^+ &= \{+m\cdot 2^{-1022-52} \mid m\in[1\lldot 2^{52})\cap\Z\}. \end{aligned} $$ また、考えたい数全体からなる集合 $F$ を下記で定義します。 $$ F = \left(\bigsqcup_{k=-1023}^{1023} {(F_k^-\sqcup F_k^+)}\right) \sqcup \{-\infty, 0, +\infty\}. $$ たとえば下記が成り立ちます。 $$ \begin{aligned} 1 &\in F_0^+, \\ -2^{-1074} &\in F_{-1023}^-, \\ 2^{1024} &\notin F, \\ 1+2^{-53} &\notin F, \\ 0.1 &\notin F. \end{aligned} $$ さらに、下記の性質も成り立ちます。
- $(F\smallsetminus \{-\infty, +\infty\})\subset \Q$
- 任意の実数 $x$ に対して $x\in F \iff -x \in F$
- 任意の実数 $x$ に対して $x\in F \implies x\equiv 0\pmod{2^{-1074}}$
- 任意の $x, y \in (F_{-1023}^-\sqcup F_{-1023}^+)$ に対して $x\pm y\in F$
Lemma 3: 任意の $k\in[-1022\lldot 1023]\cap\Z$ および $m\in[1\lldot 2^{53})\cap\N$ に対して、$$\{-m\cdot 2^{k-52}, +m\cdot 2^{k-52}\}\subset F$$が成り立つ。
Proof
任意の $k\in[-1022\lldot 1023]$ および $m\in[1\lldot 2^{53})\cap\N$ に対して $m\cdot 2^{k-52} \in F$ を示せば十分。下記で $P(k)$ を定義し、$k$ に関する数学的帰納法で示す。 $$ P(k) \iff \Forall{m\in[1\lldot 2^{53})\cap\N}{m\cdot 2^{k-52}\in F}. $$
To-be-proved 1: $P(-1022)$
任意の $m\in[1\lldot 2^{-53})\cap\N$ に対し、$m\cdot 2^{k-52} \in F_{-1023}^+\sqcup F_{-1022}^+$ が成り立つ。
To-be-proved 2: $k\lt 1023 \wedge P(k) \implies P(k+1)$
$m\in[2^{52}\lldot 2^{53})\cap\N$ のとき $m\cdot 2^{(k+1)-52}\in F_{k+1}^+$ が成り立つ。
$m\in[1\lldot 2^{52})\cap\N$ のとき、$2m\in[2\lldot 2^{53})\cap\N$ かつ $m\cdot 2^{(k+1)-52} = (2m)\cdot 2^{k-52}$ が成り立ち、$P(k)$ より従う。$\qed$
tiesToEven による丸めについても定義しましょう。$\circ\colon \R\to F$ を次のように定義します。 ただし、閾値として $\Theta_{\infty} = (2-2^{-53})\cdot 2^{1023}$ および $\Theta_0 = 2^{-53}\cdot 2^{-1022}$ としておきます。 $$ \roundcirc{x} = \begin{cases} +\infty, & \text{if}\: x \ge \Theta_{\infty}; \\ \roundcirc{x\cdot\hfloor{x}^{-1}\cdot 2^{52}}\cdot\hfloor{x}\cdot 2^{-52}, & \text{if}\: x\in[2^{53}\lldot\Theta_{\infty}); \\ \floor{x}, & \text{if}\: x\in[2^{52}\lldot 2^{53}) \wedge (x\bmod 1)\lt \tfrac12; \\ \floor{\frac{2x+1}4}\cdot 2, & \text{if}\: x\in[2^{52}\lldot 2^{53}) \wedge (x\bmod 1) = \tfrac12; \\ \ceil{x}, & \text{if}\: x\in[2^{52}\lldot 2^{53}) \wedge (x\bmod 1)\gt \tfrac12; \\ \roundcirc{x\cdot\hfloor{x}^{-1}\cdot 2^{52}}\cdot\hfloor{x}\cdot 2^{-52}, & \text{if}\: x\in[2^{-1022}\lldot 2^{52}); \\ \roundcirc{x+2^{-1022}}-2^{-1022}, & \text{if}\: x\in(\Theta_0\lldot 2^{-1022}); \\ 0, &\text{if}\: x\in[0\lldot \Theta_0]; \\ \roundcirc{-x}, & \text{if}\: x\lt 0. \end{cases} $$ 特に、$x\in F\iff \roundcirc x=x$ です。また $x\bmod 1=\tfrac12$ のとき $\floor{\frac{2x+1}4}\cdot 2 \in \{\floor x, \ceil x\}$ です。
また、$(x, y)\in \R\times F$ かつ $y\le x$ ならば $y\le \roundcirc x$（$\le$ を $\ge$ で置き換えた命題も同様）です。
任意の実数 $x$ に対して、$|\roundcirc x|\lt \infty$ ならば下記が成り立ちます。 $$ |\roundcirc x-x| \le \begin{cases} \hfloor{|x|}\cdot 2^{-53}, & \text{if}\: |x|\ge 2^{-1022}; \\ 2^{-1022-53}, & \text{if}\: |x|\lt 2^{-1022}. \end{cases} $$
さらに、任意の実数 $x$ と整数 $k$ に対して $x\equiv 0\pmod{2^k}$ ならば $\roundcirc x\equiv 0\pmod{2^k}$ が成り立ちます。
補題たち
頻出なので、改めて触れておきます。
Lemma 4 (Sterbenz): $x, y\in F\smallsetminus{\{-\infty, +\infty\}}$ について、$\frac x2\le y\le 2x$ ならば $x\ominus y=x-y$ が成り立つ。
Proof
$x=y$ のときは明らか（特に $x=y=0$ も含む）。一般性を失わず $x\gt y\gt 0$ とする*5。また、$x-y\le \frac x2$ に注意する。
Case 1: $y\in F_{-1023}^+$
$x\le 2y$ より、$x\in F_{-1023}^+\sqcup F_{-1022}^+$ である。
Case 1-1: $x\in F_{-1023}^+$
ある整数 $m_x$, $m_y$ ($0\lt m_y\lt m_x\lt 2^{52}$) に対して $x = m_x\cdot 2^{-1074}$ および $y = m_y\cdot 2^{-1074}$ が成り立つ。 $$ x - y = (m_x-m_y)\cdot 2^{-1074} $$ であり、$m_x-m_y\in[1\lldot 2^{52})\cap\N$ であるから、$x-y\in F_{-1023}^+$ である。
Case 1-2: $x\in F_{-1022}^+$
ある整数 $2^{52}\le m_x\lt 2^{53}$, $1\le m_y\lt 2^{52}$ に対して $x = m_x\cdot 2^{-1074}$ および $y = m_y\cdot 2^{-1074}$ が成り立つ。 $$ x - y = (m_x-m_y)\cdot 2^{-1074} $$ である。ここで、$m_x-m_y\in[1\lldot 2^{53})\cap\N$ であるから、$x-y\in F_{-1023}^+\sqcup F_{-1022}^+$ である。
Case 2: ある整数 $-1022\le k\le 1023$ に対して $y\in F_k^+$
Case 1 同様に $x\in F_k^+\sqcup F_{k+1}^+$ である。
Case 2-1: $x\in F_k^+$
ある整数 $2^{52}\le m_y\lt m_x\lt 2^{53}$ に対して $x = m_x\cdot 2^{k-52}$ および $y = m_y\cdot 2^{k-52}$ が成り立つ。 $$ x - y = (m_x-m_y)\cdot 2^{k-52} $$ であり、かつ $m_x-m_y\in[1\lldot 2^{52})\cap\N$ であるから、Lemma 3 より $x-y\in F$ が従う。
Case 2-2: $x\in F_{k+1}^+$
ある整数 $2^{52}\le m_x\lt 2^{53}$ および $2^{52}\le m_y\lt 2^{53}$ に対して $x = m_x\cdot 2^{k+1-52}$ および $y = m_y\cdot 2^{k-52}$ が成り立つ。 $$ \begin{aligned} x-y &= m_x\cdot 2^{k+1-52} - m_y\cdot 2^{k-52} \\ &= (2m_x-m_y)\cdot 2^{k-52} \end{aligned} $$ と書ける。$x-y\le\tfrac x2$ より、$2m_x-m_y\le m_x\lt 2^{53}$ が従う。Case 2-1 同様にして、$x-y\in F$ が従う。$\qed$
Lemma 5: 整数 $m\in[1\lldot 2^{106})$ および $k\in[-1022\lldot 1023-53]$ に対して $x = m\cdot 2^{k-52}$ と表せるとき、ある $(x_h, x_l)\in F^2$ が存在して $x=x_h+x_l$ が成り立つ。
Proof
$m_h = \floor{m/{2^{53}}}$ および $m_l = (m\bmod 2^{53})$ とすると $m = m_h\cdot 2^{53} + m_l$ が成り立つ。 $$ \begin{aligned} x_h &= (m_h\cdot 2^{53}) \cdot 2^{k-52} \\ &= m_h\cdot 2^{(k+53)-52}, \\ x_l &= m_l\cdot 2^{k-52} \end{aligned} $$ とすると、$(x_h, x_l)\in F^2$ かつ $x = x_h + x_l$ が成り立つ。$\qed$
$x, y\in F\smallsetminus{\{-\infty, +\infty\}}$ に対して、$(x\oplus y)-(x+y)$ すなわち $x\oplus y$ の誤差について考えます。
Lemma 6: $x, y\in F$ が $x\gt 0$ かつ $\hfloor x\ge \hfloor{|y|}$ かつ $x+y\lt \Theta_{\infty}$ を満たすとする。このとき、$(x\oplus y)-x \in F$ かつ $(x\oplus y)-(x+y)\in F$ が成り立つ。
Proof
Case 1: $\hfloor{|y|} \le \hfloor x\cdot 2^{-55}$
このとき $|y|\lt \hfloor x\cdot 2^{-54}$ であり、$x\oplus y=x$ であるから、$(x\oplus y)-x = 0$ かつ $(x\oplus y)-(x+y)=-y$ である。 明らかに $0, -y\in F$ である。
Case 2: $\hfloor{|y|} = \hfloor x\cdot 2^{-54}$
このとき $\hfloor x\cdot 2^{-54}\le |y|\lt \hfloor x\cdot 2^{-53}$ である。また、$x\in F$ より $\hfloor{x}\le 2^{1023}$ であるから、$\hfloor{|y|}\le 2^{1023-54}$ である。また、$x\notin F_{-1023}^+$ である。
Case 2-1: $x = \hfloor x$ かつ $-\hfloor x\cdot 2^{-54}\lt y\lt -\hfloor x\cdot 2^{-53}$
このとき $x\oplus y = x-\hfloor x\cdot 2^{-53}$ であるから、 $$ \begin{aligned} (x\oplus y)-x &= -\hfloor x\cdot 2^{-53} \\ &= -2\hfloor{|y|}, \\ (x\oplus y)-(x+y) &= -y-2\hfloor{|y|} \\ &= |y| - 2\hfloor{|y|} \end{aligned} $$ である。$\hfloor{|y|}\le 2^{969}$ より $2\hfloor{|y|}\in F$ であり、$|y|-2\hfloor{|y|}\in F$ も Sterbenz lemma から従う。
Case 2-2: $x = \hfloor x$ かつ $y = -\hfloor x\cdot 2^{-54}$
$x\oplus y = x$ であるから、Case 1 同様にして従う。
Case 2-3: $x = \hfloor x$ かつ $y\gt 0$
$x\oplus y = x$ であるから、Case 1 同様にして従う。
Case 2-4: $x \gt \hfloor x$
$x\oplus y = x$ であるから、Case 1 同様にして従う。
Case 3: $\hfloor{|y|} = \hfloor x\cdot 2^{-53}$
このとき $\hfloor x\cdot 2^{-53}\le |y|\lt \hfloor x\cdot 2^{-52}$ である。
$x=\hfloor x$ のとき、$y\lt 0$ ならば $x\oplus y\in\{x-\hfloor x\cdot 2^{-52}, x-\hfloor x\cdot 2^{-53}\}$ であり、$y\gt 0$ ならば $x\oplus y = x$ である。
$x\gt \hfloor x$ のとき、$y\lt 0$ ならば $x\oplus y\in\{x-\hfloor x\cdot 2^{-52}, x\}$ であり、$y\gt 0$ ならば $x\oplus y\in\{x, x+\hfloor x\cdot 2^{-52}\}$ である。
Sterbenz lemma に気をつけつつ、総括して、$y\lt 0$ ならば $$ \begin{aligned} (x\oplus y)-(x+y) &\in \{-y-\hfloor x\cdot 2^{-52}, -y-\hfloor x\cdot 2^{-53}, -y\} \\ &= \{-y-2\hfloor{|y|}, -y-\hfloor{|y|}, -y\} \\ &= \{|y|-2\hfloor{|y|}, |y|-\hfloor{|y|}, |y|\} \subset F \end{aligned} $$ であり、$y\gt 0$ ならば $$ \begin{aligned} (x\oplus y)-(x+y) &\in \{-y, -y+\hfloor x\cdot 2^{-52}\} \\ &= \{-y, -y+2\hfloor{|y|}\} \\ &= \{-y, -y+2\hfloor{y}\} \subset F \end{aligned} $$ である。
Case 4: $\hfloor x\cdot 2^{-52} \le \hfloor{|y|} \le \hfloor x$
$|y|\gt 0$ に注意する。
Case 4-1: $y\lt 0$
Case 4-1-1: $y\in F_{-1023}^-$ かつ $x\in F_{-1023}^+$
ある整数 $m_x, m_y\in[1\lldot 2^{52})$ に対して $x = m_x\cdot 2^{-1022-52}$ かつ $y = -m_y\cdot 2^{-1022-52}$ が成り立つ。 $$ x + y = (m_x-m_y)\cdot 2^{-1022-52} $$ であり、$|m_x-m_y| \in [0\lldot 2^{53})\cap\N$ より $x+y = x\oplus y$ が従う。 よって、$(x\oplus y)-x = y$ かつ $(x\oplus y)-(x+y) = 0$ が成り立つ。
Case 4-1-2: $y\in F_{-1023}^-$ かつ、ある整数 $k\ge -1022$ に対して $x\in F_k^+$
$$ \begin{aligned} \hfloor x &\le \hfloor{|y|}\cdot 2^{52} \\ &\lt 2^{52}\cdot 2^{-1022-52}\cdot 2^{52} \\ &= 2^{-1022+52} \end{aligned} $$ であるから、$x\lt 2^{-1022+52}$ となる。また、$x+y\equiv 0\pmod{2^{-1022-52}}$ が成り立つ。
$x+y\lt 2^{-1022}$ のとき $x+y\in F_{-1023}^+$ より $x\oplus y=x+y$ が従うため、以降は $x+y\ge 2^{-1022}$ とする。 このとき、 $$ \begin{aligned} |(x\oplus y)-(x+y)| &\le (x+y)\cdot 2^{-53} \\ &\lt x\cdot 2^{-53} \\ &\le 2^{-1022+52}\cdot 2^{-53} \\ &= 2^{-1022-1} \end{aligned} $$ より $(x\oplus y)-(x+y)\in F_{-1023}^+$ が成り立つ。 また、$y\in F_{-1023}^-$ より $(x\oplus y)-x = ( (x\oplus y)-(x+y) ) + y \in F$ が従う。
Case 4-1-3: ある整数 $k\ge -1022$ に対して $y\in F_k^-$
このとき、ある整数 $k'\ge k$ に対して $x\in F_{k'}^+$ である。Case 4-1-2 同様 $x+y\ge 2^{-1022}$ の場合について考える。
また、$\frac x2\le |y|\le 2x$ のときは Sterbenz lemma より $x\oplus y=x+y$ が従うから、$-y \lt \frac x2$ とする。 このとき $\tfrac x2\lt x+y$ である。よって $\tfrac x2\le x\oplus y$ であるから、Sterbenz lemma より $(x\oplus y)-x\in F$ である。
さて、 $$ \begin{aligned} \hfloor x &\le \hfloor{|y|}\cdot 2^{52} \\ &= 2^{52}\cdot 2^{k-52}\cdot 2^{52} \\ &= 2^{k+52} \end{aligned} $$ より $x\lt 2^{k+53}$ である。 $x+y\equiv 0 \pmod{2^{k-52}}$ から、ある整数 $0\le m\lt 2^{105}$ が存在して $x+y = m\cdot 2^{k-52}$ が成り立つ。
$x\oplus y\equiv 0\pmod{2^{k-52}}$ かつ $|(x\oplus y)-(x+y)| \le (x+y)\cdot 2^{-53}$ であるから、ある整数 $0\le m'\lt 2^{52}$ が存在して $(x\oplus y)-(x+y) = m'\cdot 2^{k-52}$ が成り立つため、$(x\oplus y)-(x+y)\in F$ が従う。
Case 4-2: $y\gt 0$
Case 4-2-1: $y\in F_{-1023}^+$ かつ $x\in F_{-1023}^+$
Case 4-1-1 同様にして従う。
Case 4-2-2: ある整数 $k\ge -1022$ に対して $x\in F_k^+$
このとき、$x+y\ge 2^{-1022}$ が成り立つ。
$|(x\oplus y)-(x+y)|\le (x+y)\cdot 2^{-53}$ より Sterbenz lemma が従い、$(x\oplus y)-(x+y)\in F$ となる。
$x+y\le 2x$ のとき $x\oplus y\le 2x$ であるから、Sterbenz lemma より $(x\oplus y)-x\in F$ が従う。 以降はそうでない場合を考える。すなわち $y\gt x$ とする。$\hfloor x=\hfloor y$ かつ $\hfloor{x+y} = 2\hfloor x$ が成り立つ。
よって、ある $m_x, m_y\in [2^{52}\lldot 2^{53})\cap\N$ に対して $x = m_x\cdot 2^{k-52}$ かつ $y = m_y\cdot 2^{k-52}$ が成り立つ。 $$ x\oplus y \in \{(m_x+m_y-1)\cdot 2^{k-52}, (m_x+m_y)\cdot 2^{k-52}, (m_x+m_y+1)\cdot 2^{k-52}\} $$ であり、 $$ (x\oplus y)-x \in \{(m_y-1)\cdot 2^{k-52}, m_y\cdot 2^{k-52}, (m_y+1)\cdot 2^{k-52}\} \subset F_k^+ $$ が成り立つ。$\qed$
Error-free transformation
さて、実際にやっていきます。$a+b$ や $a\times b$ などに対して、二つの浮動小数点数 $h$, $l$ であって（誤差なく）$h+l$ と表せるものを求めます。
今後の応用を見据えると $a\oplus b$ や $a\otimes b$ が非正規化数になる場合は除外してもいい気がしますが、非正規化数ちゃんを仲間外れにして忌避感を強めるのも嫌なので、できる限りケアしていきます。
和
最も基本的な操作ですが、Lemma 6 相当の前提がないと証明は面倒です。
- 入力: $(a, b)\in F\times F$
- 出力: $(s_h, s_l)\in F^2$
- 事前条件:
- $|a+b| \lt \Theta_{\infty}$, and
- $\hfloor{|a|}\ge \hfloor{|b|}$
- 事後条件:
- $s_h = a\oplus b$, and
- $s_h + s_l = a + b$
- 手続き:
- $s_h \gets a\oplus b$ で初期化する。
- $z \gets s_h\ominus a$ で初期化する。
- $s_l \gets b\ominus z$ で初期化する。
- $(s_h, s_l)$ を出力する。
Claim 7: 上記の手続きの実行後、事後条件が成り立つ。
Proof
$s_h = a\oplus b$ は明らかなので、$s_h+s_l = a+b$ を示す。
Case 1: $a\lt 0$
$a\gt 0$ で成り立つならば $a\lt 0$ でも成り立つことを示す。 $$ \begin{aligned} s_h^{\pm} &= (\pm a)\oplus (\pm b), \\ z^{\pm} &= s_h^{\pm}\ominus (\pm a), \\ s_l^{\pm} &= (\pm b)\ominus z^{\pm} \end{aligned} $$ とする（複号同順）。このとき、 $$ \begin{aligned} s_h^- &= (-a)\oplus(-b) \\ &= -(a\oplus b) \\ &= -s_h^+, \\ z^- &= s_h^- \ominus (-a) \\ &= -(s_h^+ \ominus a) \\ &= -z^+, \\ s_l^- &= (-b)\ominus z^- \\ &= -(b\ominus z^+) \\ &= -s_l^+ \end{aligned} $$ であるから、$s_h^- + s_l^- = -(s_h^+ + s_l^+)$ が成り立つ。 よって、$s_h^+ + s_l^+ = a+b$ ならば $s_h^- + s_l^- = -(a + b) = (-a) + (-b)$ が成り立つ。
Case 2: $a = 0$
$$ \begin{aligned} s_h &= 0\oplus b = b, \\ z &= b\ominus 0 = b, \\ s_l &= b\ominus b = 0 \end{aligned} $$ より、$s_h + s_l = b = a + b$ が成り立つ。
Case 3: $a\gt 0$
Lemma 6 より、$z = s_h\ominus a = s_h - a$ が従う。また、 $$ \begin{aligned} s_l &= b\ominus z \\ &= b - (s_h - a) \\ &= (a+b)-s_h \\ &= (a+b)-(a\oplus b) \end{aligned} $$ である。$\qed$
$\hfloor{|a|}\ge \hfloor{|b|}$ が事前に保証できない場合は $|a|\ge |b|$ で分岐すればよいです。 わざわざ $\hfloor{\bullet}$ の部分を実装したりする必要はありません。
$a$, $b$ の大小によらない亜種なども考案されていますが、ここでは割愛します。
丸め
次以降で使うサブルーチンです。$\textsc{Round}(a, k)$ として使います。 「$a$ の仮数部を $53-k$ bits に丸めた値」とその際の誤差に分離します。
- 入力: $(a, k)\in F \times \N$
- 出力: $(a_h, a_l)\in F^2$
- 事前条件:
- $|a| \le 2^{1023-k}$, and
- $k\in[1\lldot 52]$
- 事後条件:
- $a_h + a_l = a$,
- $a_h \equiv 0 \pmod{\hfloor{|a|}\cdot 2^{-(52-k)}}$, and
- $|a_l| \le \hfloor{|a|}\cdot 2^{-(53-k)}$.
- 手続き:
- $c \gets 2^k+1$ で初期化する。
- $a_c \gets a\otimes c$ で初期化する。
- $a_h \gets (a\ominus a_c)\oplus a_c$ で初期化する。
- $a_l \gets a \ominus a_h$ で初期化する。
- $(a_h, a_l)$ を出力する。
Claim 8: 上記の手続きの実行後、事後条件が成り立つ。
Proof
$a = 0$ のときは明らか。$a\gt 0$ で成り立てば $a\lt 0$ で成り立つことも明らかなので、以降 $a\gt 0$ とする。
Case 1: $a\ge 2^{-1022}$
$a\cdot (2^k+1) \ge 2^{-1022}$ に注意し、$\hfloor{(2^k+1)\cdot a} = 2^{52}$ とする。 そうでない場合は $a$ を $\hfloor{(2^k+1)\cdot a}^{-1}\cdot 2^{52}$ 倍することで帰着できる。 特に、$|a\cdot(2^k+1)| \le 2^{1023-1}\cdot (2^1+1)\lt \Theta_{\infty}$ に注意する。
まず $2^{51}+\tfrac12 \le 2^k\cdot a$ を示す。$2^{51}\le \tfrac12(2^k+1)\cdot a$ であるから、$\tfrac12(2^k+1)\cdot a\le 2^k\cdot a-\tfrac12$ を示せば十分。 $$ \begin{aligned} &\phantom{{}\iff{}} \tfrac12(2^k+1)\cdot a\le 2^k\cdot a-\tfrac12 \\ % &\iff \tfrac12\cdot 2^k\cdot a + \tfrac12\cdot a \le 2^k\cdot a-\tfrac12 \\ &\iff \tfrac12\cdot a \le \tfrac12\cdot 2^k\cdot a-\tfrac12 \\ &\iff a\le 2^k\cdot a-1 \\ &\iff 1\le (2^k-1)\cdot a \\ &\iff 1\le \tfrac{2^k-1}{2^k+1}\cdot (2^k+1)\cdot a \end{aligned} $$ $k\ge 1$ のとき $\tfrac{2^k-1}{2^k+1}\ge \tfrac13$ であり、$(2^k+1)\cdot a\ge 2^{52}$ より従う。
さて、 $$ \begin{aligned} (2^k+1)\cdot a &= 2^k\cdot a + a \\ &= \floor{2^k\cdot a + a} + ( (2^k\cdot a + a)\bmod 1) \end{aligned} $$ である。
$\hfloor{2^k\cdot a} \ge 2^{51}$ かつ $2^k\cdot a\in F$ であるから、$( (2^k\cdot a)\bmod 1)\in\{0, \tfrac12\}$ である。
また、$\hfloor a\ge 2^{51-k}$ であるから、$a\equiv 0\pmod{2^{-k-1}}$ が成り立つ。特に、$a\equiv 0\pmod{2^{-53}}$ である。
Case 1-1: $( (2^k\cdot a)\bmod 1) = 0$
$$ \begin{aligned} (2^k+1)\cdot a &= 2^k\cdot a + \floor{a} + (a\bmod 1). \end{aligned} $$
$r = (a\bmod 1)$ とすると、下記が成り立つ。
- $r\le \tfrac12 \wedge (2^k+1)\otimes a = 2^k\cdot a + \floor{a}$, or
- $r\ge \tfrac12 \wedge (2^k+1)\otimes a = 2^k\cdot a + \floor{a} + 1$.
$r\le \tfrac12$ のとき、 $$ \begin{aligned} a - ( (2^k+1)\otimes a) &= (\floor{a} + r) - (2^k\cdot a + \floor{a}) \\ &= -(2^k\cdot a - r) \end{aligned} $$ である。また、$r\ge \tfrac12$ のとき、 $$ \begin{aligned} a - ( (2^k+1)\otimes a) &= (\floor{a} + r) - (2^k\cdot a + \floor{a} + 1) \\ &= -(2^k\cdot a - (r-1) ) \end{aligned} $$ である。
Case 1-1-1: $2^k\cdot a\gt 2^{52}$
remark: $( (2^k\cdot a)\bmod 1) = 0$ より $2^k\cdot a-1\ge 2^{52}$ が成り立つ。
$r\le \tfrac12$ のとき、ある $r'\in\{0, 1\}$ に対して $a\ominus ( (2^k+1)\otimes a) = -(2^k\cdot a-r')$ かつ $|r-r'|\le \tfrac12$ が成り立つ。 よって、 $$ \begin{aligned} &\phantom{{}={}} (a\ominus ( (2^k+1)\otimes a) ) + ( (2^k+1)\otimes a) \\ &= -(2^k\cdot a-r') + (2^k\cdot a + \floor a) \\ &= \floor a + r' \end{aligned} $$ となる。$\hfloor a=2^{52-k}$ より $a\le 2^{52}-\frac12$ であるから、$|\floor a + r'|\le 2^{52}$ であり、$\floor a+r'\in F$ である。 よって $$ (a\ominus ( (2^k+1)\otimes a) ) \oplus ( (2^k+1)\otimes a) = \floor a+r' $$ となる。$\floor a+r'\equiv 0\pmod 1$ かつ $1 = \hfloor{a}\cdot 2^{-(52-k)}$ である。
よって $a_h \equiv 0 \pmod{\hfloor a\cdot 2^{-(52-k)}}$ である。また、$a_l = a-(\floor a+r') = r-r'$ であるから、$|a_l|\le \hfloor a\cdot2^{-(53-k)}$ である。
また、$r\ge \tfrac12$ のとき、ある $r'\in\{0, 1\}$ に対して $a\ominus ( (2^k+1)\otimes a) = -(2^k\cdot a-r')$ かつ $|(r-1)-r'|\le \tfrac12$ が成り立ち、$r\le \tfrac12$ の場合と同様にして従う。特に、 $$ \begin{aligned} a_l &= a - a_h \\ &= (\floor a+r) - (\floor a+r'+1) \\ &= (r-1)-r' \end{aligned} $$ である。
Case 1-1-2: $2^k\cdot a\le 2^{52}$
ある $r'\in\{0, \tfrac12, 1\}$ に対して $a\ominus ( (2^k+1)\otimes a) = -(2^k\cdot a-r')$ が成り立つ。 $r\le \tfrac12$ のとき $|r-r'|\le \tfrac14$、$r\ge \tfrac12$ のとき $|(r-1)-r'|\le \tfrac14$ が成り立つ。
よって、Case 1-1-1 同様にして $$ (a\ominus ( (2^k+1)\otimes a) ) \oplus ( (2^k+1)\otimes a) = \floor a+r' $$ となる。$\floor a+r'\equiv 0\pmod{\tfrac12}$ かつ $\tfrac12 = \hfloor{a}\cdot 2^{-(52-k)}$ である。 すなわち、$a_h\equiv 0\pmod{\hfloor a\cdot 2^{-(52-k)}}$ かつ $|a_l|\le \hfloor a\cdot 2^{-(53-k)}$ が成り立つ。
Case 1-2: $( (2^k\cdot a)\bmod 1) = \tfrac12$
このとき、$\hfloor{2^k\cdot a} = 2^{51}$ である。 $$ \begin{aligned} (2^k+1)\cdot a &= 2^k\cdot a - \tfrac12 + \floor{\tfrac12+a} + ( (\tfrac12+a)\bmod 1). \end{aligned} $$ $r = ( (\tfrac12+a)\bmod 1)$ とすると、下記が成り立つ。
- $r\le \tfrac12 \wedge (2^k+1)\otimes a = 2^k\cdot a - \tfrac12 + \floor{\tfrac12+a}$, or
- $r\ge \tfrac12 \wedge (2^k+1)\otimes a = 2^k\cdot a + \tfrac12 + \floor{\tfrac12+a}$.
note: $\tfrac12+a = \floor{\tfrac12+a} + ( (\tfrac12+a)\bmod 1)$ である。
$r\le \tfrac12$ のとき、 $$ \begin{aligned} &\phantom{{}={}} a - ( (2^k+1)\otimes a) \\ &= (-\tfrac12 + \floor{\tfrac12+a} + r) - (2^k\cdot a - \tfrac12 + \floor{\tfrac12+a}) \\ &= -(2^k\cdot a-r) \end{aligned} $$ であり、 $r\ge \tfrac12$ のとき、 $$ \begin{aligned} &\phantom{{}={}} a - ( (2^k+1)\otimes a) \\ &= (-\tfrac12 + \floor{\tfrac12+a} + r) - (2^k\cdot a + \tfrac12 + \floor{\tfrac12+a}) \\ &= -(2^k\cdot a-(r-1) ) \end{aligned} $$ である。$\hfloor a=2^{51-k}$ より、$a\le 2^{51}-\tfrac14$ であることに注意して Case 1-1-2 と同様にして従う。
Case 2: $a\lt 2^{-1022}$
Case 2-1: $a\times (2^k+1)\le 2^{-1022}$
このとき、$a_c = a\times c$ かつ $(a_h, a_l) = (a, 0)$ が成り立つ。 $|a_l| \le \hfloor a\cdot 2^{-(53-k)}$ は明らかであるから、$a \equiv 0 \pmod{\hfloor a\cdot 2^{k-53}}$ を示す。
$\hfloor a=2^{i-1022-52}$ とすると、ある整数 $2^i\le m_a\lt 2^{i+1}$ が存在して $a = m_a\cdot 2^{-1022-52}$ が成り立ち、$m_a\cdot 2^{-1022-52} \equiv 0 \pmod{2^{-1022-52}}$ となる。 $a\times 2^k \lt 2^{-1022}$ より $(i-1022-52)+k\lt -1022$ が成り立つため、 $$ \begin{aligned} &\phantom{{}\implies{}} m_a\cdot 2^{-1022-52} \equiv 0 \pmod{2^{-1022-52}} \\ &\implies m_a\cdot 2^{-1022-52} \equiv 0 \pmod{2^{i-1022-52+k-1-52}} \\ &\iff m_a\cdot 2^{-1022-52} \equiv 0 \pmod{\hfloor a\cdot 2^{k-53}} \end{aligned} $$ が従う。
Case 2-2: $a\times (2^k+1)\gt 2^{-1022}$
$\hfloor{(2^k+1)\cdot a}^{-1}\cdot 2^{52}$ 倍することで、Case 1 に帰着できる。特に、$r'\in\{0, 1\}$ に対して $$ \hfloor{(2^k+1)\cdot a}\cdot 2^{-52}\cdot(\floor{\hfloor{(2^k+1)\cdot a}^{-1}\cdot 2^{52}\cdot a}+r')\in F $$ などが成り立つことに注意せよ。$\qed$
積
さて、こちらも基本的な操作です。FMA が使える前提であれば次のようにできます。
- 入力: $(a, b)\in F\times F$
- 出力: $(p_h, p_l)\in F^2$
- 事前条件:
- $|a\times b|\lt \Theta_{\infty}$, and
- $a\times b \equiv 0 \pmod{2^{-1074}}$
- 事後条件:
- $p_h = a\otimes b$, and
- $p_h + p_l = a \times b$
- 手続き:
- $p_h \gets a\otimes b$ で初期化する。
- $p_l \gets \roundcirc{a\times b + (-p_h)}$ で初期化する。
- $(p_h, p_l)$ を出力する。
Claim 9: 上記の手続きの実行後、事後条件が成り立つ。
Proof
$p_h = a\otimes b$ は明らかであり、$\roundcirc{a\times b + (-p_h)} = a\times b - p_h$ を証明すれば十分。 $ab = 0$ の場合は明らか。$ab\gt 0$ の場合で成り立てば $ab\lt 0$ の場合で成り立つのも明らかなので、$a, b\gt 0$ とする。
Case 1: $ab\ge 2^{-1022}$
$\hfloor a=\hfloor b=2^{52}$ とする。そうでない場合は、それぞれ $\hfloor a^{-1}\cdot 2^{52}$, $\hfloor b^{-1}\cdot 2^{52}$ 倍することで帰着できる。
実数 $|\varepsilon|\le 2^{-53}$ が存在して $a\otimes b = ab+\hfloor{ab}\cdot\varepsilon$ が成り立つ。 $2^{104}\le ab\lt 2^{106}$ であるから、$\hfloor{ab}\le 2^{105}$ である。 $$ \begin{aligned} |(a\otimes b) - (a\times b)| &\le \hfloor{ab}\cdot\varepsilon \\ &\le 2^{105}\cdot 2^{-53} \\ &= 2^{52} \end{aligned} $$ であり、Property 3 から $(a\otimes b)-(a\times b)\in F$ が従う。
Case 2: $ab\lt 2^{-1022}$
このとき、$a\otimes b = a\times b$ が成り立ち、$(p_h, p_l) = (ab, 0)$ となる。$\qed$
事前条件の $a\times b\equiv 0\pmod{2^{-1074}}$ がない場合の反例は、たとえば $(a, b) = (1.5, 2^{-1074})$ です。 $$ \begin{aligned} a\times b &= 1.5\cdot 2^{-1074}, \\ a\otimes b &= 2\cdot 2^{-1074}, \\ (a\otimes b)-(a\times b) &= 0.5\cdot 2^{-1074}\notin F \end{aligned} $$ です。事前条件として $|a\times b|\gt \Theta_0$ を課したとしてもうまくいかないですね。
特に、$|a\times b|\ge 2^{-1022}$ でもだめですね。$(a, b) = ( (2-2^{-52})\cdot 2^{-971}, 2-2^{-52})$ などを考えると、 $$ \begin{aligned} a\times b &= (2-2^{-51}+2^{-105})\cdot 2^{-970}, \\ a\otimes b &= (2-2^{-51})\cdot 2^{-970} \end{aligned} $$ であり、$a\times b \ge 2^{-1022}$ ですが $(a\otimes b) - (a\times b) = 2^{-1075}\notin F$ です。
FMA が使えない場合でも可能です。サブルーチン Round を用いて次のようにできます。
- 入力: $(a, b)\in F\times F$
- 出力: $(p_h, p_l)\in F^2$
- 事前条件:
- $|a\times b| \lt \Theta_{\infty}$,
- $a\times b\equiv 0\pmod{2^{-1074}}$, and
- $\max{\{|a|, |b|}\}\le 2^{1023-27}$
- 事後条件:
- $p_h = a\otimes b$, and
- $p_h + p_l = a \times b$
- 手続き:
- $(u_h, u_l) \gets \textsc{Round}(a, 27)$ で初期化する。
- $(v_h, v_l) \gets \textsc{Round}(b, 27)$ で初期化する。
- $p_h \gets a\otimes b$ で初期化する。
- $p_l \gets ( ( ( (u_h\otimes v_h) \ominus p_h)\oplus (u_h\otimes v_l) ) \oplus (u_l\otimes v_h) ) \oplus (u_l \otimes v_l)$ で初期化する。
- $(p_h, p_l)$ を出力する。
Claim 10: 上記の手続きの実行後、事後条件が成り立つ。
Proof
Claim 9 同様 $a, b\gt 0$ として、$p_h+p_l = a\times b$ を示す。
Case 1: $ab\ge 2^{-1022}$
Claim 9 同様、$\hfloor a=\hfloor b=2^{52}$ とする。Round の事前条件が成り立っていることに注意する。
Round の事後条件より $u_h \equiv v_h \equiv 0 \pmod{2^{27}}$ かつ $|u_l|, |v_l|\le 2^{26}$ が成り立つ。 $$ \begin{aligned} a\times b &= (u_h+u_l) \times (v_h+v_l) \\ &= \underbrace{u_h\times v_h}_{[2^{104}\lldot 2^{106})} + \underbrace{u_h\times v_l}_{[-2^{79}\lldot 2^{79}]} + \underbrace{u_l\times v_h}_{[-2^{79}\lldot 2^{79}]} + \underbrace{u_l\times v_l}_{[-2^{52}\lldot 2^{52}]} \end{aligned} $$ と書ける。ある実数 $|\varepsilon|\le 2^{-53}$ が存在して $|(a\otimes b)-(a\times b)| \le \hfloor{ab}\cdot \varepsilon$ が成り立つから、 $$ \begin{aligned} a\otimes b &= \underbrace{u_h\times v_h\vphantom{\hfloor x}}_{[2^{104}\lldot 2^{106})} + \underbrace{u_h\times v_l\vphantom{\hfloor x}}_{[-2^{79}\lldot 2^{79}]} + \underbrace{u_l\times v_h\vphantom{\hfloor x}}_{[-2^{79}\lldot 2^{79}]} + \underbrace{u_l\times v_l\vphantom{\hfloor x}}_{[-2^{52}\lldot 2^{52}]} + \underbrace{\hfloor{ab}\cdot \varepsilon}_{[-2^{52}\lldot 2^{52}]} \end{aligned} $$ が成り立つ。
$m_u = u_h/2^{27}$ および $m_v = v_h/2^{27}$ とすると $u_h\times v_h = (m_u\times m_v)\times 2^{54}$ となる。 $m_u, m_v\in[2^{25}\lldot 2^{26})\cap\N$ であるから、$u_h \otimes v_h = u_h\times v_h$ が成り立つ。
$u_h\times v_l = (m_u \times v_l)\times 2^{27}$ であり、$|m_u\times v_l|\le 2^{51}$ であるから、$u_h\otimes v_l = u_h\times v_l$ が成り立つ。 同様にして、$u_l\otimes v_h = u_l\times v_h$ である。 また、$|u_l\times v_l| \le 2^{52}$ であるから、$u_l\otimes v_l = u_l\times v_l$ である。
よって、Sterbenz lemma より $$ \begin{aligned} &\phantom{{}={}} (u_h\otimes v_h) \ominus p_h \\ &= (u_hv_h) - (u_hv_h + u_hv_l + u_lv_h + u_lv_l + \hfloor{ab}\cdot\varepsilon) \\ &= -(u_hv_l + u_lv_h + u_lv_l + \hfloor{ab}\cdot\varepsilon) \end{aligned} $$ が成り立つ。
ここで、$a\otimes b = p_h \equiv 0 \pmod{2^{52}}$ かつ $u_hv_h\equiv 0\pmod{2^{54}}$ であるから、$u_hv_h-p_h\equiv 0\pmod{2^{52}}$ が成り立つ。 また、$u_hv_l\equiv 0\pmod{2^{27}}$ であるから、 $$ \begin{aligned} (u_hv_h-p_h)+u_hv_l &= -(u_lv_h + u_lv_l + \hfloor{ab}\cdot\varepsilon) \\ &\equiv 0 \pmod{2^{27}} \end{aligned} $$ かつ $$ \begin{aligned} |(u_hv_h-p_h)+u_hv_l| &\le |u_lv_h| + |u_lv_l| + \hfloor{ab}\cdot|\varepsilon| \\ &\le 2^{79} + 2^{52} + 2^{52} \\ &\lt 2^{80} \end{aligned} $$ であり、ある整数 $0\le m\lt 2^{53}$ に対して $|(u_hv_h-p_h)+u_hv_l| = m\cdot 2^{27}$ と書ける。 よって、 $$ \begin{aligned} ( (u_h\otimes v_h)\ominus p_h)\oplus (u_h\otimes v_l) &= u_hv_h - p_h + u_hv_l \\ &= -(u_lv_h + u_lv_l + \hfloor{ab}\cdot \varepsilon) \end{aligned} $$ が従う。
note: $|u_hv_l|$ が上界に近いとは限らないので、Sterbenz lemma は使えないことに注意せよ。
同様にして、$-(u_lv_h + u_lv_l + \hfloor{ab}\cdot\varepsilon)\equiv 0\pmod{2^{27}}$ かつ $u_lv_h\equiv 0\pmod{2^{27}}$ であるから、 $$ \begin{aligned} ( (u_hv_h-p_h)+u_hv_l)+u_lv_h &= -(u_lv_l + \hfloor{ab}\cdot\varepsilon \\ &\equiv 0 \pmod{2^{27}} \end{aligned} $$ かつ $$ \begin{aligned} |( (u_hv_h-p_h)+u_hv_l)+u_lv_h| &= |u_lv_l| + \hfloor{ab}\cdot|\varepsilon| \\ &\le 2^{52} + 2^{52} \\ &= 2^{53} \end{aligned} $$ であり、$u_lv_l+\hfloor{ab}\cdot\varepsilon\in F$ が従う。 すなわち $$ \begin{aligned} &\phantom{{}={}} ( ( (u_h\otimes v_h) \ominus p_h)\oplus (u_h\otimes v_l) )\oplus (u_l\otimes v_h) \\ &= -(u_lv_h + u_lv_l + \hfloor{ab}\cdot\varepsilon) \oplus u_lv_h \\ &= -(u_lv_l + \hfloor{ab}\cdot\varepsilon) \end{aligned} $$ が従う。
最後に、$\hfloor{ab}\cdot\varepsilon\le 2^{52}$ かつ $\hfloor{ab}\cdot\varepsilon\equiv 1$ より $\hfloor{ab}\cdot\varepsilon\in F$ であるから、 $$ \begin{aligned} &\phantom{{}={}} ( ( ( (u_h\otimes v_h) \ominus p_h)\oplus (u_h\otimes v_l) )\oplus (u_l\otimes v_h) )\oplus (u_l\otimes v_l) \\ &= -(u_lv_l + \hfloor{ab}\cdot\varepsilon) \oplus u_lv_l \\ &= -\hfloor{ab}\cdot\varepsilon \end{aligned} $$ となる。
したがって、$p_h = a\times b + \hfloor{ab}\cdot \varepsilon$ かつ $p_l = -\hfloor{ab}\cdot\varepsilon$ であるから、$p_h+p_l = a\times b$ が成り立つ。
Case 2: $ab\lt 2^{-1022}$
$a\otimes b = a\times b$ に注意する。
$r\ge 0$ に対して $(r_h, r_l) = \textsc{Round}(r, 27)$ とすると $r_h\equiv 0\pmod{\hfloor r\cdot 2^{-25}}$ かつ $|r_l| \le \hfloor r\cdot 2^{-26}$ が成り立つ。 よって、 $$ \begin{aligned} |r_l| &\le \hfloor r\cdot 2^{-26} \\ &\le 2^{-26}\cdot r, \\ r_h &= r - r_l \\ &\le r + r\cdot 2^{-26} \\ &= (1+2^{-26})\cdot r \end{aligned} $$ が成り立つ。よって、 $$ \begin{aligned} u_hv_h &= ab - (u_hv_l + u_lv_h + u_lv_l) \\ &\le ab + |u_hv_l| + |u_lv_h| + |u_lv_l| \\ &\le ab + (1+2^{-26})\cdot a\cdot 2^{-26}\cdot b + 2^{-26}\cdot a\cdot (1+2^{-26})\cdot b + 2^{-52}\cdot ab \\ &= (1 + 2\cdot (1+2^{-26})\cdot 2^{-26} + 2^{-52})\cdot ab \\ &\lt 2ab \\ &\lt 2^{-1021} \end{aligned} $$ であり、$u_hv_h\in F_{-1023}^+\sqcup F_{-1022}^+$ が従う。
同様にして、
- $u_hv_l$,
- $u_lv_h$,
- $u_lv_l$,
- $u_hv_l + u_lv_h + u_lv_l$, and
- $u_lv_h + u_lv_l$
はすべて $F$ の元であることが示せるから、$(p_h, p_l) = (ab, 0)$ が従う。$\qed$
$a\le 2^{1023-27}$ が事前に保証できない場合は、$a\gt 2^{1023-27}$ かどうか判定して true なら $a\otimes 2^{-53}$ などを渡し、結果を $2^{53}$ 倍すればよいです（$b$ についても同様）。
参考文献
- Dekker, T.J. A floating-point technique for extending the available precision. Numer. Math. 18, 224–242 (1971). https://doi.org/10.1007/BF01397083
あとがき
今回は、double-double や triple-double を実装する際の基本操作のみを導入しました。
非正規化数についても触れましたが、「よくわからないので避けたいもの」という印象は拭えても「場合分けが面倒になるもの」という印象が強まってしまっているかもしれません？ わからないので避けたいというのよりはマシな気はしますね。
次回いきなり気が変わって double-double ではないなにかの話を始めるかもしれません？ そのときの気分でテーマが決まるので仕方ないですね。
毎度のことですが、Sterbenz lemma が大活躍すぎます。「$\roundcirc x=x$ を示したくなったらまず Sterbenz lemma を使えないか疑え」まである気がします。 もちろん、mod の性質など、それ以外のものを使う局面もしばしばありますが、Sterbenz lemma がお手軽どころの頻出パターンという感じですね。
なんか疲れてしまいました。証明の方針が全然思いつかなかったり、場合分けをがんばる方法で書き終えたと思ったらまとめられることに気づいたり、なにも考えられなくなったり、シンプルに嘘を書いていたり、ひ〜という感じです。 お気持ちパートも流れるように書けなくなってきました。不調？
さっさと元気になって続編を勉強したいですね。correct rounding に興味のあるフォロヮがどの程度いるかは知らないですが、えびちゃんは知りたがっているので。
おわり
おわりです。
