ABC 405 E についての雑談

apple, banana, cherry, durian とかにしてほしかった...

atcoder.jp

考察

コンテスト中にやった考察に沿って話を進めます。

立式

一番右の A と一番右の B に注目します。それらの間に C が含まれるケースは次のようになります。

AAAABBBABBA ++ A ++ BCCBB ++ B ++ CDCCCCDCDD

すなわち、次を連結したものになります。

• $A-1$ 個の A と、$i$ ($0\le i\lt B$) 個の B を任意の順に並べたもの
• $1$ 個の A
• $B-i-1$ 個の B と、$j$ ($1\le j\le C$) 個の C を任意の順に並べたもの
• $1$ 個の B
• $C-j$ 個の C と $D$ 個の D を任意の順に並べたもの

よって、この通り数は次の式で表せます。 $$ \begin{aligned} &\phantom{{}={}} \sum_{i=0}^{B-1} \sum_{j=1}^C \binom{A-1+i}{i} \binom{B-i-1+j}{j} \binom{C-j+D}{D} \\ &= \sum_{i=0}^{B-1} \binom{A-1+i}{i} \sum_{j=1}^C \binom{B-i-1+j}{j} \binom{C-j+D}{D} \end{aligned} $$

また、それらの間に C が含まれないケースは次のようになります。

AABABBAABAABABBA ++ CDDCDCDCCCCDCC

すなわち、次を連結したものになります。

• $A$ 個の A と、$B$ 個の B を任意の順に並べたもの
• $C$ 個の C と、$D$ 個の D を任意の順に並べたもの

特に、A と C の制約から、一番右の B よりも一番右の A の方が右に来るケースもこちらに含まれます。 この通り数は次のようになります。

$$ \binom{A+B}{B} \binom{C+D}{D} $$

計算量改善

後者のケースはいいとして、前者のケースの内側の $\sum$ を高速に計算する必要があります。 $$ \sum_{i=0}^{B-1} \binom{A-1+i}{i} \underbrace{\sum_{j=1}^C \binom{B-i-1+j}{j} \binom{C-j+D}{D}}_{b_i} $$ すなわち、この $b_i$ の部分です。とりあえず Wolfram|Alpha に投げます。 $$ b_i = \frac{C+D}C\binom{C+D-1}{D} \cdot({}_2 F_1(-C, B-i; -C-D; 1)-1) $$ 二項係数と階乗の部分を整理すると下記になります。 $$ b_i = \frac{(C+D)!}{D!\,C!} \cdot({}_2 F_1(-C, B-i; -C-D; 1)-1) $$

ここで、${}_2 F_1(a, b; c; z)$ は 超幾何関数 (hypergeometric function) と呼ばれ、下記で与えられるらしいです*1。 $$ {}_2 F_1(a, b; c; z) = \sum_{n=0}^{\infty} \frac{(a)_n\,(b)_n}{(c)_n}\cdot\frac{z^n}{n!} $$ $(a)_n$ は Pochhammer の上昇階乗冪で、$(a)_n = a\cdot(a+1)\cdot\cdots\cdot(a+n-1)$ です。Wolfram|Alpha では Hypergeometric2F1 という名前で定義されています。

困っているようす

$c\lt a\lt 0\le b$ として、 $$ \begin{aligned} {}_2 F_1(a, b; c; 1) &= \sum_{n=0}^{\infty} \frac{(a)_n\,(b)_n}{(c)_n}\cdot\frac{1}{n!} \\ &= \sum_{n=0}^{|a|} \frac{(a)_n\,(b)_n}{(c)_n}\cdot\frac{1}{n!} \\ &= \sum_{n=0}^{|a|} \frac{(a)_n}{(c)_n}\cdot\frac{(b)_n}{n!} \\ &= \sum_{n=0}^{|a|} \frac{a\cdot(a+1)\cdot\cdots\cdot(a+n-1)}{c\cdot(c+1)\cdot\cdots\cdot(c+n-1)}\cdot\frac{(b)_n}{n!} \\ &= \sum_{n=0}^{|a|} \frac{|a|\cdot(|a|-1)\cdot\cdots\cdot(|a|-(n-1) )}{|c|\cdot(|c|+1)\cdot\cdots\cdot(|c|-(n-1) )}\cdot\frac{(b)_n}{n!} \\ &= \sum_{n=0}^{|a|} \frac{|a|!}{(|a|-n)!} \frac{(b+n-1)!}{(b-1)!} \frac{(|c|-n)!}{|c|!} \frac1{n!} \\ &= \sum_{n=0}^{|a|} \frac{|a|!}{(|a|-n)!\,n!} \frac{(b+n-1)!}{(b-1)!\,n!} \frac{(|c|-n)!\,n!}{|c|!} \\ &= \sum_{n=0}^{|a|} \binom{|a|}{n} \binom{b+n-1}{n} \binom{|c|}{n}^{-1} \end{aligned} $$ ここからどうすれば...？ リンク先の式 (74) によれば $$ {}_2 F_1(a, b; c; 1) = \frac{\Gamma(c)\,\Gamma(c-a-b)}{\Gamma(c-a)\,\Gamma(c-b)} $$ らしいです。ここでは $c$ が負なのでハチャメチャに発散しそうです。困った。

$a\lt 0$ や $\Gamma(n) = (n-1)\,\Gamma(n-1)$ などから $$ \begin{aligned} {}_2 F_1(a, b; c; 1) &= \frac{\Gamma(c)\,\Gamma(c-a-b)}{\Gamma(c-a)\,\Gamma(c-b)} \\ &= \frac{\Gamma(c)}{\Gamma(c-a)}\frac{\Gamma(c-a-b)}{\Gamma(c-b)} \\ &= \frac{\Gamma(c)}{(c-a-1)\,\Gamma(c-a-1)}\frac{(c-a-b-1)\,\Gamma(c-a-b-1)}{\Gamma(c-b)} \\ &= \cdots \\ &= \frac{(c-a-b-1)\cdot\cdots\cdot(c-b)}{(c-a-1)\cdot\cdots\cdot c} \end{aligned} $$ としてよいのであれば、計算はできそうです。今回計算したいのは ${}_2 F_1(-C, B-i; -C-D; 1)$ だったので、 $$ \begin{aligned} &\phantom{{}={}} {}_2 F_1(-C, B-i; -C-D; 1) \\ &= \frac{-C-D-(-C)-(B-i)-1}{-C-D-(-C)-1}\cdot\cdots\cdot\frac{-C-D-(B-i)}{-C-D} \\ &= \frac{-D-B+i-1}{-D-1}\cdot\cdots\cdot\frac{-C-D-B+i}{-C-D} \\ &= \frac{(B-i)+(D+1)}{D+1}\cdot\cdots\cdot\frac{(B-i)+(D+C)}{D+C} \\ &= \frac{( (B-i)+(D+C) )!}{( (B-i)+D)!}\cdot\frac{D!}{(D+C)!} \end{aligned} $$ となりそうです。

たとえば、$(B, C, D, i) = (10, 5, 2, 3)$ としてみると、${}_2 F_1(-5, 7; -7; 1) = \frac{286}3$ となり*2、 $$ \begin{aligned} b_i &= \frac{7!}{2!\,5!}\cdot({}_2 F_1(-5, 7; -7; 1)-1) \\ &= 21\cdot \left(\frac{286}3-1\right) \\ &= 7\cdot(286-3) \\ &= 1981 \end{aligned} $$ となり、元々の $\sum$ で計算した値と一致します。$\eod$

コンテスト本番では、超幾何関数はどうしようもないと思っていたので別の方針を考えました。 $B$, $C$, $D$ を $1$ ずつ変えながら値を眺めることで、下記に気づきます。

Observation 1: $b_B = 0$ かつ、各 $0\le i\lt B$ に対して下記が成り立つ。 $$ b_i = b_{i+1} + \binom{(B-1)-i+C+D}{C-1} $$

$(B-1)-i$ の部分や漸化式の向きなどから、$b'_i$ を下記で定義します。 $$ b'_i = \sum_{j=1}^C \binom{i+j}{i} \binom{C-j+D}{D} $$ なお、$b'_i = b_{(B-1)-i}$ です。

Conjecture 2: $b_{-1} = 0$ かつ、各 $0\le i\le B-1$ に対して下記が成り立つ。 $$ b'_i = b'_{i-1} + \binom{i+C+D}{C-1} $$

More observations

