えびちゃんの日記

えびちゃん（競プロ）の日記です。

(1.0 / 49.0) * 49.0 < 1.0 について

浮動小数点数

仮数部の精度が 53 bits の浮動小数点数型で計算すると下記のようになります。

>>> 1 / 49 * 49
0.9999999999999999
>>> f'{1/49*49:.100g}'
'0.99999999999999988897769753748434595763683319091796875'

今回はこれについて掘り下げます。

記法・前提知識

その文脈における浮動小数点数型と丸め方向において、実数 $x$ を丸めた値を $\roundp{x}$ と書くことにします。 型としては double（IEEE 754 で言うところの binary64）を用いて、丸め方向に roundTiesToEven を採用している言語・処理系が多いので、ここではそれを考えます。

知らない人向け：表せるもののうちで一番近いものに丸めるモード。この記事の内容では、ちょうど中間にあるケースは存在しないので、そういうことは気にしなくて大丈夫。詳しくは別の記事で解説します。

また、浮動小数点数 $x$, $y$ に対して、$x\otimes y \triangleq \roundp{x\times y}$ および $x\oslash y\triangleq \roundp{x\div y}$ とします。

たとえば、上記の例より $$ \begin{aligned} \roundp{\tfrac1{49}}\otimes 49 &= {\footnotesize 0.999999999999999888}{\scriptsize 97769753748434595763}{\tiny 683319091796875} \\ \end{aligned} $$ です。また、簡単な例として $\roundp1=1$ もあります。

$\gdef\hlog#1{\log_2\hfloor{#1}}$

正の実数 $x$ に対し、整数 $i$ を用いて $2^i$ と書ける有理数のうち、$x$ 以下で最大のものを $\hfloor x$ で表します。 たとえば $\hfloor{3.5} = 2$、$\hfloor1 = 1$、$\hfloor{0.4} = 0.25$ などです。 下記では「正の実数 $x$ に対し、$2^e \le x \lt 2^{e+1}$ なる唯一の整数 $e$ を考えたい」という局面が多く、これを $\hlog{x}$ と書けて便利です。

$\hlog{x} = i$ のとき、$2^{i+1-p}$ の位までが保持され、$2^{i-p}$ 以下の桁は丸められます。 double においては $p=53$ です。

解説

まず、1 / 49 の部分を考えます。

$$ \begin{aligned} \roundp{\tfrac1{49}} &= {\footnotesize 0.020408163265306}{\scriptsize 12082046367561360966}{\tiny 64704382419586181640625} \\ &= \texttt{14e5e0a72f0539}_{(16)} \times 2^{-58} \\ &= \tfrac1{49}\cdot(1-23\cdot 2^{-58}) \end{aligned} $$ となります。10 進の小数表記・16 進の指数表記と、誤差に着目した形式を挙げました。 ここでは一番下の表記で考えるのがわかりやすいと思いますが、いくつかの表現方法があることを覚えておき、状況に応じて適切なものを選べるとよいでしょう。 数学寄りの文脈で解析をするときは一番下の形が処理しやすい気がします。

さて、$\roundp{\tfrac1{49}}\otimes49$ を考えます。 $$ \roundp{\tfrac1{49}}\otimes 49 = \roundp{\roundp{\tfrac1{49}}\times 49} = \roundp{1-23\cdot 2^{-58}}. $$ $\tfrac12 \lt 1-23\cdot 2^{-58} \lt 1$ より、$\hlog{\roundp{\tfrac1{49}}\otimes 49}=-1$ とわかります。すなわち、$2^{-53}$ の位まで丸められます。

誤差に関して、$23\cdot 2^{-58} = \tfrac{23}{16}\cdot 2^{-54}$ より、$1\cdot 2^{-54} \lt 23\cdot 2^{-58} \lt 2\cdot 2^{-54}$ とわかります。 よって、$\roundp{\tfrac1{49}}\times49$ は $1$ よりも $1-2^{-53}$ に近く、$1-2^{-53}$ に丸められます。

疑問・興味