$i=0$ のとき、Hockey-stick identity から $$ \begin{aligned} b'_0 &= \sum_{j=1}^C \binom jj\binom{C-j+D}D \\ &= \sum_{j=1}^C \binom{C-j+D}D \\ &= \sum_{j'=0}^{C-1} \binom{D+j'}D \\ &= \binom{D+C}{D+1} \end{aligned} $$ が成り立つ。一方、 $$ \begin{aligned} b'_0 &= b'_{-1} + \binom{C+D}{C-1} \\ &= 0 + \binom{D+C}{D+1} \end{aligned} $$ より、$i = 0$ のときは成り立つ。

$$ \begin{aligned} b'_i - b'_{i-1} &= \sum_{j=1}^C \binom{i+j}j\binom{C-j+D}D - \sum_{j=1}^C \binom{i-1+j}j\binom{C-j+D}D \\ &= \sum_{j=1}^C \left( \binom{i+j}j - \binom{i-1+j}j\right) \binom{C-j+D}D \\ % &= \sum_{j=1}^C \left( \frac{(i+j)!}{i!\, j!} - \frac{(i+j-1)!}{(i-1)!\, j!}\right) \binom{C-j+D}D \\ % &= \sum_{j=1}^C \left( \frac{(i+j)!}{i!\, j!} - \frac{i}{i+j}\frac{(i+j)!}{i!\, j!}\right) \binom{C-j+D}D \\ % &= \sum_{j=1}^C \frac{j}{i+j}\frac{(i+j)!}{i!\, j!} \binom{C-j+D}D \\ % &= \sum_{j=1}^C \frac{(i+j-1)!}{i!\, (j-1)!} \binom{C-j+D}D \\ &= \sum_{j=1}^C \binom{i-1+j}{j-1} \binom{C-j+D}D \\ &= \sum_{j=0}^{C-1} \binom{i+j}{j} \binom{C-1-j+D}D \\ &= \sum_{j=0}^{C-1} \binom{i+j}{i} \binom{C-1-j+D}D \\ % &= \sum_{j=0}^{C-1} \binom{i+j}{i} \binom{C-1-j+D}D \\ % &= \binom{i+C+D}{i+D+1} \\ % &= \binom{i+C+D}{C-1}. \quad\qed \end{aligned} $$

よって、下記が成り立てばよい。 $$ \sum_{j=0}^{C-1} \binom{i+j}{i} \binom{C-1-j+D}D = \binom{i+D+C}{i+D+1}. $$

Lemma 3: 下記が成り立つ。

$$ \sum_{i=0}^{k-1} \binom{n+i}{n} \binom{m+k-1-i}{m} = \binom{n+m+k}{n+m+1}. $$

note: Pascal の三角形の $n$ 列目と $m$ 列目から $k$ 個持ってきたものの畳み込みに相当する。

                           .  
                         .   .  
                       .   .   A  
                     .   .   B   .  
                   .   .   C   .   .  
                 .   .   D   .   .   .  
               .   .   E   .   .   .   e  
             .   .   .   .   .   .   d   .  
           .   .   .   .   .   .   c   .   .  
         .   .   .   .   .   .   b   .   .   .  
       .   .   .   .   .   .   a   .   .   .   .  
     .   .   .   .   .   .   .   .   .   .   .   .  
   .   .   .   .   .   .   .   .   .   .   .   .   .  
 .   .   .   .   .   .   .   .   .  [*]  .   .   .   .  

Proof

$m$, $k$ に関する帰納法で示す。

$$ P(m, k) \iff \ForallL{n}{\sum_{i=0}^{k-1} \binom{n+i}{n} \binom{m+k-1-i}{m} = \binom{n+m+k}{n+m+1}} $$ とする。

To-be-proved 1: $P(m, 1)$

$$ \begin{aligned} \sum_{i=0}^{1-1} \binom{n+i}{n} \binom{m+1-1-i}{m} &= \binom{n}{n} \binom{m}{m} \\ &= 1, \\ \binom{n+m+1}{n+m+1} &= 1. \end{aligned} $$

To-be-proved 2: $P(m, k+1) \wedge P(m+1, k) \implies P(m+1, k+1)$

$$ \begin{aligned} &\phantom{{}={}} \sum_{i=0}^{(k+1)-1} \binom{n+i}{n} \binom{(m+1)+(k+1)-1-i}{m+1} \\ &= \sum_{i=0}^{(k+1)-1} \binom{n+i}{n} \left(\binom{(m+1)+(k+1)-1-i-1}{(m+1)-1} \right. \\ & \qquad\qquad\qquad\qquad\quad {} + \left. \binom{(m+1)+(k+1)-1-i-1}{m+1}\right) \\ &= \sum_{i=0}^{(k+1)-1} \binom{n+i}{n} \binom{(m+1)+(k+1)-1-i-1}{(m+1)-1} \\ &\qquad\qquad {} + \sum_{i=0}^{k-1} \binom{n+i}{n} \binom{(m+1)+(k+1)-1-i-1}{m+1} \\ &\qquad\qquad {} + \binom{n+k}{n} \binom{(m+1)+(k+1)-1-k-1}{m+1} \\ &= \sum_{i=0}^{(k+1)-1} \binom{n+i}{n} \binom{m+(k+1)-1-i}{m} + \sum_{i=0}^{k-1} \binom{n+i}{n} \binom{(m+1)+k-1-i}{m+1} + 0 \\ &= \binom{n+m+(k+1)}{n+m+1} + \binom{n+(m+1)+k}{n+(m+1)+1} \\ &= \binom{n+(m+1)+(k+1)}{n+(m+1)+1}. \quad\qed \end{aligned} $$

なんかふつうにお出しされた。最初から教えてよ〜 www.wolframalpha.com

変数を $k$ から $d$ に変えたら出してくれなくなり、これマジ？となりました。 www.wolframalpha.com

比較のために見ておくと、Vandermonde identity は下記のような畳み込みです。うまいこと帰着できたりするのかな？

                               .  
                             .   .  
                           .   .   .  
                         .   .   .   .  
                       A   B   C   D   .  
                     .   .   .   .   .   .  
                   .   .   .   .   .   .   .  
                 d   c   b   a   .   .   .   .  
               .   .   .   .   .   .   .   .   .  
             .   .   .   .   .   .   .   .   .   .  
           .   .   .   .   .   .   .   .   .   .   .  
         .   .   .  [*]  .   .   .   .   .   .   .   .  

Corollary 4: 任意の $n\ge 1$ に対して下記が成り立つ。

$$ \sum_{i=0}^{k-1} \binom{n+i}{n} \binom{m+k-1-i}{m} = \sum_{i=0}^{k-1} \binom{n-1+i}{n-1} \binom{m+k-i}{m+1} $$

Proof: Lemma 3 より明らか。

これを示して $n = 0$ に帰着させることで Hockey-stick identity から Lemma 3 を示す方針を考えていましたが、こちらが Corollary になりました。

おまけ

ここまでの考察で、$b_i$ を incremental に求めることで $i$ ごとに定数時間で求められることがわかりました。コンテスト中はそれで十分でしたが、もう少し進めてみます。

下記が成り立ちます。 $$ \begin{aligned} b_i &= \sum_{j=1}^C \binom{B-i-1+j}{j} \binom{C-j+D}{D} \\ &= \sum_{j=1}^C \binom{B-i+j-1}{B-i-1} \binom{C-j+D}{D} \\ &= \sum_{j=0}^{C-1} \binom{B-i+j}{B-i-1} \binom{C-1-j+D}{D} \\ &= \left( \sum_{j=0}^{C} \binom{B-i-1+j}{B-i-1} \binom{C-j+D}{D} - \binom{B-i-1}{B-i-1} \binom{C+D}{D} \right) \\ &= \binom{(B-i-1)+D+(C+1)}{(B-i-1)+D+1} - \binom{C+D}{D} \\ &= \binom{B+C+D-i}{B+D-i} - \binom{C+D}{D} \\ &= \binom{B+C+D-i}{C} - \binom{C+D}{D} \\ \end{aligned} $$

よって、冒頭に挙げた前者のケースの通り数は $$ \begin{aligned} &\phantom{{}={}} \sum_{i=0}^{B-1} \binom{A-1+i}{i} \sum_{j=1}^C \binom{B-i-1+j}{j} \binom{C-j+D}{D} \\ &= \sum_{i=0}^{B-1} \binom{A-1+i}{i} \left( \binom{B+C+D-i}{C} - \binom{C+D}{D} \right) \\ &= \sum_{i=0}^{B-1} \binom{A-1+i}{A-1} \binom{B+C+D-i}{C} - \binom{C+D}{D} \sum_{i=0}^{B-1} \binom{A-1+i}{A-1} \\ &= \sum_{i=0}^{B-1} \binom{A-1+i}{A-1} \binom{B+C+D-i}{C} - \binom{A-1+B}{A} \binom{C+D}{D} \end{aligned} $$ となります。後者のケースと合わせて、全体の答えは $$ \begin{aligned} &\phantom{{}={}} \sum_{i=0}^{B-1} \binom{A-1+i}{A-1} \binom{B+C+D-i}{C} \\ &\qquad\quad {} - \binom{A-1+B}{A} \binom{C+D}{D} + \binom{A+B}{A} \binom{C+D}{D} \\ &= \sum_{i=0}^{B-1} \binom{A-1+i}{A-1} \binom{B+C+D-i}{C} + \binom{A+B-1}{A-1} \binom{C+D}{D} \\ &= \sum_{i=0}^{B} \binom{A-1+i}{A-1} \binom{B+C+D-i}{C} \end{aligned} $$ とできますね。

また、対称性から、入力が $(A, B, C, D)$ のときの答えと $(D, C, B, A)$ のときの答えが一致することもわかります。実際、 $$ \begin{aligned} &\phantom{{}={}} \sum_{i=0}^{B} \binom{A-1+i}{A-1} \binom{B+C+D-i}{C} \\ &= \sum_{i=0}^{C} \binom{D-1+i}{D-1} \binom{C+B+A-i}{B} \\ &= \sum_{i=0}^{C} \binom{D-1+C-i}{D-1} \binom{B+A+i}{B} \\ &= \sum_{i=0}^{C} \binom{A+B+i}{B} \binom{D-1+C-i}{D-1} \\ \end{aligned} $$ となり、公式解説 と同じ式が得られました。

ということで、この畳み込みも高速化できないものかね？という気持ちになります。下記が成り立つらしいです。 $$ \begin{aligned} &\phantom{{}={}} \sum_{i=0}^{k-1} \binom{n+r+i}{n} \binom{m+k-1-i}{m} \\ &= \binom{k+m-1}{m}\binom{n+r}{r} {}_3 F_2 (1, 1-k, n-r+1; -k-m+1, r+1; 1). \end{aligned} $$ なんやかんやで答えは下記になります。 $$ \begin{aligned} \binom{A+B}{A} \binom{C+D-1}{C} {}_3 F_2(1, -C, A+B+1; A+1, -C-D+1; 1) \end{aligned} $$ ${}_3 F_2$ を高速に有理数表示で計算する方法（の存在）がわからず、おわりです...

定義通りの $$ {}_p F_q(a_1, \dots, a_p; b_1, \dots, b_q; z) = \sum_{n=0}^{\infty} \frac{(a_1)_n\cdot\cdots\cdot(a_p)_n}{(b_1)_n\cdot\cdots\cdot(b_q)_n}\frac{z^n}{n!} $$ の $n$ 項目の各 factor を持っておけば、次の項は $O(p+q)$ 時間程度で計算できます。$n\ge C$ のとき $(-C)_n = 0$ なので先頭 $C$ 項だけ見ればよいことから $O(C)$ 時間で求められはしますが、なんだかな〜という感じですね。

もしかして二項係数関連のどうにもならない積和を超幾何関数として追いやっているだけですか？

おまけのおまけ

$$ \begin{aligned} {}_3 F_2 (1, 1-k, n+1; -k-m+1, 1; 1) &= {}_2 F_1 (1-k, n+1; -k-m+1; 1) \end{aligned} $$ が成り立つらしいので、$r=0$ を考えて $$ \begin{aligned} {}_2 F_1 (1-k, n+1; -k-m+1; 1) = \binom{n+m+k}{n+m+1} \binom{k+m-1}{m}^{-1} \end{aligned} $$ であり、$(k-1, n+1, m-a) = (a, b, c)$ として $$ \begin{aligned} {}_2 F_1 (-a, b; -c; 1) = \binom{b+c}{b+c-a} \binom{c}{c-a}^{-1} \end{aligned} $$ を得ます。

あとがき

たぶん失敗方針ではありそうな気もしつつ、おもしろい等式を知ることができたのでよかったな〜という感じです。 超幾何関数とも知り合いになれました。

この手の問題で、「とりあえず愚直を書いて、差分を睨んでいたら二項係数が出てくるだろうから、エスパーして合わせればなんとかなる」みたいなやり方ばかりやっているのはやっぱりよくないと思います。でも頭を使わずに済むから楽なんですよね。頭を使いたくないなら ABC をやめればいいのに... とも思いますし、ABC で頭を使っている時点で負けという気もします。

気づいたら日曜が終わろうとしているのは納得いきません。本題の浮動小数点数シリーズの方は進捗がありません。えびちゃんだって進捗を出したいと思っています。

おわり

おわりです。

*1:$,$ と $;$ は誤植ではなく、一般には ${}_p F_q(a_1, \dots, a_p; b_1, \dots, b_q; z)$ という形をしているみたいです。

*2:実際、Wolfram|Alpha に計算させてもその値になりました。また、整数からちょっとだけずらして $\Gamma$ を計算すると、たしかにそれっぽい値に近づきそうではありました。

correct rounding への道 (5) table-maker’s dilemma

やりましょう。

前回、double-double や triple-double の話をしましたが、「triple-double は必要なのか？」とか「quadruple-double は不要なのか？」とかについて触れた方がいいかなと思ったので、先にそういう話をします。重要なお話です。

disclaimer: 重要なお話ではあるものの、計算するにあたってのアルゴリズムなどの話ではないです。

記法

多くの記法については前回定義したものを使います。特に、double で表せる浮動小数点数全体からなる集合を $F$ とします。

今回以降、実数の表記として主に指数表記つきの 16 進法や 2 進法を使っていきます。たとえば $$ \pi = {\small 3.14159265358979}{\footnotesize 323846264338327}{\scriptsize 950288419716939}{\tiny 937510582097494{\ldots}} $$ ではなく $$ \begin{aligned} \pi &= \texttt{1}.\texttt{{\small 921FB54442D18}{\footnotesize 469898CC51701}{\scriptsize B839A252049C1}{\tiny 114CF98E80417{\ldots}}}_{(16)}\times 2^{1} \\ &= 1. \underbrace{\footnotesize 1001}_{\texttt{9}} {\footnotesize\,} \underbrace{\footnotesize 0010}_{\texttt{2}} {\footnotesize\,} \underbrace{\footnotesize 0001}_{\texttt{1}} {\footnotesize\,} \underbrace{\footnotesize 1111}_{\texttt{F}} {\footnotesize\,} \underbrace{\footnotesize 1011}_{\texttt{B}} {\footnotesize\,} \underbrace{\footnotesize 0101}_{\texttt{5}} {\footnotesize\,} \underbrace{\footnotesize 0100}_{\texttt{4}} {\footnotesize\,} \underbrace{\footnotesize 0100}_{\texttt{4}} {\footnotesize\,} \underbrace{\footnotesize 0100}_{\texttt{4}} {\footnotesize\,} \underbrace{\footnotesize 0010}_{\texttt{2}} {\footnotesize\,} \underbrace{\footnotesize 1101}_{\texttt{D}} {\footnotesize\,} \underbrace{\footnotesize 0001}_{\texttt{1}} {\footnotesize\,} \underbrace{\footnotesize 1000}_{\texttt{8}} {\footnotesize\,} \\ & \qquad \underbrace{\footnotesize 0100}_{\texttt{4}} {\footnotesize\,} \underbrace{\footnotesize 0110}_{\texttt{6}} {\footnotesize\,} \underbrace{\footnotesize 1001}_{\texttt{9}} {\footnotesize\,} \underbrace{\footnotesize 1000}_{\texttt{8}} {\footnotesize\,} \underbrace{\footnotesize 1001}_{\texttt{9}} {\footnotesize\,} \underbrace{\footnotesize 1000}_{\texttt{8}} {\footnotesize\,} \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 0101}_{\texttt{5}} {\footnotesize\,} \underbrace{\footnotesize 0001}_{\texttt{1}} {\footnotesize\,} \underbrace{\footnotesize 0111}_{\texttt{7}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0001}_{\texttt{1}} {\footnotesize\,} \\ & \qquad \underbrace{\footnotesize 1011}_{\texttt{B}} {\footnotesize\,} \underbrace{\footnotesize 1000}_{\texttt{8}} {\footnotesize\,} \underbrace{\footnotesize 0011}_{\texttt{3}} {\footnotesize\,} \underbrace{\footnotesize 1001}_{\texttt{9}} {\footnotesize\,} \underbrace{\footnotesize 1010}_{\texttt{A}} {\footnotesize\,} \underbrace{\footnotesize 0010}_{\texttt{2}} {\footnotesize\,} \underbrace{\footnotesize 0101}_{\texttt{5}} {\footnotesize\,} \underbrace{\footnotesize 0010}_{\texttt{2}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0100}_{\texttt{4}} {\footnotesize\,} \underbrace{\footnotesize 1001}_{\texttt{9}} {\footnotesize\,} \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 0001}_{\texttt{1}} {\footnotesize\,} \\ & \qquad \underbrace{\footnotesize 0001}_{\texttt{1}} {\footnotesize\,} \underbrace{\footnotesize 0001}_{\texttt{1}} {\footnotesize\,} \underbrace{\footnotesize 0100}_{\texttt{4}} {\footnotesize\,} \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 1111}_{\texttt{F}} {\footnotesize\,} \underbrace{\footnotesize 1001}_{\texttt{9}} {\footnotesize\,} \underbrace{\footnotesize 1000}_{\texttt{8}} {\footnotesize\,} \underbrace{\footnotesize 1110}_{\texttt{E}} {\footnotesize\,} \underbrace{\footnotesize 1000}_{\texttt{8}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0100}_{\texttt{4}} {\footnotesize\,} \underbrace{\footnotesize 0001}_{\texttt{1}} {\footnotesize\,} \underbrace{\footnotesize 0111}_{\texttt{7}} {\footnotesize\,} {\footnotesize \,\dots\,}_{(2)} \times 2^1 \end{aligned} $$ のように書きます。指数の底は $2$ です。2 進法で書く場合は 4 桁ごとに区切って 16 進法での値を添えるようにします。

言語によっては、浮動小数点型の 16 進法リテラルがサポートされていて、0x1.921FB54442D18p+1 のように書くこともできます。

C ではリテラルとして可能で、Python でも float.fromhex() を使うことで可能です。Rust では標準ではできず、hexfloat2 などのクレートを使う必要があります*1。また、出力に関しても、C では printf("%a", x) としたり、Python では x.hex() とすることで可能です。Rust はお察しです。

å°Žå…¥

例

$\sin(x)$ を計算したくなったことにしましょう。

ということで、次の手続き $\textsc{Fr54-Sin}$ を考えます*2。

• 任意の $x\in F\smallsetminus\{-\infty, 0, +\infty\}$ を入力とする。
• 下記を満たす $y\in \Q$ を出力する。
  • ある組 $(s, e, m) \in \{0, 1\}\times\Z\times([2^{53}\lldot 2^{54})\cap\N)$ が一意に存在して ${y = (-1)^s\cdot m\cdot 2^{e-53}}$ が成り立つ。
  • 上記の $e$ に対して $|y-\sin(x)|\le 2^{e-53}$ が成り立つ。

$\roundcirc{\sin(x)} = \roundcirc{\textsc{Fr54-Sin}(x)}$ として求めようというわけです。

たとえば $x = \texttt{1}.\texttt{D037CB27EE6DF}_{(16)}\times 2^{-3}\in F$ を考えます。 実際の計算結果は下記の通りです。 $$ \begin{aligned} &\phantom{{}={}} \textsc{Fr54-Sin}(x) \\ &= \texttt{1}.\texttt{CC40C3805229A8}_{(16)}\times 2^{-3} \\ &= 1. \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 0100}_{\texttt{4}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 0011}_{\texttt{3}} {\footnotesize\,} \underbrace{\footnotesize 1000}_{\texttt{8}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0101}_{\texttt{5}} {\footnotesize\,} \underbrace{\footnotesize 0010}_{\texttt{2}} {\footnotesize\,} \underbrace{\footnotesize 0010}_{\texttt{2}} {\footnotesize\,} \underbrace{\footnotesize 1001}_{\texttt{9}} {\footnotesize\,} \underbrace{\footnotesize 1010}_{\texttt{A}} {\footnotesize\,} \underbrace{\footnotesize 1\hphantom{000}}_{\texttt{8}} {}_{(2)} \times 2^{-3} \end{aligned} $$ 最後の部分について、$\underbrace{1010}_{\texttt{A}}\,100\ldots01\ldots$ と $\underbrace{1010}_{\texttt{A}}\,011\ldots10\ldots$ のどちらを丸めた結果なのかを判別できません。前者であれば、$\roundcirc{\sin(x)} = \texttt{1}.\texttt{CC40C3805229B}_{(16)}\times 2^{-3}$ ですし、後者であれば $\roundcirc{\sin(x)} = \texttt{1}.\texttt{CC40C3805229A}_{(16)}\times 2^{-3}$ なので、$\textsc{Fr54-Sin}(x)$ を用いて求めることはできません。

note: $x\in\Q\smallsetminus\{0\}$ より $\sin(x)\notin\Q$ なので、$\sin(x) = \texttt{1}.\texttt{CC40C3805229A8}_{(16)}\times 2^{-3}$ となることはありません。

ということで、もう少し高い精度で計算してみましょう。同様に $\textsc{Fr55-Sin}(x)$ を定義して計算します。 $$ \begin{aligned} &\phantom{{}={}} \textsc{Fr55-Sin}(x) \\ &= \texttt{1}.\texttt{CC40C3805229A8}_{(16)}\times 2^{-3} \\ &= 1. \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 0100}_{\texttt{4}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 0011}_{\texttt{3}} {\footnotesize\,} \underbrace{\footnotesize 1000}_{\texttt{8}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0101}_{\texttt{5}} {\footnotesize\,} \underbrace{\footnotesize 0010}_{\texttt{2}} {\footnotesize\,} \underbrace{\footnotesize 0010}_{\texttt{2}} {\footnotesize\,} \underbrace{\footnotesize 1001}_{\texttt{9}} {\footnotesize\,} \underbrace{\footnotesize 1010}_{\texttt{A}} {\footnotesize\,} \underbrace{\footnotesize 10\hphantom{00}}_{\texttt{8}} {}_{(2)} \times 2^{-3}. \end{aligned} $$ まだ同様の理由でだめですね。さらに高い精度でやりましょう。

〜しばらく繰り返し〜

$$ \begin{aligned} &\phantom{{}={}} \textsc{Fr105-Sin}(x) \\ &= \texttt{1}.\texttt{CC40C3805229A8000000000000}_{(16)}\times 2^{-3} \\ &= 1. \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 0100}_{\texttt{4}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 0011}_{\texttt{3}} {\footnotesize\,} \underbrace{\footnotesize 1000}_{\texttt{8}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0101}_{\texttt{5}} {\footnotesize\,} \underbrace{\footnotesize 0010}_{\texttt{2}} {\footnotesize\,} \underbrace{\footnotesize 0010}_{\texttt{2}} {\footnotesize\,} \underbrace{\footnotesize 1001}_{\texttt{9}} {\footnotesize\,} \underbrace{\footnotesize 1010}_{\texttt{A}} {\footnotesize\,} \underbrace{\footnotesize 1000}_{\texttt{8}} \\ & \qquad \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} {}_{(2)} \times 2^{-3}. \end{aligned} $$ ま、まだだめです...。

$$ \begin{aligned} &\phantom{{}={}} \textsc{Fr106-Sin}(x) \\ &= \texttt{1}.\texttt{CC40C3805229A7FFFFFFFFFFFF8}_{(16)}\times 2^{-3} \\ &= 1. \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 0100}_{\texttt{4}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 0011}_{\texttt{3}} {\footnotesize\,} \underbrace{\footnotesize 1000}_{\texttt{8}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0101}_{\texttt{5}} {\footnotesize\,} \underbrace{\footnotesize 0010}_{\texttt{2}} {\footnotesize\,} \underbrace{\footnotesize 0010}_{\texttt{2}} {\footnotesize\,} \underbrace{\footnotesize 1001}_{\texttt{9}} {\footnotesize\,} \underbrace{\footnotesize 1010}_{\texttt{A}} {\footnotesize\,} \underbrace{\footnotesize 0111}_{\texttt{7}} {\footnotesize\,} \\ & \qquad \underbrace{\footnotesize 1111}_{\texttt{F}} {\footnotesize\,} \underbrace{\footnotesize 1111}_{\texttt{F}} {\footnotesize\,} \underbrace{\footnotesize 1111}_{\texttt{F}} {\footnotesize\,} \underbrace{\footnotesize 1111}_{\texttt{F}} {\footnotesize\,} \underbrace{\footnotesize 1111}_{\texttt{F}} {\footnotesize\,} \underbrace{\footnotesize 1111}_{\texttt{F}} {\footnotesize\,} \underbrace{\footnotesize 1111}_{\texttt{F}} {\footnotesize\,} \underbrace{\footnotesize 1111}_{\texttt{F}} {\footnotesize\,} \underbrace{\footnotesize 1111}_{\texttt{F}} {\footnotesize\,} \underbrace{\footnotesize 1111}_{\texttt{F}} {\footnotesize\,} \underbrace{\footnotesize 1111}_{\texttt{F}} {\footnotesize\,} \underbrace{\footnotesize 1111}_{\texttt{F}} {\footnotesize\,} \underbrace{\footnotesize 1\hphantom{000}}_{\texttt{8}} {\footnotesize\,} {}_{(2)} \times 2^{-3}. \end{aligned} $$ わ、ようやく区別がつきました。 これにより $$ \sin(x) \lt \texttt{1}.\texttt{CC40C3805229A8}_{(16)}\times 2^{-3} $$ が言えるので、 $$ \roundcirc{\sin(x)} = \texttt{1}.\texttt{CC40C3805229A}_{(16)}\times 2^{-3} $$ であるとわかったことになります。めでたしめでたし*3。実際のところ、 $$ \sin(x) = \texttt{1}.\texttt{{\small CC40C3805229A}{\footnotesize 7FFFFFFFFFFFF}{\scriptsize 83E76BBF0FCBE}{\tiny 19AFB5AC38007{\ldots}}} {\,}_{(16)} \times 2^{-3} $$ です。

用語について

さて、上記の例では 106-bit での計算で十分でしたが、これが「106-bit で十分ですよ」ということを予め知っておくことは可能でしょうか？ 平方根など一部の例では可能ですが、三角関数を始めとする多くの関数ではこれは難しいということになっています。

William M. Kahan*4は次のように記しています^[6]。

Why can’t Y^W be rounded within half an ulp like SQRT ? Because nobody knows how much computation it would cost to resolve what I long ago christened “The Table-Maker’s Dilemma”.

拙訳：なぜ $Y^W$ は、平方根のように 0.5 ULP 以内（の誤差で収まるよう）に丸められないのか？ それは、以前私が 数表作成者のジレンマ (Table-Maker’s Dilemma) と名づけた事象を解決するためにどれだけの計算が必要かを誰も知らないからだ。

note: table-maker’s dilemma は、しばしば TMD と略されます。

この後、$\tfrac18\cosh(\pi\sqrt{163}) - (2^{53}-1)$ を 53 ビットに正確に丸めた値を計算しようとして、精度を上げつつ ${\small 7401389035307055}.{\footnotesize 49999978}$ や ${\small 7401389035307055}.{\footnotesize 5000000000015}$ を得る話が続きます。先ほどの $\textsc{Fr106-Sin}(x)$ の例と同じような感じです。

No general way exists to predict how many extra digits will have to be carried to compute a transcendental expression and round it correctly to some preassigned number of digits.

拙訳：超越式（超越関数を含む式？）を計算して予め決められた桁数に正確に丸めるにあたり、どれだけの桁数を余分に持っておく必要があるかを予測する一般的な方法は存在しない。

Table と言っているのは、たぶん教科書の巻末に載っている三角関数表のようなもののことだと思います。下記のようなイメージですね。

$x$	$\roundcirc{\sin(x)}$
$\texttt{1}.\texttt{D037CB2700000}_{(16)}\times 2^{-3}$	$\texttt{1}.\texttt{CC40C37F69D51}_{(16)}\times 2^{-3}$
$\vdots$	$\vdots$
$\texttt{1}.\texttt{D037CB27EE6DF}_{(16)}\times 2^{-3}$	$\texttt{1}.\texttt{CC40C3805229A}_{(16)}\times 2^{-3}$
$\vdots$	$\vdots$
$\texttt{1}.\texttt{D037CB2800000}_{(16)}\times 2^{-3}$	$\texttt{1}.\texttt{CC40C3806348B}_{(16)}\times 2^{-3}$

これの各行について、正確に丸めた値を計算するにあたってどれだけの精度が必要なのかを予測できず、徐々に桁数を上げながら繰り返す*5にしてもどの程度待てばいいかわからないねえという感じでしょう。

たとえば、たかが 53-bit での correctly-rounded な値を得るための途中計算に $10^{3000}$ bits の精度が必要ないとどうして言い切れるでしょうか？ 実際、先の $\sin(x)$ のような例が存在することも予想していなかった人が多いのではないかと思います。 下記のような例もあります。 $$ \begin{aligned} &\phantom{{}={}} \tan(\texttt{1}.\texttt{50486B2F87014}_{(16)}\times 2^{-5}) \\ &= 1. \underbrace{\footnotesize 0101}_{\texttt{5}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0111}_{\texttt{7}} {\footnotesize\,} \underbrace{\footnotesize 1000}_{\texttt{8}} {\footnotesize\,} \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 1110}_{\texttt{E}} {\footnotesize\,} \underbrace{\footnotesize 1011}_{\texttt{B}} {\footnotesize\,} \underbrace{\footnotesize 1111}_{\texttt{F}} {\footnotesize\,} \underbrace{\footnotesize 1111}_{\texttt{F}} {\footnotesize\,} \underbrace{\footnotesize 1001}_{\texttt{9}} {\footnotesize\,} \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 0111}_{\texttt{7}} {\footnotesize\,} \underbrace{\footnotesize 0010}_{\texttt{2}} {\footnotesize\,} \underbrace{\footnotesize 1000}_{\texttt{8}} {\footnotesize\,} \\ & \qquad \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0010}_{\texttt{2}} {\footnotesize\,} {\footnotesize \, \ldots} {\,}_{(2)} \times 2^{-5}. \end{aligned} $$

本題？ 雑談？

double-double だけだと 106-bit 弱になってしまうので、先のような例では負けてしまいます。triple-double ではどうでしょうか？ 足りますか？ わかりません... いかがでしたか？

「100 台程度のワークステーションを使って 4 年かけて最悪ケースの探索をしました」という論文^[2]が 2000 年に出ていますが、下記のように書かれています。