ここで、いくつかの疑問が自然に出てくると思います。

  • $\roundp{\tfrac1x}\otimes x \gt 1$ となるような $x$ はある？
    • ただし $\infty$ になるケースは除く
  • $x\in[1\ldots 2)$ なる $x$ のうちで $\roundp{\tfrac1x}\otimes x\lt 1$ となる最小のものは？
    • すなわち $\hlog x=0$ の範囲
    • $49$ の例から $x = \tfrac{49}{32}$ でも成り立つことはわかるので、存在性はわかる

ネタバレ

  • $\texttt{10000004e62386}_{(16)}\times 2^{970}$ など
    • $\roundp{\tfrac1x}$ が正規化数になる範囲では存在しない
  • $\texttt{1000000f5cbf2a}_{(16)}\times 2^{-52}$

これについて考察していきます。

疑問 1 の解答

$\gdef\emin{e_{\text{min}}}$ $\gdef\emax{e_{\text{max}}}$

まず、正規化数の最小値を $1\times 2^{\emin}$、正規化数の最大値を $(2-2^{1-p})\times 2^{\emax}$ とおきます。 double においては、$1\times 2^{-1022}$ および $(2-2^{-52})\times 2^{1023}$ です。 一般に、$|\emin|+1 = \emax$ が成り立ちます。

証明などは折りたたんでいます。

Claim 1: $x\in(1\ldots 2)$ のとき、$\roundp{\tfrac1x}\otimes x\in \{1-2^{-p}, 1\}$ が成り立つ。

Proof

このとき $\tfrac1x\in(\tfrac12\ldots1)$ であり、$\hlog{\tfrac1x}=-1$ である。 そのため、ある $|\varepsilon|\le 2^{-1-p}$ を用いて $$ \roundp{\tfrac1x} = \tfrac1x + \varepsilon $$ と書ける。よって、 $$ \begin{aligned} \roundp{\tfrac1x}\times x &= (\tfrac1x + \varepsilon)\times x \\ &= 1 + \varepsilon x \end{aligned} $$ となる。$\varepsilon$ と $x$ の範囲から、$|\varepsilon x|\lt 2^{-p}$ とわかる。

$1$ と隣り合う二数は $1-2^{-p}$ と $1+2^{1-p}$ であり、$1$ に丸められる範囲は $$[1-2^{-1-p}\ldots 1+2^{-p}]$$ であるから、$1+\varepsilon x\lt 1+2^{-p}$ より $\roundp{\tfrac1x}\otimes x\le 1$ となる。 同様にして $\roundp{\tfrac1x}\otimes x\ge 1-2^{-p}$ もわかる。$\qed$

Lemma 2: 浮動小数点数 $x\in(1\ldots2)$ に対し、$\roundp{\tfrac1x}\in(\tfrac12\ldots1)$ が成り立つ。

Proof

$1\oslash(1+2^{1-p})\lt 1$ および $1\oslash(2-2^{1-p})\gt\tfrac12$ を示せば十分。

$$ \begin{aligned} 1\div(1+2^{1-p}) &= 1 - 2^{1-p} + 2^{2(1-p)}\cdot (1+2^{1-p})^{-1} \\ &= 1 - 2^{-p} - 2^{-p} + 2^{-p}\cdot 2^{2-p}\cdot (1+2^{1-p})^{-1} \\ &= (1-2^{-p}) - 2^{-p}\cdot (1 - 2^{2-p}\cdot (1+2^{1-p})^{-1}) \\ &= (1-2^{-p}) - 2^{-p}\cdot (1 - 2\cdot (2^{p-1}+1)^{-1}) \\ &\lt 1-2^{-p}; \\ 1\div(2-2^{1-p}) % &= \tfrac12 \div (1-2^{-p}) \\ &= \tfrac12 \cdot (1+2^{-p}+2^{-2p}\cdot(1-2^{-p})^{-1}) \\ &\gt \tfrac12 \cdot (1+2^{-p}) \end{aligned} $$ より従う。$\qed$