The results [...] give the worst cases for functions $\sin$, $\arcsin$, $\cos$, $\arccos$, $\tan$ and $\arctan$. For those functions, we have worst cases in some bounded domain only, because trigonometric functions are more difficult to handle than the other functions.

拙訳・要約：三角関数は他の関数たちより扱いが難しく、いくつかの区間における最悪ケースしか得られておりません。

そんなあ。

ということで、次からの記事は log や exp の correct rounding な実装の話になるかもしれないし、三角関数関連の話をするかもしれないし、サーベイの旅に出たっきり帰ってこないかもしれません。

Kahan ã‚‚ [6] で触れていますが、Bessel 関数 $J_n(x)$ のような名前のついた非初等関数や、$y^w$ のように複数の引数がある関数についての最悪ケースの探索はほぼ望み薄でしょう。身近なところだと atan2 とかはどうなるんでしょうね。hypot は超越関数ではないのでそこまで悲惨なことにはならないかもしれません。

これはあくまで感情側の話なのですが、標準ライブラリの実装が「ま〜 1–10 ULP くらいの誤差に収まっている気がするしいいでしょ（未証明）」という感じだったり規格側が「別に精度については規定しないよ」という感じだったりするのは困ります。一方で、規格側が「（たとえば）3 ULP までの誤差は許容する」と言ってきても「その値の根拠はどこから...？」となりますし、そもそも「可搬性は...？」となります。値自体の可搬性はないけど事後条件に関しては担保するよ*6みたいな？ 可搬性を第一目標にするなら「この数学関数はこういうアルゴリズムで計算する」とでもすればよいですが、よりよい精度・速度のアルゴリズムが提案されたときにどうするの？という話にもなってきます。

というような旨のことが [7] の文献で書かれていて、次のような主張*7の節があります。

3.2 $\quad$ The only specification that makes sense is correct rounding

なかなか強めの思想な気もしないではないですが、まぁ否定もできないかなという感じです。

といった感じで、三角関数に関してはちゃんとした最悪ケースがまだわかっていないと思っているのですが、correct rounding を謳うライブラリたちはどうやっているのでしょう。Ziv’s rounding test をしながら最大 768-bit くらいの精度まで試すみたいな話はどこかで見たような気はしつつ、それで十分である証明なしには correct と主張できないはずですし。

次数 $d\ge 2$ の代数的数 $\alpha$ に対しては、ある定数 $C_{\alpha}$ が存在して任意の整数 $p$, $q$ ($q\ge 1$) に対して $|\alpha-\tfrac pq| \ge \tfrac{C_{\alpha}}{q^d}$ が成り立つことが Liouville によって示されているらしいですが、超越数は...？という感じで、参っています。どうしましょう？ やはり枝刈り全探索するしか？

たとえば何らかのアルゴリズムによって高々 $2^{-768}$ 程度の相対誤差しかない値を得られるのであれば 53-bit の correct rounding にこだわる必要はないのでは？という主張は真っ当な気がしつつ、まぁそれはそれかなという気持ちでやっています。

参考文献

• [1] Vincent Lefèvre (1997). “An Algorithm that Computes a Lower Bound on the Distance Between a Segment and $\Z^2$.” ftp://ftp.lip.ens-lyon.fr/pub/LIP/Rapports/RR/RR97/RR97-18.ps.Z
• [2] Vincent Lefèvre, Jean-Michel Muller (2000). “Worst Cases for Correct Rounding of the Elementary Functions in Double Precision.” https://inria.hal.science/inria-00072594/document
• [3] Vincent Lefèvre (2000). “The Table Maker’s Dilemma.” http://www.vinc17.org/research/slides/aoc2000-11.pdf
• [4] Vincent Lefèvre, Jean-Michel Muller (2003). “The Table Maker's Dilemma: our search for worst cases.” https://perso.ens-lyon.fr/jean-michel.muller/Intro-to-TMD.htm
• [5] Florent de Dinechin, David Defour, Christophe Lauter (2004). “Fast correct rounding of elementary functions in double precision using double-extended arithmetic.” https://inria.hal.science/inria-00071446/file/RR2004-10.pdf
• [6] William M. Kahan (2004). “A Logarithm Too Clever by Half.” https://people.eecs.berkeley.edu/~wkahan/LOG10HAF.TXT
• [7] Nicolas Brisebarre, Guillaume Hanrot, Jean-Michel Muller, Paul Zimmermann (2025). “Correctly-rounded evaluation of a function: why, how, and at what cost?” https://hal.science/hal-04474530/document

Vincent Lefèvre が大活躍すぎ。また、最後のサーベイ論文に関しては 2025 年（初版は 2024 年）で、びっくりです。

おわり

ということでえびちゃんはこれからお勉強です。フォロヮもごきげんよう。

おわりです。

*1:標準入りしてくれ〜〜〜。

*2:$\sin(x)$ の 54-bit 仮数部への faithful-rounded な値を返すという意味合いでの命名です。

*3:普通の人にとっては $\sin(x) \approx \texttt{1}.\texttt{CC40C3805229B}_{(16)}\times 2^{-3}$ で十分めでたくない？という気もします。普通の人のことは知りませんが。

*4:IEEE 754 の策定にも関わっておられる方で、とてもえらい方ですね。

*5:Ziv’s strategy とか Ziv’s rounding test とか呼ばれるものがあるようです。

*6:人々にはそもそも浮動小数点型ってそういうものだと思われていますか？ う〜〜〜〜ん。

*7:いわゆる数学の意味合いでの Claim ではなくて、お気持ち表明に近いと思います。

浮動小数点数のお手軽チェックリスト

小さな誤差を積み重ねる操作と、積み重なった誤差に影響される操作を区別して考えましょう。

同じ符号の浮動小数点数を加算や乗算を計 $k$ 回繰り返しても、$(1+2^{-53})^k-1 \approx 2^{-53}\cdot k$ 程度の相対誤差しか生まれません。 競プロにおける演算回数はせいぜい $2^{30}$ 回程度でしょうから、だいたい $2^{-23}\approx 1.19\times 10^{-7}$ 程度の相対誤差に収まります。 $2^{20}$ 回なら $2^{-33}\approx 1.16\times 10^{-10}$ 程度です。

これらの操作が問題になることは基本的にないでしょう。 問題になる操作は主に三つです。整数への丸め、比較、引き算です。

具体例を交えながら見ていきます。

前提

IEEE 754 の規格に準拠した処理系を考えます。 C における double のデフォルトの挙動について考え、overflow ã‚„ underflow については起きないこととします。 すなわち、考える浮動小数点数 $x$ は有限であり、下記を満たします。

• $|x|\gt 0$ ならば、ある整数 $2^{52}\le m\lt 2^{53}$ と $-1022\le k\le 1023$ が存在し、$|x| = m\cdot 2^{k-52}$ が成り立つ。

デフォルトの丸め方は tiesToEven と呼ばれるもので、「double で表せる値のうち、与えられた実数に最も近いものを返す」というものです。表せる 2 数のちょうど中間の場合の tiebreak は厳密に定められていますが、込み入らない誤差評価においては問題になってこないので、今回は割愛します。 実数 $x$ を丸めたものを $\roundcirc x$ と書きます。$|x|\gt 0$ のとき下記が成り立ちます。

• 任意の実数 $x$ に対し、実数 $|\varepsilon|\le 2^{-53}$ が存在し、$\roundcirc x = x\cdot(1+\varepsilon)$ が成り立つ。
• 任意の実数 $x$ に対し、実数 $|\varepsilon|\le 2^{-53}$ が存在し、$x = \roundcirc x\cdot(1+\varepsilon)$ が成り立つ。

浮動小数点数における四則演算を $\oplus$, $\ominus$, $\otimes$, $\oslash$ と書き、$x \oplus y = \roundcirc{x+y}$ などが成り立ちます（他の各演算についても同様）。 すなわち、四則演算においては「無限精度で計算したものを正確に丸めたものが得られる」ということです。

たとえば $\tfrac43 = 1.333{\ldots}$ で、これと隣り合う double の値は $$ \begin{aligned} 6004799503160661\times 2^{-52} &= {\small 1.333333333333333}{\footnotesize 259318465024989}{\scriptsize 563971757888793}{\tiny 9453125}, \\ 6004799503160662\times 2^{-52} &= {\small 1.333333333333333}{\footnotesize 481363069950020}{\scriptsize 872056484222412}{\tiny 109375} \end{aligned} $$ です。近いのは前者ですから、$4\oslash 3 = 6004799503160661\times 2^{-52}$ となります。丸めモードを変えない限り、勝手に後者になってしまうことはありません。 特に、厳密に $$ \begin{aligned} 4\oslash 3 &= \frac{4-2^{-52}}3 \\ &= \frac43\cdot(1-2^{-54}) \end{aligned} $$ が成り立ちます。相対誤差が $2^{-53}$ 以下になっていることがわかりますね。

入門編

さて、問題になる二つの操作について見ていきましょう。

整数への丸め

床関数 $\floor x$ を考えます。次のような問題を見てみましょう。

何らかの式で表される実数 $x^{\ast}$ があるとします。これに対して $\floor{x^{\ast}}$ を出力してください。

なお、double で表せる任意の値 $x$ に対して $\roundcirc{\floor x} = \floor x$ が成り立ちます。すなわち、床関数によって誤差が生じることはありません*1。

簡単な例として、$x^{\ast} = 0.1 \times 10$ を考えてみます。 ここで、$\roundcirc{0.1} = \tfrac1{10}\cdot(1+2^{-54})$ かつ $\roundcirc{10} = 10$ です。 $x^{\ast}$ の近似値に相当する浮動小数点数を $x$ として書きます。 $$ \begin{aligned} x &= \roundcirc{0.1} \otimes \roundcirc{10} \\ &= \left(\tfrac1{10}\cdot (1+2^{-54})\right) \otimes 10 \\ &= \roundcirc{1+2^{-54}} \\ &= 1 \end{aligned} $$ です。一方、$0.1\times 10 = 0.1+0.1+\dots+0.1$ として計算した場合を考えてみます。 $$ \begin{aligned} x &= \underbrace{\roundcirc{0.1} \oplus \roundcirc{0.1} \oplus \dots \oplus \roundcirc{0.1}}_{10} \\ &= \left(\tfrac1{10}\cdot (\hphantom{0}1+\hphantom{0}1\cdot 2^{-54})\right) \oplus \underbrace{\roundcirc{0.1} \oplus \roundcirc{0.1} \oplus \dots \oplus \roundcirc{0.1}}_9 \\ &= \left(\tfrac1{10}\cdot (\hphantom{0}2+\hphantom{0}2\cdot 2^{-54})\right) \oplus \underbrace{\roundcirc{0.1} \oplus \roundcirc{0.1} \oplus \dots \oplus \roundcirc{0.1}}_8 \\ &= \left(\tfrac1{10}\cdot (\hphantom{0}3+\hphantom{0}8\cdot 2^{-54})\right) \oplus \underbrace{\roundcirc{0.1} \oplus \roundcirc{0.1} \oplus \dots \oplus \roundcirc{0.1}}_7 \\ &= \left(\tfrac1{10}\cdot (\hphantom{0}4+\hphantom{0}4\cdot 2^{-54})\right) \oplus \underbrace{\roundcirc{0.1} \oplus \roundcirc{0.1} \oplus \dots \oplus \roundcirc{0.1}}_6 \\ &= \left(\tfrac1{10}\cdot (\hphantom{0}5+\hphantom{0}0\cdot 2^{-54})\right) \oplus \underbrace{\roundcirc{0.1} \oplus \roundcirc{0.1} \oplus \dots \oplus \roundcirc{0.1}}_5 \\ &= \left(\tfrac1{10}\cdot (\hphantom{0}6-\hphantom{0}4\cdot 2^{-54})\right) \oplus \underbrace{\roundcirc{0.1} \oplus \roundcirc{0.1} \oplus \dots \oplus \roundcirc{0.1}}_4 \\ &= \left(\tfrac1{10}\cdot (\hphantom{0}7-\hphantom{0}8\cdot 2^{-54})\right) \oplus \roundcirc{0.1} \oplus \roundcirc{0.1} \oplus \roundcirc{0.1} \\ &= \left(\tfrac1{10}\cdot (\hphantom{0}8-12\cdot 2^{-54})\right) \oplus \roundcirc{0.1} \oplus \roundcirc{0.1} \\ &= \left(\tfrac1{10}\cdot (\hphantom{0}9-16\cdot 2^{-54})\right) \oplus \roundcirc{0.1} \\ &= \left(\tfrac{1}{10}\cdot(10-20\cdot 2^{-54})\right) \\ &= 1 - 2^{-53}. \end{aligned} $$ すなわち、$x^{\ast}$ との相対誤差は $-2^{-53}\approx -1.11\times 10^{-16}$ です。（無駄な誤差が出ているとはいえ）十分に近い値が得られていると言って概ね差し支えないでしょう。

しかし、$\floor{x^{\ast}}$ を計算するとなると話は別です。 整数について考えたいときは相対誤差が興味の対象であることはあまりなく、絶対誤差が気になることが多いです。すなわち、$1$ 以上の絶対誤差は許容しないということです。 $\floor{1-2^{-53}} = 0$ ですから、これは許容できません。 相対誤差がどれだけ小さく抑えられていたとしても、$x$ が $\floor{x^{\ast}}$ より少しでも小さいと破綻してしまいます。

$\floor{x^{\ast}}-x$ の上界を $\delta$ とします。$x\oplus\delta\lt \floor{x^{\ast}}+1$ が常に成り立つのであれば、$\floor{x\oplus\delta}$ を計算すればよいです*2。

同符号の計算 $k$ 回程度の場合であれば、$\delta\approx kx^{\ast}\cdot 2^{-53}$ 程度にとってみて、問題ないかを調べるのがよいでしょうか。

一応、「わざわざ掛け算で済むところを足し算にしたからおかしくなった」という主張に対する反例も挙げておきます*3。 $1\oslash 49 = \tfrac1{49}\cdot(1-23\cdot 2^{-58})$ を使います。 $$ \begin{aligned} (1\oslash 49)\otimes 49 &= 1-2^{-53}, \\ \underbrace{(1\oslash 49)\oplus\dots\oplus(1\oslash 49)}_{49} &= 1+6\cdot 2^{-53} \end{aligned} $$ です。これは、足し算をした方のみが（たまたま）うまくいっていますね。

いずれにせよ、「たまたまうまくいくのを祈る」のではなくて、誤差評価をして補正項を入れるなりして「うまくいくことを示してやっている」にするのが重要です。 実際の具体例は後半の章で紹介します。

比較演算

本質的には整数への丸めとほぼ同様だと思います。たとえば本来 $x^{\ast} = 1$ かつ $y^{\ast} = 1$ となるべきところで $x = 1+2^{-52}$ かつ $y = 1$ となってしまっただけで $x\le y$ は成り立たないですね。

いわゆる eps を足し引きして調整するのが小手先テクニックとして広がっていますが、適切な値の決め方は前項の $\delta$ と同様です。

実数 $x^{\ast}$, $y^{\ast}$ のつもりの値として浮動小数点数 $x$, $y$ を得たとします。$x^{\ast} \le y^{\ast}$ の判定のために $x\le y$ とするのは妥当かな？という問題です。なにも考えないと下記のようになり得ます。

	$x\ll y$	$x\le y$	$x\gt y$	$x\gg y$
$x^{\ast} \le y^{\ast}$	true	true	false (?)	N/A
$x^{\ast} \gt y^{\ast}$	N/A	true (?)	false	false

N/A と書いているのは、その欄の条件を満たすような入力が存在しないという意図です。 $x$ と $y$ が十分に離れているときは期待通り判定できますが、そうでないときは誤っちゃうかも？ということですね。

そこで、ある $\delta$ が存在して下記のようになったらうれしいです。

	$x\le y+\delta$	$x\gt y+\delta$
$x^{\ast} \le y^{\ast}$	true	N/A
$x^{\ast} \gt y^{\ast}$	N/A	false

$x$ が $y$ より少しだけ大きかったときは $x^{\ast}\le y^{\ast}$ として判定したいということです。 右上や左下が N/A にならず、下記のようになってしまう場合はうれしくないです。

	$x\le y+\delta$	$x\gt y+\delta$
$x^{\ast} \le y^{\ast}$	true	false (?)
$x^{\ast} \gt y^{\ast}$	N/A	false

$x\le y+\delta$ の判定は $x\ominus y\le \delta$ として行えば、$\tfrac y2\le x\le 2y$ ならば $x\ominus y=x-y$ が成り立つのでよいですね。 これは Sterbenz lemma と呼ばれる有名な補題から従います。 $x\gt 2y$ や $x\lt \tfrac y2$ のときは $x\gg y$ や $x\ll y$ の扱いになっているとうれしいです。

評価の方法としては、整数への丸めの際の $\delta$ と同様に行えばよいですね。 もちろん適切な $\delta$ は状況によります。「テンプレート用におすすめの eps はいくらですか？」のような質問*4はナンセンスであることがわかると思います。

桁落ち

はい、本題の中の本題です。桁落ち (catastrophic cancellation (of significant digits)) と呼ばれるものについて話します。 数が打ち消し合うことを指して cancellation と呼んでいます*5。 主要な桁同士が打ち消し合ったときに誤差の部分が無視できなくなるのが問題となります。 もちろん、浮動小数点数側は「自分は誤差をこのくらい含んでいます」というのを自覚しているわけではないので、「この項はこのくらい誤差を含んでいるねえ」というのは人間側が解析するしかないですね。

実数 $x^{\ast}$, $y^{\ast}$ の近似値として計算した浮動小数点数を $x$, $y$ とします。$x^{\ast} = x\cdot(1+\varepsilon_x)$ かつ $y^{\ast} = y\cdot(1+\varepsilon_y)$ としましょう。$|\varepsilon_x|$, $|\varepsilon_y|$ の上界は状況によるものです。

このとき、$x^{\ast}-y^{\ast} = (x-y) + (x\varepsilon_x - y\varepsilon_y)$ となります。$x-y$ がとても小さいと、$x\varepsilon_x - y\varepsilon_y$ が無視できなくなってきます。

Exercise 1: $x = 2^{52}+2$ かつ $y = 2^{52}+1$ に対して $\sqrt x-\sqrt y$ を計算せよ。

note: $\sqrt{1+x} = 1 + \tfrac12x - \tfrac18x^2 + \tfrac1{16}x^3 + O(x^4)$ が成り立つ。

Bad solution

$$ \begin{aligned} \sqrt{2^{52}+2} &= 2^{26}\cdot\sqrt{1+2^{-25}} \\ &= 2^{26}\cdot(1 + 2^{-26} - 2^{-28} + 2^{-29} + \cdots), \\ \sqrt{2^{52}+1} &= 2^{26}\cdot\sqrt{1+2^{-26}} \\ &= 2^{26}\cdot(1 + 2^{-27} - 2^{-29} + 2^{-30} + \cdots) \\ \end{aligned} $$ なので $$ \begin{aligned} \roundcirc{\sqrt{2^{52}+2}} &= 2^{26} + 2^{-26}, \\ \roundcirc{\sqrt{2^{52}+1}} &= 2^{26} \end{aligned} $$ が成り立つ。よって $$ \begin{aligned} \roundcirc{\sqrt{2^{52}+2}} \ominus \roundcirc{\sqrt{2^{52}+1}} &= (2^{26}+2^{-26})\ominus 2^{-26} \\ &= 2^{-26} \end{aligned} $$ である。よって $\sqrt x-\sqrt y \approx 2^{-26}$ である (?)。$\eod$

Better solution

$$ \begin{aligned} &\phantom{{}={}} \sqrt{2^{52}+2} - \sqrt{2^{52}+1} \\ &= \frac{(\sqrt{2^{52}+2} - \sqrt{2^{52}+1})\cdot(\sqrt{2^{52}+2} + \sqrt{2^{52}+1})}{\sqrt{2^{52}+2} + \sqrt{2^{52}+1}} \\ &= \frac{(2^{52}+2) - (2^{52}+1)}{\sqrt{2^{52}+2} + \sqrt{2^{52}+1}} \\ &= \frac1{\sqrt{2^{52}+2} + \sqrt{2^{52}+1}} \end{aligned} $$ であり、 $$ \begin{aligned} 1\oslash (\roundcirc{\sqrt x}\oplus \roundcirc{\sqrt y}) &= 1\oslash ( (2^{26}+2^{-26})\oplus 2^{26} ) \\ &= 1\oslash 2^{27} \\ &= 2^{-27} \end{aligned} $$ である。よって $\sqrt x-\sqrt y \approx 2^{-27}$ である。また、相対誤差は高々 $4\cdot 2^{-53}$ 程度である。$\eod$

Exercise 2: $x =2^{-53}$ に対して $\tfrac1x\cdot(e^x-1)$ を計算せよ。

Bad solution

$$ \exp(2^{-53}) = {\small 1.000000000000000}{\footnotesize 111022302462515}{\scriptsize 660205338988848}{\tiny 236989105051344{\dots}} $$ であり、 $$ \roundcirc{\exp(2^{-53})} = 1 + 2^{-52} $$ である*6。実際に計算し、 $$ \begin{aligned} &\phantom{{}={}} (\roundcirc{\exp(2^{-53})}\ominus 1)\oslash 2^{-53} \\ &= 2^{-52}\oslash 2^{-53} \\ &= 2 \end{aligned} $$ を得る。よって $\tfrac1x\cdot(e^x-1) \approx 2$ である (?)。$\eod$

Better solution（むずかしい）

$e^x-1$ を計算する際の桁落ちを避けたい。