Claim 3: $e \in[\emin\ldots\emax-2]$ に対して、$x\in(2^e\ldots2^{e+1})$ のとき、$\roundp{\tfrac1x}\otimes x\in \{1-2^{-p}, 1\}$ が成り立つ。

Proof

ある浮動小数点数 $\alpha\in(1\ldots2)$ を用いて $x = \alpha\times 2^e$ と書ける。 また、Lemma 2 より $\roundp{\tfrac1\alpha}\in(\tfrac12\ldots1)$ が成り立つことに注意しつつ $$ \begin{aligned} \roundp{\tfrac1x} &= \roundp{\tfrac1\alpha\times 2^{-e}} \\ &= \roundp{\tfrac1\alpha}\times 2^{-e} \\ &= \roundp{\tfrac2\alpha}\times 2^{-1-e} \end{aligned} $$ が成り立つ。 ここで $\emin\le -1-e\le \emax$ が成り立っていることに注意せよ。さて、 $$ \begin{aligned} \roundp{\tfrac1x}\times x &= (\roundp{\tfrac1\alpha}\times 2^{-e})\times (\alpha\times 2^e) \\ &= \roundp{\tfrac1\alpha}\times \alpha \end{aligned} $$ が成り立つので、Claim 1 から従う。$\qed$

Conjecture 4: $x\in(2^{\emax-1}\ldots 2^{\emax+1})$ のとき、$\roundp{\tfrac1x}\otimes x\in\{1-2^{-p}, 1\}$ が成り立つ。

Disproof

double において、$x = \texttt{10000004e62386}_{(16)}\times 2^{970}$ とすると $$\roundp{\tfrac1x}\otimes x = 1+2^{-52}\gt 1$$ が成り立つ。

$\tfrac1x\lt 2^{-(\emax-1)} = 2^{\emin}$ より、正規化数の最小値以下に丸められるので、ある $|\varepsilon|\le 2^{\emin-p}$ を用いて $$ \roundp{\tfrac1x} = \tfrac1x + \varepsilon $$ と書ける。よって、 $$ \begin{aligned} \roundp{\tfrac1x}\times x &= (\tfrac1x+\varepsilon)\times x \\ &= 1 + \varepsilon x \end{aligned} $$ となる。$|\varepsilon x| \lt 2^{\emax+1+\emin-p} = 2^{2-p}$ でしか抑えられず、適切な $x$ を選ぶことで反例を作れそうである。

Claim 5: $e\in[-(\emax+1)\ldots\emin)$ に対して、$x\in(2^e\ldots 2^{e+1})$ のとき、$\roundp{\tfrac1x}\otimes x\in\{1-2^{-p}, 1\}$ が成り立つ。

Proof

このとき $\tfrac1x\in(2^{-(e+1)}\ldots2^{-e})$ であり、$\hlog{\tfrac1x}=-(e+1)$ である。 そのため、ある $|\varepsilon|\le 2^{-(e+1)-p}$ を用いて $$ \roundp{\tfrac1x} = \tfrac1x + \varepsilon $$ と書ける。よって、

$$ \begin{aligned} \roundp{\tfrac1x}\times x &= (\tfrac1x + \varepsilon)\times x \\ &= 1 + \varepsilon x \end{aligned} $$ となる。$\varepsilon$ と $x$ の範囲から、$|\varepsilon x|\lt 2^{-p}$ とわかる。 Claim 1 同様にして、成り立つことがわかる。$\qed$

note: Claim 3 では $x$ が正規化数の範囲をまとめて扱ったが、Claim 5 の $e$ の範囲を広げることでも示せそうである。

Claim 6: $e\in[\emin\ldots\emax]$ に対して、$x=1\times2^e$ のとき $\roundp{\tfrac1x}\otimes x=1$ が成り立つ。

Proof

$\roundp{\tfrac1x}=2^{-e}$ より従う。$\qed$

Claim 7: $\roundp{\tfrac10}\otimes0$ および $\roundp{\tfrac1{\infty}}\otimes\infty$ は NaN となる。また、$x\in(0\ldots 2^{-(\emax+2)}]$ では $\roundp{\tfrac1x}\otimes x = \infty$ となる。

Proof

前半は、$\roundp{\tfrac10}=\infty$ および $\roundp{\tfrac1{\infty}}=0$ であることと $0\otimes \infty$ が NaN であることから従う。 後半は、該当の範囲で $\roundp{\tfrac1x}=\infty$ であることと $x\gt 0$ であることから従う。 $\qed$

よって、$\roundp{\tfrac1x}$ が正規化数となる範囲では $\roundp{\tfrac1x}\otimes x\in\{1-2^{-p}, 1\}$ であることがわかりましたが、そうでない範囲においては反例があることがわかりました。

たとえば、$x = \texttt{10000004e62386}_{(16)}\times 2^{970}$ に対して $\roundp{\tfrac1x}\otimes x = 1 + 2^{1-p}$ となります。 また、最小性は示していませんが、$x = \texttt{18000003000000}_{(16)}\times 2^{970}$ に対して $\roundp{\tfrac1x}\otimes x = 1 - 2\cdot 2^{-p}$ となります。

疑問 2 の解答

$\roundp{\tfrac1x}\otimes x \lt 1$ なる $x\in(1\ldots2)$ の最小値について考えます。

Claim 8: $\varepsilon^2\lt 2^{-(p+1)}$ に対して $x=1+\varepsilon$ と書けるとき、$\roundp{\tfrac1x}\otimes x=1$ が成り立つ。

Proof

等比級数の和の公式から、 $$ \begin{aligned} \tfrac1x &= \tfrac1{1+\varepsilon} \\ &= 1 - \varepsilon + \tfrac{\varepsilon^2}{1+\varepsilon} \end{aligned} $$ である。 ここで、$\tfrac{\varepsilon^2}{1+\varepsilon} \le \varepsilon^2 \lt 2^{-(p+1)}$ より、 $1\oslash x = 1-\varepsilon$ となる*1。 この範囲では $$ \begin{aligned} \roundp{\tfrac1x}\times x &= (1-\varepsilon)\times(1+\varepsilon) \\ &= 1-\varepsilon^2 \end{aligned} $$ となり、$\varepsilon^2\lt 2^{-(p+1)}$ より $\roundp{\tfrac1x}\otimes x=1$ となる。 $\qed$

Claim 8 の証明から、$\tfrac{\varepsilon^2}{1+\varepsilon}$ による誤差が十分に大きくなれば $\roundp{\tfrac1x}\otimes x\lt 1$ が成り立ちそうです。 $1+2^{-(p+1)/2}$ 付近に解があることに期待して、仮数部を昇順に全探索してもよいかもしれません。

実際、仮数部が 53 bits である double や、long double のうち仮数部が 64 bits であるものについては比較的すぐ見つかりました。一方、仮数部が 113 bits である __float128（あるいはそういう処理系の long double）ではうまくいきませんでした。そこで、もう少し考察を進めます。

まず、整数 $m\in(2^{p-1}\ldots 2^p)$ に対して $x = m\times 2^{1-p}$ と表すことができ、 $$ \begin{aligned} \tfrac1x &= \tfrac{2^p}x\times 2^{-p} \\ &= \tfrac{2^{2p-1}}m\times 2^{-p} \end{aligned} $$ です。ここで、$m$ の範囲から $\tfrac{2^{2p-1}}m\in(2^{p-1}\ldots 2^p)$ であることがわかります。 よって、 $$ \roundp{\tfrac1x}\in\left\{\floor{\tfrac{2^{2p-1}}m}\times 2^{-p}, \ceil{\tfrac{2^{2p-1}}m}\times 2^{-p}\right\} $$ となります。 また、 $$ \Floor{\frac{2^{2p-1}}m} = \frac{2^{2p-1} - (2^{2p-1}\bmod m)}m $$ です。

丸めの際は近い方が採用されるため、$(2^{2p-1}\bmod m) \lt \tfrac m2$ のときは $\floor{\ldots}\times 2^{-p}$ の方が採用され、$(2^{2p-1}\bmod m)\gt \tfrac m2$ のときは $\ceil{\ldots}\times 2^{-p}$ の方が採用されます。$m$ の範囲から $\tfrac{2^{2p-1}}m$ が（2 進法で）無限小数となることがわかり、$(2^{2p-1}\bmod m)=\tfrac m2$ となるケースは存在しないことが従います。