$y = \roundcirc{\exp(x)}$ とする。ここで、ある実数 $x'$ に対して $y = \exp(x')$ が成り立つ。 これにより、$\roundcirc{\exp(x)}\ominus 1 = \exp(x')-1$ が成り立つため、$\tfrac1x\cdot(e^x-1)$ の代わりに $\tfrac1{x'}\cdot(e^{x'}-1)$ を計算することにする。

数値計算により $\roundcirc{\log(y)} = 2^{-52} - 2^{-105}$ を得る。よって、 $$ \begin{aligned} (y\ominus 1)\oslash \roundcirc{\log(y)} &= 2^{-52} \oslash (2^{-52} - 2^{-105}) \\ &= 1+2^{-52} \end{aligned} $$ となる。

さて、誤差評価をがんばる。 $$ \begin{aligned} \exp(x) &= \sum_{i=0}^{\infty} \frac{x^i}{i!}, \\ \frac{\exp(x)-1}x &= \sum_{i=0}^{\infty} \frac{x^i}{(i+1)!} \end{aligned} $$ であり、 $$ \begin{aligned} &\phantom{{}={}} \frac{\exp(x+h)-1}{x+h} - \frac{\exp(x)-1}x \\ &= \sum_{i=0}^{\infty} \frac{(x+h)^i}{(i+1)!} - \sum_{i=0}^{\infty} \frac{x^i}{(i+1)!} \\ &= \sum_{i=0}^\infty \frac1{(i+1)!}\sum_{j=1}^i \binom ij h^j x^{i-j} \\ % &= \frac1{2!} \left(\binom 11 h^1 x^0 \right) \\ % & \qquad\quad {} + \frac1{3!} \left(\binom 21 h^1 x^1 + \binom 22 h^2 x^0 \right) \\ % & \qquad\quad {} + \frac1{4!} \left(\binom 31 h^1 x^2 + \binom 32 h^2 x^1 + \binom 33 h^3 x^0 \right) \\ % & \qquad\quad {} + \cdots \\ % &= \left(\frac1{2!} \binom11 h^1 + \frac1{3!} \binom 22 h^2 + \frac1{4!} \binom 33 h^3 + \cdots \right) x^0 \\ % & \qquad\quad {} + \left(\frac1{3!} \binom21 h^1 + \frac1{4!} \binom 32 h^2 + \cdots \right) x^1 \\ % & \qquad\quad {} + \cdots \\ &= \sum_{i=0}^{\infty} \sum_{j=0}^{\infty} \frac1{(i+j+2)!} \binom{i+j+1}{j+1} h^{j+1} x^i \\ &= \sum_{j=0}^{\infty} \frac1{(j+2)!} h^{j+1} + \sum_{j=0}^{\infty} \frac{j+2}{(j+3)!} h^{j+1} x + O(x^2) \\ &= \left(\frac{\exp(h)-(1+h)}h\right) + \left(\frac{\exp(h)\cdot(h-1)+1}{h^2}-\frac12\right)x + O(x^2) \\ % &= \left(\frac{\tfrac12h^2 + O(h^3)}h\right) + \left(\frac{(1+h+\tfrac12h^2+O(h^3) )\cdot(h-1)+1}{h^2}-\frac12\right)x + O(x^2) \\ % &= \left(\frac12h + O(h^2)\right) + \left(\frac{(h+h^2+\tfrac12h^3 - 1-h-\tfrac12h^2-\tfrac16h^3 + O(h^4) )+1}{h^2}-\frac12\right)x + O(x^2) \\ % &= \left(\frac12h + O(h^2)\right) + \left(\frac{(\tfrac12h^2+\tfrac13h^3 + O(h^4) )}{h^2}-\frac12\right)x + O(x^2) \\ &= \left(\tfrac12h + O(h^2)\right) + \left(\tfrac13h + O(h^2)\right)x + O(x^2) \\ &= \tfrac h2 + O(h(x+h) ) \end{aligned} $$ になると思う。また、$\tfrac1{x'}(y-1)$ と $(y\ominus 1)\oslash \roundcirc{\log(y)}$ の相対誤差は $2\cdot 2^{-53}$ 程度である。真の値は $1$ 程度であるから、絶対誤差も $2\cdot 2^{-53}$ 程度である。 $h$ は高々 $2^{-53}$ 程度なので、合わせて $2.5\cdot 2^{-53}$ 程度の絶対誤差であることがわかる。$\eod$

Sterbenz lemma は、$x$ と $y$ の比が $\tfrac12$ 以上 $2$ 以下のときは $x\ominus y = x-y$ になる（差による誤差が生じない）ということを主張しています。これはあくまで「入力の浮動小数点数からの誤差がない」というだけで、入力の浮動小数点数が誤差を含む場合は「出力もそのままその誤差を含む」ということになります。 逆に、入力の誤差がない場合は誤差なく引き算できるので、「近い値の引き算は怖いから...」などと変な怯え方をする必要はありません。

このような、入力に誤差がないことがわかっている状況での cancellation は、benign cancellation と呼ばれることがしばしばあります。 日本語での定訳があるかはわかりません。

実践編

解析のやり方とかも含めて追っていきましょう。コンテスト中でもできるように、手軽な感じの手法を心がけます。

整数への丸め

reference: AGC 047 A - Integer Product

整数部分が 4 桁以下、小数部分が 9 桁以下であるような実数 $A_i$ が与えられます。解法の詳細には触れませんが、$10^9\cdot A_i$ を計算したくなります。 しかし、$\floor{10^9\otimes \roundcirc{A_i}}$ として計算すると失敗するような例はいくらでもあります。 $$ \begin{aligned} \roundcirc{1024.000000006} &= {\small 1024.000000005999936}{\footnotesize 547596007585525}{\scriptsize 5126953125} \\ &= 4503599627396884\cdot 2^{-42}, \\ \roundcirc{1024.000000006} \otimes 10^9 &= {\small 1024000000005.9998779296875} \\ &= 8388608000049151 \cdot 2^{-13}. \end{aligned} $$

ある $\delta$ を導入して $\floor{(10^9\otimes \roundcirc{A_i})\oplus\delta}$ を考える方法を考えましょう。

以下、$x = \roundcirc{A_i}$ とします。 $$ \begin{aligned} |x - A_i| &\le x\cdot 2^{-53}, \\ |10^9\otimes x - 10^9\cdot x| &\le 10^9\cdot x\cdot 2^{-53} \end{aligned} $$ が成り立つので、三角不等式より $$ \begin{aligned} |10^9\otimes x - 10^9\cdot A_i| &= |10^9\otimes x - 10^9\cdot x| + |10^9\cdot x - 10^9\cdot A_i| \\ &\le 10^9\cdot x\cdot 2^{-53} + 10^9\cdot x\cdot 2^{-53} \\ &\lt 2\cdot 10^9\cdot 10^5\cdot 2^{-53} \\ &\lt 0.022204\\ &\lt \roundcirc{0.02221} \end{aligned} $$ が成り立ちます。よって、 $$ 10^9\cdot A_i - 0.022204 \lt 10^9\otimes x \lt 10^9\cdot A_i + 0.022204 $$ であり、 $$ \begin{aligned} 10^9\cdot A_i \lt 10^9\otimes x + \roundcirc{0.02221} &\lt 10^9\cdot A_i + 0.022204 + \roundcirc{0.02221} \\ &\lt 10^9\cdot A_i + 1 \end{aligned} $$ です。また、$10^9\cdot A_i\lt 2^{52}$ から、$(10^9\otimes x) \oplus \roundcirc{0.02221} \lt 10^9\cdot A_i + 1$ が成り立ちます。 よって、$\floor{(10^9\otimes x) \oplus \roundcirc{0.02221}} = 10^9\cdot A_i$ が従います。わざわざ細かい桁まで指定せず、$\delta = \roundcirc{0.03}$ として問題ないでしょう。

note: たとえば $2^{52} + \roundcirc{0.6} \lt 2^{52}+1$ ですが、$2^{52}\oplus\roundcirc{0.6} = 2^{52}+1$ なので、そのような部分が問題ないかも気をつける必要があります。

note: 差が $0.5$ 未満で抑えられていることから、.round() を用いる解法も正当であることがわかりますね。

提出コード：提出 #65357163

---

次の問題です。

reference: ABC 169 C - Multiplication 3

これは double だけで解けるでしょうか？

たとえば $(A, B) = (999689151469700, 9.01)$ などのケースで $\floor{AB}$ 自体を double で表すことができないので、「答えを double で返す」という前提では厳しそうですね。 $100B$ ã‚’ $901$ 以上 $999$ 以下の奇数として固定し、 $$ A = \left(2\Floor{\frac12\Ceil{\frac{2^{53}}{100B}}}+1\right)\cdot 100 $$ とすることで、所望の $(A, B)$ を構成することができます。

このように、「そもそも答えを double で表せるのか？」と考えるのは重要です。$2^{53}+1$ 以上の整数は double で表せるとは限らないので、「答えがそういう値になる入力を構成できるか？」という視点で考えてみるのがよいですね。

先の例と同様にして xx.yy を補正しつつ 100 倍するなどして xxyy にしてから整数型として扱うのが賢明でしょう*7。

仮数部が 64-bit ã‚„ 113-bit の浮動小数点型を使える言語であれば、それを使えば可能な気はします。もちろん素朴に floor(a * b) とするだけではだめですが、前述のように補正項を入れるなどで対応することはできますね。「誤差とかよくわかんないけど long double を使って強引に通した」ではなくて「long double だと仮数部のビット長が足りるので long double を使いました」と言えるようになってほしいものです*8。

ということで、ここでは、64-bit 仮数部の型を使いましょう。 求めたいのは $\floor{AB}$ で、計算可能なのは $A\otimes\roundcirc B$ です。

実数 $|\varepsilon_B|, |\varepsilon_{\otimes}|\le 2^{-64}$ が存在して $$ A\otimes\roundcirc B = \left(A\cdot (B\cdot (1+\varepsilon_B) )\right)\cdot (1+\varepsilon_{\otimes}) $$ が成り立つので、 $$ \begin{aligned} |A\otimes\roundcirc B - AB| &= AB\cdot ( (1+2^{-64})^2-1) \\ &\lt 10^{15}\cdot 9.99\cdot 1.085\times 10^{-19} \\ &\lt 1.084\times 10^{-3} \end{aligned} $$ が成り立ちます。

今回は、$AB$ が整数とは限らないため、先のものよりも上界はぎりぎりです。 $$A\otimes \roundcirc B - \floor{AB} \in (-0.001084\lldot 0.991084)$$ であり、 $$(A\otimes \roundcirc B)\oplus \roundcirc{0.002} - \floor{AB} \in (0.00092\lldot 0.99309)$$ なので、$\floor{(A\otimes \roundcirc B)\oplus \roundcirc{0.002}}$ を計算すればよいですね。

提出コード：提出 #65358049

比較演算

大まかな考え方は前項と同じでしょうから、ここで特に話すことはないような気がしています。気が変わったら追記します。

桁落ち

二つの実数を異なる式で計算し、「本来それらは差が $0$ なので問題ないはずだ」というのを前提としているような例では、桁落ちで撃墜できることが多いです。

減算をしない形*9を意識して式変形を行うのが重要でしょう。 分数であれば、分子 $n$ と分母 $d$ を整数で計算し、$\roundcirc n\oslash \roundcirc d$ として計算するのがよさそうです。

実数 $|\varepsilon_n|, |\varepsilon_d|, |\varepsilon_{\oslash}|\le 2^{-53}$ が存在して $$ \roundcirc n\oslash \roundcirc d = \frac{n\cdot(1+\varepsilon_n)}{d\cdot(1+\varepsilon_d)}\cdot(1+\varepsilon_{\oslash}) $$ と書けるので、 $$ \begin{aligned} |(\roundcirc n\oslash \roundcirc d) - (n\div d)| &= \frac nd\cdot \left| \frac{(1+\varepsilon_n)(1+\varepsilon_{\oslash})}{1+\varepsilon_d} -1 \right| \\ &\le \frac nd\cdot \left| \frac{(1+2^{-53})^2}{1-2^{-53}} -1 \right| \\ % &= \frac nd\cdot \left| 1 - (1+2^{-53})^2 \cdot \sum_{i=0}^{\infty} {(2^{-53})^i} \right| \\ % &= \frac nd\cdot \left| 1 - (1+2^{-53})^2 \cdot \left(1+2^{-53}+\sum_{i=2}^{\infty} {(2^{-53})^i}\right) \right| \\ % &= \frac nd\cdot \left| 1 - (1+2^{-53})^3 - (1+2^{-53})^2\cdot \sum_{i=2}^{\infty} {(2^{-53})^i} \right| \\ % &= \frac nd\cdot \left| 1 - (1+2^{-53})^3 - (1+2^{-53})^2\cdot \frac{(2^{-53})^2}{1-2^{-53}} \right| \\ % &= \frac nd\cdot \left| 1 - (1+2^{-53})^3 - (2^{-53})^2\cdot \frac{(1+2^{-53})^2}{1-2^{-53}} \right| \\ &= \frac nd\cdot \left( (1+2^{-53})^3 - 1 + (2^{-53})^2\cdot \frac{(1+2^{-53})^2}{1-2^{-53}} \right) \\ &\lt \frac nd\cdot (3+10^{-15})\cdot 2^{-53} \end{aligned} $$ となり、$3\cdot 2^{-53}$ 程度の相対誤差で抑えられます。

分子や分母に根号が含まれる場合は、分子と分母の適切な方の有理化を行うことで減算を避けるのがよさそうです。それ以外の無理数がある場合は、たぶん困ります。 三角関数であれば和積公式などが使える場合もあるかも。

具体的に撃墜するケースの構築や式変形の方法の例については、別の記事で過去に書いたのでそちらを参考にしてください。

rsk0315.hatenablog.com

同符号のたくさんの積和

例として期待値の DP とかを考えてみましょう。

reference: ABC 402 E - Payment Required

まず最初に、整数 $p$ に対しての $(100 - p) \oslash 100$ と $1 \ominus (p\oslash 100)$ の誤差の違いについて考えてみたりしますか？ 前者は誤差の出る操作が 1 回なので $2^{-53}$ 程度の相対誤差なのに対し、後者は catastrophic cancellation の影響でもっと大きな誤差が出そうですね。実際*10、 $$ \begin{aligned} % 1\oslash 100 &= \frac{1+3\cdot 2^{-57}}{100}, \\ 99\oslash 100 &= \frac{99-2^{-50}}{100} \end{aligned} $$ なので、Sterbenz lemma を使いつつ $$ \begin{aligned} 1\ominus (99\oslash 100) &= 1 - \frac{99-2^{-50}}{100} \\ &= \frac{1+2^{-50}}{100} \\ &= 0.01 \cdot (1+8\cdot 2^{-53}) \end{aligned} $$ となり、本来の上界の $8$ 倍程度の誤差が無駄に出ています。実際にこの誤差が問題になることは少ないでしょうが、わざわざ誤差の出る方の形を好む必要もないでしょうね。

さて、DP で計算される値（ここでは期待値）の誤差はどうなるでしょう？ 雑に見積もっても、遷移は $2^8\cdot 5000$ 回で抑えられます、遷移は $$ \mathrm{dp}[i] = (p_0\otimes \mathrm{dp}[i_0]) \oplus (p_1\otimes \mathrm{dp}[i_1]) $$ のような式になり、$p_0$ や $p_1$ にそれぞれ計算 $1$ 回分の相対誤差があるので、遷移あたり計算 $5$ 回分の相対誤差で抑えられます*11。

ということで、計 $2^8\cdot 5000\cdot 5 = 6.4\times 10^6$ 回分になり、$2^{-53}\cdot 6.4\times 10^6 \approx 7.1\times 10^{-10}$ 程度の相対誤差で抑えられることがわかります。許容誤差は $10^{-6}$ なので問題なさそうですね。

おまけ

二分探索

よく $x_U - x_L \ge \roundcirc{10^{-9}}$ などをループ継続条件にして失敗する人がいますね。 たとえば $(x_L, x_U) = (10^8, x_L+2^{-26})$ を考えてみます。 $x_L \lt x\lt x_U$ を満たすような $x$ は存在しませんが、$2^{-26}\gt 10^{-9}$ が成り立ちます。 このとき $(x_L\oplus x_U)\oslash 2 = x_L$ なので、更新式の条件的に $x_U \gets x_L$ で更新されない場合は無限ループとなります。

そこで、よく「浮動小数点数で二分探索をするときは、差で判定するんじゃなくて回数でループするといいんだよ」ということが言われています。 果たして回数でのループで問題ないことは明らかでしょうか？ その回数はどの程度？ 得られる精度はどのくらい？ $[x_L\lldot x_U]$ の範囲でない値がうっかり返ってしまうことはない？ など、「それっぽいから大丈夫そう」でやっている人はやっぱり多いのではないでしょうか。

$|x\oplus y| \lt \infty$ ならば $(x\oplus y)\oslash 2 = \roundcirc{\tfrac{x+y}2}$ が成り立ちます。 よって、$x\le y$ ならば $x \le (x\oplus y)\oslash 2\le y$ となりますね。

たとえば $(x_L\oplus x_U)\oslash 2\notin \{x_L, x_U\}$ をループ継続条件にしてもよさそうな気がしていますが、どうでしょう。 どのくらいの回数で収束してくれるでしょうか。

今回の記事のスコープとは微妙に外れている気がするので、次回以降の記事で気が向いたら解析してみることにします。

練習問題

浮動小数点数や誤差が関連している問題を挙げておきます。必ずしも「double で解けますよ〜」という意味ではありません。その見極めも含めて使ってもらえたらいいんじゃないかなと思います。

ABC 300 から ABC 403 で、問題文中に「誤差」が含まれている問題のリストです。

• ABC313 F
• ABC314 E
• ABC314 Ex
• ABC315 F
• ABC319 C
• ABC324 F
• ABC327 E
• ABC332 E
• ABC334 F
• ABC341 G
• ABC342 F
• ABC350 E
• ABC356 G
• ABC374 D
• ABC375 B
• ABC382 E
• ABC385 F
• ABC392 D
• ABC393 G
• ABC401 G
• ABC402 E

739 問中 21 問しかないらしいです。まじ？ 2.84 % くらいです。

「小数」のリストです。これはあまり関係ないですね。

• ABC314 A
• ABC333 G
• ABC354 D
• ABC367 B
• ABC397 A

また、それ以外で、手元のコードで f64 が含まれていたものに基づくリストです。あまり素直じゃないものが多かったです。

• ABC308 C
• ABC330 C
• ABC361 F
• ABC390 B
• ABC397 A
• ABC400 B
• ABC400 C
• ABC400 E
• ABC402 D

上記に当てはまらないが「初心者が浮動小数点数で解こうとしてハマりそうなもの」とかがあれば募集中です。

あとがき

最近は correct rounding モノばかり書いていたので、「えっ　今日は全員誤差出していいのか‼️」「ああ…しっかり出せ」「桁落ちもいいぞ！（？？）」のような気持ちになりました。「ただ今より整数訓練を開始する‼️」（？？？）

浮動小数点数に限らないことだとは思うのですが、土壇場で必要になったときに日頃から訓練していないと誤ってしまうことが多いと思います。 コンテスト本番でわざわざ浮動小数点数による解法を好んで使おうとする理由はあまりないとも思いますが、復習の際には敢えて「浮動小数点数で解けないもんかな？」とやってみることも重要かなとも思います。

初心者層の優先度としてはもっと他にやるべきことがあるとは思いつつ、「誤差とかよくわかんないけどなんか通った」「なんかわかんないけど落ちた」と言い続けるよりはいいんじゃないかな？とも思います。Twitter で検索してみるとそういう人はしばしばいます。

う〜〜〜ん、少し考えてみても、「このくらいの層であれば他の勉強よりも浮動小数点数の勉強を優先するのがよい」と勧めたくなる状況が思いつかないですね。 気分転換としてやるのがよいのではないでしょうか。

初心者が望んでいる浮動小数点数の解説記事ってどんなのですか？ 正直あまりわからなくなってきました。もちろん初心者以外が望んでいるものでもいいんですが、あまりわからないです。わからないのでえびちゃんが好き勝手に書いちゃいます。

いま GW 何日目ですか？ あまり現実逃避に記事を書いて時間を溶かさない方がいいですよ。correct rounding から逃げてはいけません。

おわり

おわりです。

*1:もちろんライブラリの実装にはよるのですが、容易に実現できるため、まともな実装であることを仮定するのは十分現実的でしょう。

*2:$\floor{x+\delta}$ ではなく $\floor{x\oplus\delta}$ を考える必要があることに注意。

*3:$x=1$ とすればいいのにわざわざ四則演算をしているからおかしくなった ← それはそう。

*4:こういう質問をする人が昔しばしばいたような気もするし、架空の人だったような気もします。老人の記憶はアテになりません。

*5:打ち消しの操作を主語とする文では、たとえば “In this situation, the cancellation is catastrophic.” とか言われたりします。他にも “the terms are cancelled catastrophically.” と書かれていたとしても桁落ちのことを指すでしょう。catastrophic の代わりに severe と言われることもあります。

*6:correct rounding な実装でないとここで $1$ を返すかもしれません。correct rounding でない実装を前提にしても仕方がないですし、correct rounding な実装を前提にしていても桁落ちでこわれますよ〜というのを話したいので、大勢には影響がないと思います。

*7:浮動小数点数を避けるのが賢明？ まぁそういう立場もあるでしょうね。

*8:答えを表すのに必要なビット長があれば十分という意味ではなくて、計算途中の精度なども当然加味した上での主張です。

*9:もちろん加算だけの形であっても異符号になりうるなら無意味です。

*10:こういうのって読者への課題として残した方がいいですか？

*11:$\oplus$ の部分をちゃんと考えると $3$ 回分で抑えられる気がしますが、ここでは不要だと思うので大雑把にいきます。

correct rounding への道 (4) error-free transformation

倍精度浮動小数点型であるところの 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
  • https://link.springer.com/article/10.1007/BF01397083

あとがき

今回は、double-double や triple-double を実装する際の基本操作のみを導入しました。

非正規化数についても触れましたが、「よくわからないので避けたいもの」という印象は拭えても「場合分けが面倒になるもの」という印象が強まってしまっているかもしれません？ わからないので避けたいというのよりはマシな気はしますね。

次回いきなり気が変わって double-double ではないなにかの話を始めるかもしれません？ そのときの気分でテーマが決まるので仕方ないですね。

毎度のことですが、Sterbenz lemma が大活躍すぎます。「$\roundcirc x=x$ を示したくなったらまず Sterbenz lemma を使えないか疑え」まである気がします。 もちろん、mod の性質など、それ以外のものを使う局面もしばしばありますが、Sterbenz lemma がお手軽どころの頻出パターンという感じですね。

なんか疲れてしまいました。証明の方針が全然思いつかなかったり、場合分けをがんばる方法で書き終えたと思ったらまとめられることに気づいたり、なにも考えられなくなったり、シンプルに嘘を書いていたり、ひ〜という感じです。 お気持ちパートも流れるように書けなくなってきました。不調？

さっさと元気になって続編を勉強したいですね。correct rounding に興味のあるフォロヮがどの程度いるかは知らないですが、えびちゃんは知りたがっているので。

おわり

おわりです。

*1:もちろん、考えなしにそういうことをすると二重丸めなどのうれしくないことが発生しえますが、それはまた今後の回で扱う予定です。

*2:別に標準で用意されていることにこだわる必要もないにはないですね。

*3:当初は今回で終わらせる予定だったのですが、思っていたよりボリュームが大変なことになって収集がつかなくなったためです。

*4:$k=-1023$ のような見た目をしているの微妙かな？ 後の回で変えるかもしれないです。

*5:$x\gt 0\gt y$ のときは $\frac x2\le y\le 2x$ は偽である。

5.0 * x == x + x + x + x + x

#浮動小数点数 に関するクイズ（マジメなやつ）をつくってみた😃 pic.twitter.com/2aYU8zcguU

— 某ZR（ざんねん🙃） (@zr_tex8r) 2025å¹´4月8æ—¥

「反例がない気はしないんだが、具体的にぱっと構成するのは難しい」みたいなシリーズ

— えびちゃん🍑🍝🦃 (@rsk0315_h4x) 2025å¹´4月8æ—¥

自分は
a + a + a + a + a == a * 5.0
が判断できなくて困ってる（ざんねん🙃） https://t.co/b32AzZcZQ8

— 某ZR（ざんねん🙃） (@zr_tex8r) 2025å¹´4月8æ—¥

証明

binary64, tiesToEven で考えます。以下、$x$ は有限の浮動小数点数とします。特に、$x$ は NaN ではありません*1。

Claim 1: 任意の $x$ に対し、$x \oplus x = 2 \otimes x = 2x$ が成り立つ。

Proof: 明らか。$\qed$

Claim 2: 任意の $x$ に対し、$(x \oplus x)\oplus x = 3 \otimes x$ が成り立つ。

Proof

$$ \begin{aligned} (x\oplus x)\oplus x &= 2x \oplus x \\ &= \roundcirc{2x+x} \\ &= \roundcirc{3x} \\ &= 3\otimes x. \quad\qed \end{aligned} $$

Claim 3: 任意の $x$ に対し、$( (x\oplus x)\oplus x)\oplus x = 4\otimes x = 4x$ が成り立つ。

Proof

$4\otimes x = 4x$ は明らか。$(3\otimes x) \oplus x = 4x$ を示せばよい。

$x = 0$ の場合は明らかなので、$x\gt 0$ について考える。一般性を失わず $2^{52}\le x\lt 2^{53}$ とする。 ある実数 $|\varepsilon|\le 2^{-53}$ が存在して $3\otimes x = 3x + \hfloor{3x}\cdot\varepsilon$ が成り立つ。

$1.5\cdot 2^{53}\le 3x\lt 1.5\cdot 2^{54}$ より、$\hfloor{3x} = 2^{53}$ または $\hfloor{3x} = 2^{54}$ が成り立つ。

Case 1: $\hfloor{3x} = 2^{53}$

$\varepsilon' = \hfloor{3x}\cdot\varepsilon$ とすると $|\varepsilon'|\le 1$ が成り立つ。

$$ \begin{aligned} (3\otimes x)\oplus x &= (3x + \hfloor{3x}\cdot\varepsilon) \oplus x \\ &= (3x + \varepsilon')\oplus x \\ &= \roundcirc{3x + \varepsilon' + x} \\ &= \roundcirc{4x + \varepsilon'}. \end{aligned} $$

Case 1-a: $x = 2^{52}$

$4x$ 未満で最大の浮動小数点数は、$4x-2$ である。$\roundcirc{4x-1} = 4x$ であるから、$\roundcirc{4x+\varepsilon'}\ge 4x$ である。 また、$4x$ より大きい最小の浮動小数点数は $4x+4$ であるから、$\roundcirc{4x+\varepsilon'}\le 4x$ である。 よって、$\roundcirc{4x+\varepsilon'} = 4x$ である。

Case 1-b: $x \gt 2^{52}$

$4x$ 未満で最大の浮動小数点数は $4x-4$、$4x$ より大きい最小の浮動小数点数は $4x+4$ であるから、Case 1-a 同様にして $\roundcirc{4x+\varepsilon'} = 4x$ である。

Case 2: $\hfloor{3x} = 2^{54}$

$\varepsilon' = \hfloor{3x}\cdot\varepsilon$ とすると $|\varepsilon'|\le 2$ が成り立つ。 また、$2^{54}\le 3x\lt 4x\lt 2^{55}$ であることに注意する。

任意の整数 $2^{54}\le m\lt 2^{55}$ について、 $$ \roundcirc m = \begin{cases} m, & \text{if}\: m \equiv 0 \pmod 8; \\ m-1, & \text{if}\: m \equiv 1 \pmod 8; \\ m-2, & \text{if}\: m \equiv 2 \pmod 8; \\ m+1, & \text{if}\: m \equiv 3 \pmod 8; \\ m, & \text{if}\: m \equiv 4 \pmod 8; \\ m-1, & \text{if}\: m \equiv 5 \pmod 8; \\ m+2, & \text{if}\: m \equiv 6 \pmod 8; \\ m+1, & \text{if}\: m \equiv 7 \pmod 8 \\ \end{cases} $$ が成り立つ。よって $$ 3\otimes x = \roundcirc{3x} = \begin{cases} 3x, & \text{if}\: 3x \equiv 0 \pmod 8; \\ 3x-1, & \text{if}\: 3x \equiv 1 \pmod 8; \\ 3x-2, & \text{if}\: 3x \equiv 2 \pmod 8; \\ 3x+1, & \text{if}\: 3x \equiv 3 \pmod 8; \\ 3x, & \text{if}\: 3x \equiv 4 \pmod 8; \\ 3x-1, & \text{if}\: 3x \equiv 5 \pmod 8; \\ 3x+2, & \text{if}\: 3x \equiv 6 \pmod 8; \\ 3x+1, & \text{if}\: 3x \equiv 7 \pmod 8 \\ \end{cases} $$ であり、$3^{-1} \equiv 3 \pmod 8$ を踏まえて整理して、 $$ \begin{aligned} 3\otimes x &= \begin{cases} 3x, & \text{if}\: x \equiv 0 \pmod 8; \\ 3x-1, & \text{if}\: x \equiv 3 \pmod 8; \\ 3x-2, & \text{if}\: x \equiv 6 \pmod 8; \\ 3x+1, & \text{if}\: x \equiv 1 \pmod 8; \\ 3x, & \text{if}\: x \equiv 4 \pmod 8; \\ 3x-1, & \text{if}\: x \equiv 7 \pmod 8; \\ 3x+2, & \text{if}\: x \equiv 2 \pmod 8; \\ 3x+1, & \text{if}\: x \equiv 5 \pmod 8 \\ \end{cases} \\ &= \begin{cases} 3x, & \text{if}\: x \equiv 0 \pmod 8; \\ 3x+1, & \text{if}\: x \equiv 1 \pmod 8; \\ 3x+2, & \text{if}\: x \equiv 2 \pmod 8; \\ 3x-1, & \text{if}\: x \equiv 3 \pmod 8; \\ 3x, & \text{if}\: x \equiv 4 \pmod 8; \\ 3x+1, & \text{if}\: x \equiv 5 \pmod 8; \\ 3x-2, & \text{if}\: x \equiv 6 \pmod 8; \\ 3x-1, & \text{if}\: x \equiv 7 \pmod 8. \\ \end{cases} \end{aligned} $$

よって、 $$ (3\otimes x) + x = \begin{cases} 4x, & \text{if}\: x \equiv 0 \pmod 8; \\ 4x+1, & \text{if}\: x \equiv 1 \pmod 8; \\ 4x+2, & \text{if}\: x \equiv 2 \pmod 8; \\ 4x-1, & \text{if}\: x \equiv 3 \pmod 8; \\ 4x, & \text{if}\: x \equiv 4 \pmod 8; \\ 4x+1, & \text{if}\: x \equiv 5 \pmod 8; \\ 4x-2, & \text{if}\: x \equiv 6 \pmod 8; \\ 4x-1, & \text{if}\: x \equiv 7 \pmod 8 \\ \end{cases} $$ である。さて、 $$ \begin{aligned} x\equiv 0 \pmod 8 &\implies 4x \equiv 0 \pmod 8, \\ x\equiv 1 \pmod 8 &\implies 4x+1 \equiv 5 \pmod 8, \\ x\equiv 2 \pmod 8 &\implies 4x+2 \equiv 2 \pmod 8, \\ x\equiv 3 \pmod 8 &\implies 4x-1 \equiv 3 \pmod 8, \\ x\equiv 4 \pmod 8 &\implies 4x \equiv 0 \pmod 8, \\ x\equiv 5 \pmod 8 &\implies 4x+1 \equiv 5 \pmod 8, \\ x\equiv 6 \pmod 8 &\implies 4x-2 \equiv 6 \pmod 8, \\ x\equiv 7 \pmod 8 &\implies 4x-1 \equiv 3 \pmod 8 \end{aligned} $$ であるから、 $$ \begin{aligned} (3\otimes x) \oplus x &= \roundcirc{(3\otimes x)+x} \\ &= \begin{cases} 4x, & \text{if}\: x \equiv 0 \pmod 8; \\ (4x+1)-1, & \text{if}\: x \equiv 1 \pmod 8; \\ (4x+2)-2, & \text{if}\: x \equiv 2 \pmod 8; \\ (4x-1)+1, & \text{if}\: x \equiv 3 \pmod 8; \\ 4x, & \text{if}\: x \equiv 4 \pmod 8; \\ (4x+1)-1, & \text{if}\: x \equiv 5 \pmod 8; \\ (4x-2)+2, & \text{if}\: x \equiv 6 \pmod 8; \\ (4x-1)+1, & \text{if}\: x \equiv 7 \pmod 8 \\ \end{cases} \\ &= 4x. \quad\qed \end{aligned} $$

Claim 4: 任意の $x$ に対し、$( ( (x \oplus x)\oplus x)\oplus x)\oplus x = 5 \otimes x$ が成り立つ。

Proof

$$ \begin{aligned} ( ( (x \oplus x)\oplus x)\oplus x)\oplus x &= 4x \oplus x \\ &= \roundcirc{4x+x} \\ &= \roundcirc{5x} \\ &= 5\otimes x. \quad\qed \end{aligned} $$

Claim 5: 任意の $x$ に対し、$( ( ( (x \oplus x)\oplus x)\oplus x)\oplus x)\oplus x = 6 \otimes x$ が成り立つ。

Disproof

$x = 2^{52}+1$ とする。下記が成り立つ。 $$ \begin{aligned} (5\otimes x)\oplus x &= (5x-1)\oplus x \\ &= 6x-2, \\ 6\otimes x &= 6x+2. \quad\qed \end{aligned} $$

おまけ

note: 非正規化数に関しては、frexp を実装するときと同様にして、上記のケースに帰着できるはずです。非正規化数は実質的には固定小数点数なので、なんかいい感じになるはずです。

note: オーバーフローに関しても、なんやかんやでうまくいくはずです。

note: x == -0.0 のケースでも問題ないです。

おわり

「反例がない気はしないんだが」（笑）

*1:“floating-point number” は、無限大は含みますが NaN は含まない用語として定義されていますね。

ABC 400 B の浮動小数点型解法の正当性の証明

ずんだもんなのだ
今日はこの解法の正当性の証明をやっていくのだ pic.twitter.com/UeXkq4BG5O

— えびちゃん🍑🍝🦃 (@rsk0315_h4x) 2025å¹´4月5æ—¥

fn solve_f64(n: f64, m: i32) -> Option<u32> {
    if n == 1.0 || m == 1 {
        (n + (m as f64) <= 1.0e9).then(|| n as u32 + m as u32)
    } else {
        let res = (n.powi(m + 1) - 1.0) / (n - 1.0);
        (res <= 1.0e9).then(|| res as u32)
    }
}

20 分くらいかけて全部のケースを回して正しい答えを出力することがわかったのだ
今回も簡単だったのだ
見てくれてありがとうなのだ

— えびちゃん🍑🍝🦃 (@rsk0315_h4x) 2025å¹´4月5æ—¥

バ、
馬鹿もーんっ・・・・！
なんだっ・・！ この証明はっ・・・・・・・・！
通るかっ・・・・・・！
こんなもん・・！*1

まぁ通るんですが、つまらないかなということで、やっていきます。

記事化を待つ

— ねぼこ (@nebocco27) 2025å¹´4月6æ—¥

問題リンク：ABC 400 B

証明

下記が成り立つことは既知とします。 $$ \sum_{i=0}^M N^i = \begin{cases} M+1, & \text{if}\: N=1; \\ \displaystyle\frac{N^{M+1}-1}{N-1}, & \text{if}\: N\ne 1. \end{cases} $$

powi について

'llvm.powi.*' Intrinsic には下記の記述があります。

Semantics: This function returns the first value raised to the second power with an unspecified sequence of rounding operations.

たとえば、$\texttt{powi}(n, 7)$ は $( (n\otimes n)\otimes (n\otimes n) ) \otimes ( (n\otimes n)\otimes n)$ かもしれないし、$( ( (n\otimes n)\otimes n)\otimes ( (n\otimes n)\otimes n) )\otimes n$ かもしれないし、$( ( ( ( ( (n\otimes n)\otimes n)\otimes n)\otimes n)\otimes n)\otimes n)\otimes n$ かもしれないし、みたいな感じだと思います。

コード (src) は繰り返し二乗法のようですが、累乗のアルゴリズムはそれだけではないですし、それに限定した仕様にする理由はないのでこういう感じになっているのでしょうかね。$\roundcirc{n^{3.5}}\otimes\roundcirc{n^{3.5}}$ とか $\roundcirc{n^8}\oslash n$ のようにされると困りますが、そういうことはないと仮定しておきます。

上記の前提であれば、非負整数 $n\le 2^{53}$ と $i$ に対して、$n^i \le 2^{53}$ であれば $\texttt{powi}(n, i) = n^i$ が成り立ちます。 また、$n^i\gt 2^{53}$ の場合でも、$\texttt{powi}(n, i)\lt\infty$ ならば、ある実数 $|\theta|\le (1+2^{-53})^{i-1}-1$ に対して $\texttt{powi}(n, i) = (1+\theta)n^i$ が成り立ちます。

場合分け

整数 $1\le N\le 10^9$, $1\le M\le 100$ について考えます。

Claim 1: $N = 1 \wedge M+1 \le 10^9 \implies M\oplus 1 = M+1$ が成り立つ。

Proof: 明らか。$\qed$

Claim 2: $N = 1 \wedge M+1 \gt 10^9 \implies M\oplus 1 \gt 10^9$ が成り立つ。

Proof: 明らか。$\qed$

Claim 3: $N \gt 1 \wedge M = 1 \wedge N+1 \le 10^9 \implies N\oplus 1 = N+1$ が成り立つ。

Proof: 明らか。$\qed$

Claim 4: $N \gt 1 \wedge M = 1 \wedge N+1 \gt 10^9 \implies N\oplus 1\gt 10^9$ が成り立つ。

Proof: 明らか。$\qed$

$f(N, M) = (\texttt{powi}(N, M+1)\ominus 1)\oslash(N\ominus 1)$ とします。

Lemma 5: $M\gt 1 \implies \frac{\log_2(10^9)}M \le \frac{53}{M+1}$ が成り立つ。

Proof

$$ \begin{aligned} \frac{\log_2(10^9)}M \le \frac{53}{M+1} &\iff 1+\frac1M \le \frac{53}{\log_2(10^9)} \\ &\iff \frac1M \le \frac{53-\log_2(10^9)}{\log_2(10^9)} \\ &\iff M\ge \frac{\log_2(10^9)}{53-\log_2(10^9)} \gt 1.294. \end{aligned} $$

$M$ の整数性より従う。$\qed$

Claim 6: $N \gt 1 \wedge M\gt 1 \wedge \frac{N^{M+1}-1}{N-1}\le 10^9 \implies f(N, M) = \frac{N^{M+1}-1}{N-1}$ が成り立つ。

Proof

$N^{M+1}\le 2^{53}$ について考える。 $$ \begin{aligned} N^{M+1}\le 2^{53} % &\iff (M+1)\log_2(N) \le 53 \\ % &\iff M \le \frac{53}{\log_2(N)}-1 \\ % &\iff M \le \frac{53-\log_2(N)}{\log_2(N)}. &\iff N \le 2^{\frac{53}{M+1}}. \end{aligned} $$ また、$\frac{N^{M+1}-1}{N-1} \le 10^9$ について考える。 $$ \begin{aligned} \frac{N^{M+1}-1}{N-1} \le 10^9 &\implies N^M \le 10^9 \\ % &\iff M\log_2(N) \le 9\log_2(10) \\ % &\iff M \le \frac{9\log_2(10)}{\log_2(N)}. &\iff N \le 2^{\frac{\log_2(10^9)}M}. \end{aligned} $$ Lemma 5 より、$\frac{N^{M+1}-1}{N-1} \le 10^9 \implies N^{M+1}\le 2^{53}$ が成り立つ。

よって、$f(N, M) = (N^{M+1}\ominus)\oslash(N\oslash 1) = \frac{N^{M+1}-1}{N-1}$ となる。$\qed$

Lemma 7: $N \gt 1 \wedge M\gt 1 \wedge N^{M+1}\gt 2^{53} \implies \frac{N^{M+1}-1}{N-1}\gt 2^{53\cdot 2/3}$ が成り立つ。

Proof: Lemma 5, Claim 6 と同様にして示せる。$\qed$

Claim 8: $N \gt 1 \wedge M\gt 1 \wedge \frac{N^{M+1}-1}{N-1}\gt 10^9 \implies f(N, M) \gt 10^9$ が成り立つ。

Proof

Case 1: $N^{M+1}\le 2^{53}$

Claim 6 同様に $f(N, M) = \frac{N^{M+1}-1}{N-1}$ が成り立ち、$f(N, M)\gt 10^9$ となる。

Case 2: $N^{M+1}\gt 2^{53}$

$\texttt{powi}(N, M+1)=\infty$ のとき $f(N, M)=\infty\gt 10^9$ となる。

そうでないとき、実数 $|\theta_{\otimes}|\le (1+2^{-53})^M-1$, $|\theta_{\ominus}|\le 2^{-53}$, $|\theta_{\oslash}|\le 2^{-53}$ が存在し、下記が成り立つ。 $$ \begin{aligned} f(N, M+1) &= (\texttt{powi}(N, M+1)\ominus 1)\oslash (N\oslash 1) \\ &= ( (1+\theta_{\otimes})N^{M+1}\ominus 1)\oslash (N-1) \\ &= (1+\theta_{\ominus})( (1+\theta_{\otimes})N^{M+1}-1)\oslash (N-1) \\ &= (1+\theta_{\oslash})\cdot \frac{(1+\theta_{\ominus})( (1+\theta_{\otimes})N^{M+1}-1)}{N-1} \\ &= (1+\theta_{\oslash})(1+\theta_{\ominus})\cdot \frac{(1+\theta_{\otimes})N^{M+1}-1}{N-1} \\ &= (1+\theta_{\oslash})(1+\theta_{\ominus})\cdot \frac{(1+\theta_{\otimes})(N^{M+1}-1)+\theta_{\otimes}}{N-1} \\ &= (1+\theta_{\oslash})(1+\theta_{\ominus})(1+\theta_{\otimes})\cdot \frac{N^{M+1}-1}{N-1} \\ &\qquad + (1+\theta_{\oslash})(1+\theta_{\ominus})\cdot\frac{\theta_{\otimes}}{N-1} \\ &\gt (1-2^{-53})^{M+2}\cdot 2^{53\cdot 2/3} + (1+2^{-53})^2\cdot\frac{-( (1+2^{-53})^M -1)}{N-1} \\ &\ge (1-2^{-53})^{102}\cdot 2^{53\cdot 2/3} - (1+2^{-53})^2\cdot\frac{(1+2^{-53})^{100}-1}{2-1} \\ &\gt 4.329\times 10^{10} \\ &\gt 10^9. \quad\qed \end{aligned} $$

よって、各 $N$, $M$ について、$\sum_{i=0}^M N^i$ が $10^9$ 以下ならば $f(N, M) = \sum_{i=0}^M N^i$ であり、$10^9$ より大きいならば $f(N, M)\gt 10^9$ であることが示せました。

おまけ

$f(N, 1)$ がうまくいかない例について考えます。なお、結合の仕方は一通りしかないため $\texttt{powi}(N, 2) = N\otimes N$ となります。

たとえば $N=189812534$ のとき、下記のようになります。 $$ \begin{aligned} f(189812534, 1) &= ( (189812534\otimes 189812534)\ominus 1)\oslash(189812534\ominus 1) \\ &= ( (189812534^2-4)\ominus 1)\oslash(189812534\ominus 1) \\ &= (189812534^2-4)\oslash(189812534-1) \\ &= {\small 189812534.999999970197677}{\footnotesize 6123046875}. \end{aligned} $$ よって、これを切り捨てて整数にすると $1+N$ とは異なる値が得られます。

他にも、$189812538$ や $999999992$ など、たくさんの $N$ がこのようになります。

おわり

おわりです。全探索の方が早く済みました。いかがでしたか？

*1:出典：https://youtu.be/YPaJ30gY1ic?t=61。

correct rounding への道 (3) Newton’s method

今回は整数を使わずに、いわゆる「数値計算」っぽいアルゴリズムで平方根を求めていきます。そういうアルゴリズムであっても 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)$ のグラフ

$\hat p(a^{[2]}, 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}$ の倍数であるため。