よって、次のようになります。

Claim 9: $r = 2^{2p-1}\bmod m$ が $r\gt\tfrac m2$ を満たすとき、$\roundp{\tfrac1x}\otimes x=1$ が成り立つ。

Proof

該当の条件のとき、 $$ \begin{aligned} \roundp{\tfrac1x}\times x &= \ceil{\tfrac{2^p}x}\times 2^{-p} \\ &\gt \tfrac{2^p}x\times 2^{-p} \\ &= \tfrac1x \end{aligned} $$ であることと、Claim 1 より従う。$\qed$

すなわち、$\floor{\ldots}\times 2^{-p}$ が採用される方のみ考えればよいです。

Claim 10: $r = 2^{2p-1}\bmod m$ が $2^{p-2}\lt r\lt \tfrac m2$ を満たすとき、$\roundp{\tfrac1x}\otimes x=1-2^{-p}$ となる。

Proof

$r\lt \tfrac m2$ の範囲で考える。 $$ \begin{aligned} \roundp{\tfrac1x}\times x &= \left(\floor{\tfrac{2^{2p-1}}m}\times 2^{-p}\right)\cdot (m\times 2^{1-p}) \\ &= \left(\tfrac{2^{2p-1}-r}m\times 2^{-p}\right)\cdot (m\times 2^{1-p}) \\ % &= (2^{2p-1}-r)\cdot 2^{-(2p-1)} \\ &= 1 - r\cdot 2^{-(2p-1)} \end{aligned} $$ より、$r\cdot 2^{-(2p-1)} \gt 2^{-(p+1)}$ であれば $1-2^{-p}$ に丸められる。 すなわち、$r\gt 2^{p-2}$ であることが必要十分である。$\qed$

Observation 11: $r = 2^{2p-1}\bmod m$ は次のようになっている。

$p=53$ とする。また、$m = 2^{p-1}+i$ とする。

$i$ $0$ $1$ $2$ $3$ $\cdots$ $47453132$ $47453133$ $47453134$ $47453135$
$r$ $0$ $2$ $8$ $18$ $\cdots$ ${\tiny 4503599473218848}$ ${\tiny 4503599663031378}$ ${\scriptsize 178020282}$ ${\scriptsize 367832819}$

ある $k$ に対し、$r = 2i^2 - k\cdot(2^{p-1}+i)$ となっている。

Proof

$$ \begin{aligned} 2^{2p-1} &= (2^p-2i)\times(2^{p-1}+i) + 2i^2 \\ &= (2^p-2i)\times(2^{p-1}+i) + \Floor{\tfrac{2i^2}{2^{p-1}+i}}\cdot (2^{p-1}+i) + (2i^2\bmod (2^{p-1}+i)) \\ &= \left(2^p-2i+\Floor{\tfrac{2i^2}{2^{p-1}+i}}\right)\times(2^{p-1}+i) + \left(2i^2-\Floor{\tfrac{2i^2}{2^{p-1}+i}}\cdot(2^{p-1}+i)\right) \end{aligned} $$ である。$k$ は $2i^2\div(2^{p-1}+i)$ の商に他ならない。$\qed$

そこで、次のような方針が考えられます。

  1. $k$ を固定する。
  2. その範囲で $r\gt 2^{p-2}$ を満たす $i$ の最小値 $i^{\star}$ を求める。
    • ここは二分探索できる。
  3. $i^{\star}$ で計算したあまり $r^{\star}$ が $r^{\star}\lt\tfrac m2$ を満たしていればそれが答え。
    • 満たしていなければ次の $k$ を試す。

$k$ の上界はよくわからないものの、いくつか試したところでは十分小さい $k$ で停止しました。今回の目的は計算量改善ではないのでよしとします。

Example 12:

メジャーな処理系における _Float16, float, double, long double, __float128 を想定し、上記の計算をした。

type minimum $x\in(1\ldots2)$ subject to $\roundp{\tfrac1x}\otimes x\lt 1$ $p$ $k$
_Float16 $\texttt{41c}_{(16)}\times2^{-10}$ $11$ $1$
float $\texttt{801467}_{(16)}\times2^{-23}$ $24$ $6$
double $\texttt{1000000f5cbf2a}_{(16)}\times2^{-52}$ $53$ $29$
long double $\texttt{800000005a82799a}_{(16)}\times2^{-63}$ $64$ $0$
__float128 $\texttt{1000000000000022df0668c89bf13}_{(16)}\times2^{-112}$ $113$ $9$

long double 以下の型については、$1$ から昇順に全探索したものと値が一致することを確認した。 また、各型の値について、実際に $\roundp{\tfrac1x}\otimes x\lt 1$ となることを確認した。

あとがき

「浮動小数点数型」に対して、「謎の誤差が出て変な感じになる」というような感覚を持っている人はおそらく多いのではないかと思っていますが、こうした解析で遊んでいるとある程度「数学側のやり方で（手間は多少かかっても）ある程度まともに扱える」という気持ちになってくるのではないでしょうか。少なくともえびちゃんの中ではそういう感じになってきました。コンテスト中になにかできそうな気はまだしません。

とはいえまだ基礎パートの解説記事を書いていないので、取っかかりが掴めていない段階の人にはまだかもという気はしています。 もう少し待っていてもらえるとうれしいです。

というのも、ある程度自分が親しんでいない段階では基礎用の記事になにを書くべきか判断しかねたからです。 ボツになった下書きが二つか三つくらいあります。

その練習のために書いたのが、Sqrt Inequality の記事 と floor(n * 0.1) の記事 と今回の記事でした（最初にやるにしては重くないか？）。 もう一つ練習に書こうとしているのがあるので、それを終えたら基礎用の記事を書こうと思っています。

こういう細々とした詰め詰めな解析も楽しいんですが、競プロ文脈で（コンテスト中に計算量のオーダーを見積もるくらいの手軽さで）誤差の見積もりができるようにもなれたらうれしいので、そういうところのお勉強もしていきたいです。

あまり多くの競プロ er が興味を持たないところに興味を持ってしまうところ、えびちゃんのよいところである気はしていますが、自分としてはさみしいところでもあります。

おまけ

#include <cstdio>
#include <string>
#include <quadmath.h>

using f16 = _Float16;
using f32 = float;
using f64 = double;
using f80 = long double;
using f128 = __float128;

std::string q(f128 x) {
    char buf[1000];
    quadmath_snprintf(buf, sizeof buf, "%.999Qg", x);
    return std::string(buf);
}

int main() {
    int x = 27355;
    printf("%.99g\n", f16(f16(f16(1.0) / f16(x)) * f16(x)));
    printf("%.99g\n", f32(f32(1.0F / f32(x)) * f32(x)));
    printf("%.99g\n", 1.0 / f64(x) * f64(x));
    printf("%.99Lg\n", 1.0L / f80(x) * f80(x));
    printf("%s\n", q(f128(1.0) / f128(x) * f128(x)).data());
}

0.99951171875
0.999999940395355224609375
0.99999999999999988897769753748434595763683319091796875
0.9999999999999999999457898913757247782996273599565029144287109375
0.99999999999999999999999999999999990370350278063820734720110287075363407309491758923059023800306022167205810546875

wandbox.org

5 つの型すべてで $\roundp{\tfrac1x}\otimes x\lt 1$ にさせる最小の非負整数です。実際には _Float16 では $\roundp{x}\ne x$ なのでちょっとずるいです。

自分で仮数部 $p$ bits の浮動小数点数型をエミュレートする型を作り、これら以外の型について遊ぶこともできます（遊びました）。

おわり

おわりです。

*1:$1+\varepsilon\lt 2$ が該当の浮動小数点数型で表せていることから、$\roundp{1-\varepsilon} = 1-\varepsilon$ です。