えびちゃんの日記

えびちゃん(競プロ)の日記です。

浮動小数点型の算術とお近づきになりたい人向けの記事

お近づきになりたい人向けシリーズです。

いろいろなトピックを詰め込みましたが、「これら全部を知らないといけない」のようなつもりではなく、いろいろなことを知るきっかけになったらいいなという気持ちなので、あまり身構えずにちょっとずつ読んでもらえたらうれしい気がします。

まえがき

浮動小数点型といえば、競プロ er がよくわからずに使っているものの代表格でございましょう。 初心者のうちは素朴に「実数を扱うときはこれを使う」とだけ思っていて、そのうち誤差の概念を知って「誤差があるらしいが祈れば大丈夫」と思う人がたくさんいます。中級者になってくると「浮動小数点型を使うのはよくなくて、基本的には整数型で処理するべき」のような考え方になっていきがちです。誤差のことを知っている中級者には、「浮動小数点型では、なにをやっても誤差でめちゃくちゃになる」のように(初心者とは対極の)悲観的な意識を持っている人も多そうです。

いずれにせよ、これによる計算でどの程度の誤差が出るとか、こういうケースでは誤差が出ないとか、このように補正すれば意図通りの計算がある程度の精度でできるとか、そうしたことを論ずることができる人はかなり限られているのが現状なのではないでしょうか。 コンテスト中にこうした評価を強いられることは(昨今の傾向では)あまりないですが、できるようになっても損はないでしょう。

証明や詳細な部分に関しては、適宜折りたたんで記載しています。スタイルの都合上、折りたたみの終端がわかりにくいので、$\eod$ の記号によって明示しています。

折りたたみの例

ここに詳細を書く。$\eod$

ここは折りたたみの外です。

予備知識

規格

浮動小数点の算術は、IEEEアイトリプルイー 754 (ISO/IEC 60559) という規格で定められています[Std754]。最新のバージョンは IEEE 754-2019 (ISO/IEC 60559:2020) です。

よく見る説明の中に「この型は 53-bit の仮数部を持つことが IEEE 754 の規格で定められている」のようなものがあったりして、「IEEE 754 というのは型のフォーマットを定めている規格なのだな」と思う人もいそうですが、実際にはそれだけを定めている規格ではありません。 規格のタイトルが “IEEE Standard for Floating-Point Arithmetic” であるように、浮動小数点の算術(各演算の結果はどうとか、特殊なケースではこういう振る舞いをするとか、言語側で提供するべき関数とか)についての規格です。

よく見る浮動小数点型は 2 進法ですが、規格では 10 進法についても触れられています。簡単のため、ここでは 2 進法の方のみ触れることにします。 メジャーな言語(の多くの処理系)では、この規格中で binary64 と呼ばれている型を使うことができます。実際には、型の表現が binary64 であることしか定めていない言語(算術までは厳密に従っていないかも?)や、言語仕様としては binary64 であることすら決めていない言語もありますが、その調査についてはこの記事の範囲外とします。基本的な対応は次の表の通りです。

言語 型名
C double (ref, see Annex F)
C++ std::float64_t (ref)
Rust f64 (ref)
Python float (ref)
Kotlin Double (ref)
JavaScript number (ref)

用語

規格中では、floating-point datum, floating-point number, floating-point operation, floating-point representation の用語が定義されています。

  • floating-point datum
    • 該当の形式で表現可能な数 (floating-point number) または NaN。
    • 規格中では、表現やエンコード方式とは必ずしも区別されずに用いられる。
  • floating-point number
    • NaN ではない floating-point datum。
    • 該当の形式で表現可能な、有限または無限の数。
      • 実数の部分集合と $\{-\infty, +\infty\}$ の和集合。
  • floating-point operation
    • 引数または返り値が floating-point datum であるような演算。
  • floating-point representation
    • 該当の形式における(エンコードされていない)要素。
    • 有限値、(符号つきの)無限大、NaN のいずれかを表す。
      • 有限値は、符号部・指数部・仮数部の 3 つからなる。

floating-point number が NaN を含まない用語であることから、double などの総称に「浮動小数 型 (floating-point number types)」を用いるのは不適切かもしれません? これからは浮動小数点型と呼んでいくことにします*1

NaN (Not-a-Number) は、不正な演算の返り値などを表すのに使われる形式的な元です。qNaN (quiet NaN) と sNaN (signaling NaN) の二種類がありますが、下記では主に(有効な演算を前提として)誤差などについての話をするので、これらの区別はしません。

tips: 複数形には、not numbers や NaNs が使われます*2。sNaN の複数形は sNaNs で回文です。

浮動小数点型が具体的に数値をどう表すかについてはまだ触れていませんが、(固定長のビットを用いて表すことを前提にしている以上は)有限個の値しか表現できないことは明らかです。実数 $x$ に対して、浮動小数点数 $\hat x$ として表現することを 丸め (rounding) と言います。$x$ を丸めた結果が $\hat x$ です。 $\hat x$ としてどのようなものを選ぶかは 丸めモード (rounding mode) によります。これについては後述します。

丸めた結果の絶対値が、無限大になってしまうことを overflow、正規化数(後述)のうち正のものの最小値を下回ってしまうことを underflow と言います。 正規化数の最小値を下回ってもゼロにはせずに非正規化数(後述)とするような方式(およびその現象)を、graceful underflow や gradual underflow と呼びます[Go91, Hi02]

精度という語について

精度という言葉は、適切でない使われ方をすることが多そうな気がします。 英語では precisionaccuracy という異なる 2 つの単語があるのですが、日本語ではどちらも「精度」として訳されがちです。 さらに悪いことに、英語においても precision と accuracy が誤って使われることも多いようです。

precision は「どれだけ細かく表そうとしているか」を指し、accuracy は「実際にどれだけ正確か」を指しています[KD98]。 precision が道具の性能、accuracy がそれを使って得られた結果に対応しているようなイメージでしょうか。

たとえば $3.143$ は十進 4 桁、$3.25678$ は 6 桁ぶんの precision を持ちますが、円周率 $3.1415\ldots$ の近似値としては前者は 3 桁、後者は 1 桁ぶんの accuracy となります*3。accuracy は文脈(たとえばここでは円周率の近似値であるとか)がないことには意味を持たない概念です。 状況に応じて適切に用語を使えるとよいでしょう。

正しい単語とだいたい正しい単語には、稲妻 (lightning) と蛍 (lightning bug) ほどの違いがある。 — Mark Twain

ここでは必要に応じて補足しますが、上記の事情をわかっていれば多くは文脈からわかることでしょう。

他の文献を読むときも、どちらの意味で書かれているかを注意して読むとよさそうです。

他の定義について

[Hi02] においては、precision は基本演算(四則演算など)の accuracy のことだと言っています。 ややしっくり来ない気もしつつ、これは数値計算の文献なので、多少文脈が変わるのかもしれません。

[KD98] の著者には、IEEE 754 の策定にも関わっておられた Kahan 先生も含まれているので、その定義をここでは採用したいと思います。$\eod$

記法

任意の実数 $x\gt 0$ について、$2^e\le x\lt 2^{e+1}$ と表せる唯一の整数 $e$ に対して $2^e = \hfloor x$ と書くことにします[BC04]。 これを $x$ の hyperfloor と呼びます。 たとえば、$\hfloor{2.5} = 2$, $\hfloor{0.4} = 0.25$, $\hfloor{32} = 32$ です。$\hfloor 0$ は定義されません。 $\hfloor x\le x\lt2\hfloor x$ や $\tfrac x2\lt\hfloor x\le x$ などが成り立ちます。

その文脈で考えている丸めモードにおいて、実数 $x$ を丸めた結果を $\roundp x$ と表すことにします。$\roundp{\roundp x}=\roundp x$ や、$x\lt y\implies \roundp x\le\roundp y$ が成り立ちます。 これは自分で勝手に導入した記号なので、一般的なものではないと思います。

他の記法について $\mathit{fl}(x)$ や $\mathrm{fl}(x)$ のような記法も見かけますが、これは「引数の箇所を浮動小数点数として計算する」ことを意味しているようです。たとえば $$\mathit{fl}(x+y\times z) = \roundp{\roundp x + \roundp{\roundp y\times\roundp z}}$$ を意図していそうです[Hi02, OgRuOi05, etc.]。ここでは丸めが起きている箇所を明示したいので、ここでは採用しないことにします。また、$[x]$ のような記法も Kahan に提案されていた[Hi02]ようですが、これは各種記号と紛らわしいのでここでは採用しませんでした。 $\eod$

浮動小数点数 $x$, $y$ における四則演算を、$x\oplus y$, $x\ominus y$, $x\otimes y$, $x\oslash y$ と表すことにします[Go91, Kn14, etc.]浮動小数点数同士で四則演算をした場合、(それを実数で行った場合の)計算結果が該当の浮動小数点型で表せるとは限りません。しかし、IEEE 754 においては、「まず無限精度 (infinite precision) でその演算を行い、その結果を該当の丸めモードで正確に丸めた (correctly rounded) 値を返す」ということが定められています。そこで、次のことが成り立ちます。 $$ \begin{aligned} x\oplus y &= \roundp{x+y}, \\ x\ominus y &= \roundp{x-y}, \\ x\otimes y &= \roundp{x\times y}, \\ x\oslash y &= \roundp{x\div y}. \\ \end{aligned} $$ 実際には無限精度で計算をすることは不可能ですが、何らかの工夫(ここでは触れない)によってそれを再現することが可能のようです*4。 競プロ er にはおなじみの「答えはとても大きくなるので $998244353$ で割ったあまりを求めてください」が、実際の値を求めることなく計算可能なのと似ている気がします。

これらの演算は結合的ではなく、たとえば $(a\oplus b)\oplus c = a\oplus (b\oplus c)$ が成り立つとは限りませんが、簡潔さのため、$a\oplus b\oplus c$ と書いたときは $(a\oplus b)\oplus c$ を指すものとします($\ominus$, $\otimes$, $\oslash$ についても同様)。 優先順位については、通常の四則演算と同様に $\otimes$, $\oslash$ の方が $\oplus$, $\ominus$ よりも先に計算するものとします。

>>> 2**53
9007199254740992

# 結合的でないことを示す例
>>> (9007199254740992.0 + 1.0) + 1.0
9007199254740992.0
>>> 9007199254740992.0 + (1.0 + 1.0)
9007199254740994.0

また、実数の区間を次のように表します[GKP94]。 $$ \begin{aligned} [\alpha\lldot\beta] &= \{x\in\R\mid \alpha\le x\le\beta\}, \\ [\alpha\lldot\beta) &= \{x\in\R\mid \alpha\le x\lt\beta\}, \\ (\alpha\lldot\beta] &= \{x\in\R\mid \alpha\lt x\le\beta\}, \\ (\alpha\lldot\beta) &= \{x\in\R\mid \alpha\lt x\lt\beta\}, \\ \end{aligned} $$ 整数の区間を意味するときは、$[\alpha\lldot\beta)\cap\Z$ のようにして表します。たとえば、$[-1\lldot3)\cap\Z = \{-1, 0, 1, 2\}$ などです。

コードでの表記に関して、多くの言語ではたとえば 2e-3 のように書いて $\roundp{2\times 10^{-3}}$ を表すことができます。Wolfram|Alpha においては 2*^-3 のような表記も見られますが、ここでは用いません。

コード例の文脈では等幅フォント0.1 のように書き、この数値は通常のコードで書いたときのように丸めの対象になるものとします。そうでなく $0.1$ のように書いたときは正確な値を意味するものとします。たとえば、0.1 は $0.1$ のことではなく(その文脈での丸めモードに応じて)${\small 0.100000000000000005}{\scriptsize 55111512312578270211}{\tiny 81583404541015625}$ などを指し、$0.1$ と書いたときは正確に $0.1$ を指します。数式上で丸めを考えたいときは $\roundp{0.1}$ のように明示します。すなわち、地の文の数値が暗黙に丸めの対象になることはないようにします*5

表現について

浮動小数点形式での各値の表現について考えます。NaN は詳しく触れないとして、$+\infty$ と $-\infty$ は難しいことを考える必要はなさそうです(「そういう値」として別途用意しておけばよいので)。

有限値の表現について説明した後、それらを(いわゆる double などの型で持てるように)ビット列にエンコードする方法について話します。 数学的に誤差評価などをする上では、エンコードの話はわからなくても大丈夫なので、目的に応じて読み飛ばしても問題ないでしょう。 浮動小数点型で考えている数自体とそのエンコードは分けて考えてよいということです*6

浮動小数点型の各々の数は特定の実数ひとつを表していて、たとえば (double)3.0 というのは $3$ を正確 (accurate) に表しています。 表している数学的対象としては (int)3 と同じものです。 また、3.03.0 + 1e-3242.99 + 0.0099999999999999 も、全く同じ 3.0 ($=3$) を指します。「これまでどのように計算してきて、どのくらい誤差が積み上がった」のような情報は全く保持していません[KD98]。つまり、「3.0 の方が優れた $3$ で 2.99 + 0.0099999999999999 の方が劣った $3$ である」のような区別は(値からは)不可能・ナンセンスです*7

また、$3$ に丸められる実数の区間(たとえば binary64 で、後述するデフォルトの丸めモードにおいては $[3-2^{-52}\lldot3+2^{-52}]$)を表しているというイメージをしている方もおられるかもしれませんが、やや妥当さに欠ける気がします*8。 たとえば 3.0 + 1.0e-100 == 3.0 であり、浮動小数点数同士の加算を区間の(たとえば両端の)加算と見なそうとしてもうまくいかず、うれしくないためです。 もちろん、何らかの数値計算をした結果が 3.0 になった場合でも、真の値が $[3-2^{-52}\lldot3+2^{-52}]$ に収まっているという保証もありません*9

有限値の表現について

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

パラメータ $\emin$ (minimum exponent), $\emax$ (maximum exponent), $p$ (precision) を定めます。有名な型においては次のような値となります。

型名 $\emin$ $\emax$ $p$
binary32 $-126$ $+127$ $24$
binary64 $-1022$ $+1023$ $53$
binary128 $-16382$ $+16383$ $113$

一般に、$\emin = -(\emax-1)$ となるように決められています。

ゼロでない有限値のうち、浮動小数点型で表せるものは次のように分類されます。

  • $(-1)^s\times (m\cdot 2^{-(p-1)})\times 2^e$
    • $s\in\{0, 1\}$
    • $e\in[\emin\lldot\emax]\cap\Z$
    • $m\in[2^{p-1}\lldot2^p)\cap\Z$
  • $(-1)^s\times (m\cdot 2^{-(p-1)})\times 2^{\emin}$
    • $s\in\{0, 1\}$
    • $m\in(0\lldot 2^{p-1})\cap\Z$

前者を 正規化数 (normal number) と言い、後者を 非正規化数 (denormalized number, subnormal number) と言います。 $(-1)^s$ の部分を 符号部 (sign)、$e$ の部分を 指数部 (exponent)、$m\cdot 2^{-(p-1)}$ の部分が 仮数 (significand, mantissa)と呼ばれます*10

$\gdef\poszero{{0_+}}$ $\gdef\negzero{{0_-}}$

ゼロは +0.0-0.0 が区別されて存在し、正規化数にも非正規化数にも含まれません。ここでは前者を $\poszero$、後者を $\negzero$ と書くことにします。0.0 と書いた場合は $\poszero$ を表します。浮動小数点数としての値を考える文脈では $0$ ではなく $\poszero$ や $\negzero$ と符号を明示することにします。とはいえ、特殊なケースについて考えたい状況を除いては、$\negzero$ が出てくることはあまり多くないです。$-\poszero = \negzero$ が成り立ちます。

ゼロの符号に関して

ゼロを区別しない場合、$1\oslash(1\oslash(-\infty)) = 1\oslash0 = +\infty$ のようにするしかなく、うれしくなさそうです。 規格においては、$1\oslash(1\oslash(-\infty)) = 1\oslash\negzero = -\infty$ のように定められています。 $\lim_{x\to +0} f(x)$ や $\lim_{x\to -0} f(x)$ などに基づいて考えるとよい場合が多い気がします。

絶対値が大変小さい正の値と負の値がそれぞれ $\poszero$, $\negzero$ に対応すると説明している人もしばしばいます。 一方、$\sqrt{\negzero} = \negzero$ が規定されており、やや納得しにくい気もします。 $\negzero$ のべき乗に関しても別途定められており、$(\negzero)^{0.5} = \poszero$ となっています。そうした意味での一貫性はないかもしれません。

最適化に関して、$a\ominus b = -(b\ominus a)$ とすることはできません。具体的には、$b=a$ のとき $a\ominus a = \poszero$ ですが $-(a\ominus a) = \negzero$ となります。また、$\roundp{-x} = \poszero \ominus \roundp x$ とは限りません。たとえば $x = 10^{-324}$ のとき $\roundp{-x} = \negzero$, $\roundp x = \poszero$ なので、左辺は $\negzero$ で右辺は $\poszero\ominus\poszero = \poszero$ となります。

ここでは記号としての $\negzero$ と $\poszero$ を区別するために $\negzero\ne\poszero$ としていますが、-0.0 == +0.0 は true です。これについても後述しますが、記号の等価比較としての $=$ と、浮動小数点型の順序つき比較 == を別で定義することで納得するものとします。

下記の話においては、ゼロの符号が重要になることは一旦ないので、この記事においてはあまり立ち入らないことにします。 とはいえ、浮動小数点型と仲よくなるためには重要なトピックの一つである気もします。$\eod$

binary64 における(絶対値の意味での)最大の正規化数、最小の正規化数、最大の非正規化数、最小の非正規化数はそれぞれ次のようになります。 $$ \begin{aligned} (2^{53}-1)\times 2^{1023-(53-1)} &\approx 1.7976931348623157\times10^{308}, \\ 2^{53-1}\times 2^{-1022-(53-1)} &\approx 2.2250738585072014\times10^{-308}, \\ (2^{53-1}-1)\times 2^{-1022-(53-1)} &\approx 2.2250738585072009\times10^{-308}, \\ 1\times 2^{-1022-(53-1)} &\approx 4.9406564584124654\times10^{-324}. \end{aligned} $$

エンコードについて

$w$ bits の浮動小数点型は、符号部 $1$ bit、指数部 $w-p$ bits、仮数部 $p-1$ bits を持ちます。

正規化数 $(-1)^s\times (m\cdot 2^{-(p-1)})\times2^e$ については、$(s, e+\emax, m-2^{p-1})$ を順に保持します。指数部の最小値が $\emin+\emax = 1$ であることに注意しましょう。また、$m-2^{p-1}\in[0\lldot2^{p-1})$ です。

非正規化数 $(-1)^s\times (m\cdot 2^{-(p-1)})\times2^{\emin}$ については、$(s, 0, m)$ を保持します。指数部に関して $(\emin-1)+\emax = 0$、すなわち正規化数での最小値より $1$ 小さいです。 $\poszero$ については $(0, 0, 0)$ を、$\negzero$ については $(1, 0, 0)$ を保持します。

$\pm\infty$ については、$(-1)^s\times\infty$ と見なして $(s, 2\emax+1, 0)$ を保持します。NaN については、$0$ でない値 $d$ を用いて $(s, 2\emax+1, d)$ を保持します。$d$ に関してはここでは触れません。NaN の符号についても触れないことにします。

表にすると次のようになります。

分類 符号部 指数部 仮数
NaN $s$ $2\emax+1$ $\gt 0$
無限大 $s$ $2\emax+1$ $0$
正規化数 $s$ $e+\emax$ $m-2^{p-1}$
非正規化数 $s$ $0$ $m$
ゼロ $s$ $0$ $0$

整数として読んだ場合の解釈について

エンコードされた値 $(s, e', m')$ に関して、符号 $s$、指数部 $e'$、仮数部 $m'$ をこの順に連結させた整数 $(s'\cdot 2^{w-p} + e')\cdot 2^{p-1} + m'$ を考えます。 このとき、同じ符号ビットを持つもののうちで比較すると、NaN の絶対値が最も大きく、無限大、正規化数、非正規化数、ゼロの順に絶対値が小さくなっています。 NaN 以外に関しては、それらの大小関係と一致していることもわかります。 このことは、浮動小数点数を整数に対応づける場合に役立つことがあります。

rsk0315.hatenablog.com

$\eod$

例として、$\pi$ の近似値について考えます。 $$ \begin{aligned} \pi &\approx 3.141592653589793115997963468544185161590576171875 \\ &= (\text{\texttt{1921fb54442d18}}_{(16)}\cdot 2^{-(53-1)})\times2^{1} \end{aligned} $$ より、 $$ \texttt{\small 0}\: \underbrace{\tt\small 10000000000}_{1+\emax=1024}\: \underbrace{\tt\small 10010010}_{\texttt{92}}\, \underbrace{\tt\small 00011111}_{\texttt{1f}}\, \underbrace{\tt\small 10110101}_{\texttt{b5}}\, \underbrace{\tt\small 01000100}_{\texttt{44}}\, \underbrace{\tt\small 01000010}_{\texttt{42}}\, \underbrace{\tt\small 110100011000}_{\texttt{d18}} $$ がエンコードされたビット列となります。

16 進法リテラルが使える言語では、0x1.921fb54442d18p+1 のように表すこともできます。$[\text{\texttt-}]\text{\texttt{0x}}\{m\cdot 2^{-(p-1)}\}\text{\texttt p}\{e\}$ のような形式です*11。 C なら printf("%a\n", pi)Python なら pi.hex() などでこうした 16 進法の文字列を得ることができます。

bartaz.github.io

こちらのサイトで遊んでみると楽しいかもしれません。上部の入力欄に数値を入れたり、0/1 の部分をクリックしたりすることで、いろいろ楽しむことができます。

丸めについて

丸めに関しては、5 つのモードがあります。

  • 最も近いものへの丸め
    • roundTiesToEven
    • roundTiesToAway
  • 方向に基づく丸め
    • roundTowardPositive
    • roundTowardNegative
    • roundTowardZero

最も近いものへの丸めでは、ちょうど中間にあった場合の挙動により種類が分かれます。 roundTiesToEven では仮数部が偶数の方(状況によってはこれが一意とは限らず、その場合は絶対値が大きい方)に、roundTiesToAway では絶対値が大きい方に丸められます。

方向に基づく丸めでは、名前が示す通り、それぞれ $+\infty$, $-\infty$, $0$ に近い方に丸められます。

規格では、デフォルトの丸めモードは roundTiesToEven となるように定められています。 言語によってはこれを変更する方法が提供されていたりしますが、ここでは詳しくは立ち入りません。C であれば、#include <fenv.h> などに関して調べるとよいでしょう。

そうした事情から、下記の解析では roundTiesToEven を前提とします。また、現在は 64-bit の型が広く使われていることから、$\roundp x$ は binary64 の roundTiesToEven を意味することにします。 理解を深めることで、そうでない状況での解析もできるようになることと思います。

丸めの方法が定義されていることから、丸めモードを変えたり特殊な最適化(GCC-ffast-math など)をしない限りは「コードのこの箇所では $2^{53}\oplus 1$ が $2^{53}$ になったが、別の箇所では $2^{53}\oplus 1$ が $2^{53}+2$ になった。さらに別の箇所では $2^{53}\oplus 1$ が $2^{53}-1$ になった」のようなことにはならず、再現性が保証されることがわかります*12。これにより、浮動小数点数の性質を数学的に記述したり証明したりすることができます。

see also: Semantics of Floating Point Math in GCC

よくある誤差や勘違いの例

基本的な説明をしたので、ある程度込み入った話をしても許される状態になりました。 必要に応じて、戻って読み直してくださればと思います。

0.1 = 1 / 10?

さて、よくある誤りとして、次のようなものがあります。

>>> 0.1
0.1
>>> 1.0 / 10.0
0.1
>>> 0.1 == 1.0 / 10.0
True

「コンピュータでは $0.1$ を正確に表せない*13」というのを聞いて、このように試して「いや、正しく計算できてそうだよ?」と思ってしまう、というものです。実際に上記で示されているのは、$\roundp{0.1} = \roundp{1}\oslash\roundp{10}$ であるということです。

実際には $\roundp{0.1} = 0.1$ は成り立たず、$\roundp{0.1} = \tfrac1{10} \cdot(1+2^{-54})$ となります。 また、$\roundp{1}=1$, $\roundp{10}=10$ であることから、$\roundp{1}\oslash\roundp{10}=1\oslash10=\roundp{0.1}$ とわかります。 次のように出力の桁数を増やすことで確かめられます*14

>>> print(f'{0.1:.100g}')
0.1000000000000000055511151231257827021181583404541015625
>>> print(f'{1.0/10.0:.100g}')
0.1000000000000000055511151231257827021181583404541015625

$\roundp{0.1} = \tfrac1{10}\cdot(1+2^{-54})$ を求める方法については後述することにして、ここでは一旦 fact として認めてもらうことにします。

さらに「正しく計算できているっぽく見える例」として、次のようなものがあります。

>>> 0.1 * 10.0 == 1.0
True

これも、たまたま $\roundp{0.1}\otimes 10 = 1$ が成り立つだけで、$0.1$ の正確さとは関係ありません。事情を知らないと、どこで丸めが発生しているかなどのことに考えが及ばず、誤った解釈をしてしまう人が多そうです。

実際、$\roundp{\tfrac1d}\otimes d = 1$ が成り立つとは限りません。これについては別の記事に書いたので別途紹介しますが、この記事を通して浮動小数点型の算術とお近づきになった後は自力で示せるかもしれません。

exercise: overflow・underflow が起きない前提で、浮動小数点数 $d$ に対して $1-2^{-53}\le \roundp{\tfrac1d}\otimes d\le 1$ となることを示せ。

exercise: $\roundp{\tfrac1{49}}\otimes 49\lt 1$ となることを手元のコンピュータで確かめ、そうなる理由を説明せよ*15。$49$ でうまくいかなかった場合、$41$ や $43$ などについても試してみよ。

“たまたま” $\roundp{0.1}\otimes 10 = 1$ が成り立つ

実際には、非負整数 $i$, $j$ に対して $n = 2^i+2^j$ であるとき、非負整数 $m\le 2^{p-1}$ に対して $m\oslash n\otimes n = m$ が成り立つことが示せるようです[Go91]。$\eod$

0.1 + 0.2 = 0.3?

さて、「うまくいったように見える例」ではなく、「うまくいっていないように見える例」も挙げます。

>>> print(f'{0.1 + 0.2:.100g}')
0.3000000000000000444089209850062616169452667236328125
>>> print(f'{0.3:.100g}')
0.299999999999999988897769753748434595763683319091796875
>>> 0.1 + 0.2 == 0.3
False

$$ \begin{aligned} \roundp{0.1} &= \tfrac1{10}\cdot(1+2^{-54}), \\ \roundp{0.2} &= \tfrac1{10}\cdot(2+2^{-53}) = \tfrac1{10}\cdot(2+2\cdot2^{-54}), \\ \roundp{0.1}\oplus\roundp{0.2} &= \tfrac1{10}\cdot(3+2^{-51}) = \tfrac1{10}\cdot(3+8\cdot2^{-54}), \\ \roundp{0.3} &= \tfrac1{10}\cdot(3-2^{-53}) = \tfrac1{10}\cdot(3-2\cdot2^{-54}). \\ \end{aligned} $$ です。$\roundp{0.1} + \roundp{0.2} = \tfrac1{10}\cdot(3+3\cdot2^{-54})$ に最も近い正規化数は $\tfrac1{10}\cdot(3-2\cdot2^{-54})$ と $\tfrac1{10}\cdot(3+8\cdot2^{-54})$ であり、仮数部が偶数である後者が選ばれています。

もちろん、「小数が絡むと、無限精度では成り立つはずの等式は絶対成り立たない」ということはなく、たとえば $\roundp{0.3} \oplus \roundp{0.4} = \roundp{0.7} = \tfrac1{10}\cdot(7-4\cdot2^{-53})$ が成り立ちます。

整数の誤差

多少わかってきた人は「ははん、小数だと誤差が出るのね」と思うようになってきますが、まだ足りないです*16

>>> int(1e18 + 1)
1000000000000000000
>>> int(1e23)
99999999999999991611392
>>> int(1e25)
10000000000000000905969664

binary64 においては、$\roundp{10^{18}}\oplus 1 = 10^{18}$ や、$\roundp{10^{23}} \ne 10^{23}$ などになるわけですね。 $2^{53}\lt 10^{18}+1$ かつ $5^{18}\lt 2^{53}$ や、$2^{53}\lt 5^{23}$ が成り立つことに注意すると、もう納得できるのではないでしょうか。 もちろん、真の値より小さくなったり大きくなったりします。

Python においては、int 同士の / での除算が float を返すため、意図せず誤差が生じたりもします。

>>> int(2**56 / 3)
24019198012642644
>>> 2**56 // 3
24019198012642645

Rump’s Example

古くから知られている例題もあります。 下記で定義される $f(a, b)$ を考えます。 $$ f(a, b) = 333.75b^6 + a^2 (11a^2b^2 - b^6 - 121b^4 - 2) + 5.5b^8 + \frac a{2b} $$ これに対して $f(77617, 33096)$ を計算せよという問題です[Ru88]。 登場する各定数は binary64 で表せることに注意して、 $$ \begin{aligned} f(a, b) &= 333.75 \otimes b\otimes b\otimes b\otimes b\otimes b\otimes b \\ &\phantom{{}={}} \quad \oplus a\otimes a \otimes (11\otimes a\otimes a\otimes b\otimes b \\ &\phantom{{}={}} \quad\quad \ominus b\otimes b\otimes b\otimes b\otimes b\otimes b \\ &\phantom{{}={}} \quad\quad \ominus 121 \otimes b\otimes b\otimes b\otimes b \\ &\phantom{{}={}} \quad\quad \ominus 2) \\ &\phantom{{}={}} \quad \oplus 5.5\otimes b\otimes b\otimes b\otimes b\otimes b\otimes b\otimes b\otimes b \\ &\phantom{{}={}} \quad \oplus a \oslash(2\otimes b) \end{aligned} $$ となります。原典においては、pow(b, 8) のような関数は用いず、掛け算を順々に行うこと (successive multiplications) を想定しているようです*17

>>> a, b = 77617.0, 33096.0
>>> 333.75*b*b*b*b*b*b + a*a*(
...     11.0*a*a*b*b - b*b*b*b*b*b - 121.0*b*b*b*b - 2.0
...     ) + 5.5*b*b*b*b*b*b*b*b + a/(2.0*b)
1.1726039400531787

みんな大好き Decimal でも計算しておきましょうか。

>>> from decimal import Decimal
>>> a, b = Decimal('77617'), Decimal('33096')
>>> Decimal('333.75')*b*b*b*b*b*b + a*a*(
...     Decimal('11')*a*a*b*b - b*b*b*b*b*b - Decimal('121')*b*b*b*b - Decimal('2')
...     ) + Decimal('5.5')*b*b*b*b*b*b*b*b + a/(Decimal('2')*b)
Decimal('1.172603940053178631858834905')

あるいは、C++long double でも計算してみましょう。113-bit の精度の環境で試します。

#include <cfloat>
#include <cstdio>

long double f(long double a, long double b) {
  return 333.75 * b * b * b * b * b * b +
         a * a *
             (11.0 * a * a * b * b - b * b * b * b * b * b -
              121.0 * b * b * b * b - 2.0) +
         5.5 * b * b * b * b * b * b * b * b + a / (2.0 * b);
}

int main() {
  printf("LDBL_MANT_DIG: %d\n", LDBL_MANT_DIG);
  long double a = 77617.0L, b = 33096.0L;
  printf("%.20Lf\n", f(a, b));
}
LDBL_MANT_DIG: 113
1.17260394005317863186

$1.17260394$ 程度であろうという気分になったところで、真の値について考えてみましょう。 $$ f(a, b) = \frac{1335b^6 + 44a^4b^2 - 4a^2b^6 - 484a^2b^4 - 8a^2 + 22b^8}4 + \frac a{2b} $$ のようにして大部分を 128-bit 整数で計算できることを踏まえると、 $$ \begin{aligned} f(77617, 33096) &= -2 + \tfrac{77617}{66192} \\ &= -\tfrac{54767}{66192} \approx -0.827396059946821368 \end{aligned} $$ となることがわかります。あらあら。long doubleDecimal も信用ならないですね。

こうした誤差が起きる理由は、先の $\roundp{0.1}$ や $10^{18}\oplus 1$ よりは複雑ですが、これについても後述します。

基本的な誤差評価

前述の通り、数値自体は誤差に関する情報を保持していません。 そのため、「このアルゴリズムの出力はこの程度の誤差を含みうる」「真の値は、このアルゴリズムの出力からこの程度離れた範囲には存在する」などといった評価によって保証する必要があります。アルゴリズムの出力がどの程度 accurate か、ということですね。こうした営みなしには、Rump’s example で見たように「精度を上げても同じ結果が得られるからたぶんあっているだろう」のような怪しい祈りをするしかないわけです。

以下では、overflow や (gradual) underflow が起きないものとします。つまり、計算結果が有限の正規化数であるものとします(ただしちょうどゼロになる場合は許容)。 それらが起きた場合は $p$ 桁ぶんの精度 (precision) で計算できなくなり、うれしくない・解析が複雑になるためです。 多くのケースでは正規化数のみの算術が重要なので、簡単のために一旦はここについて考えることにしたいです。

実際には underflow が起きた場合でも成り立つ命題もありますが、ここではあまり触れません。

用語に関して

計算したい値を $x$、実際に計算できた値を $\hat x$ とするとき、$\varepsilon_{\text{abs}} = |\hat x-x|$ を 絶対誤差 (absolute error)、$\varepsilon_{\text{rel}} = \tfrac1x|\hat x-x| = \tfrac{\varepsilon_{\text{abs}}}x$ を 相対誤差 (relative error) と言います。

浮動小数点数 $x = m\cdot 2^{e-(p-1)}$ に対して、その値の最下位桁の大きさ $2^{e-(p-1)}$ を、(その値の)ULP (unit in the last place) と呼びます。$\hfloor x = 2^e$ より、$x$ の ULP は $\hfloor x\cdot 2^{-(p-1)}$ と表せます。 誤差の大きさを表すときに ULP を用いることがしばしばあります。

実数の丸め

任意の実数 $x\ne 0$ に対して $\roundp{-x} = -\roundp x$ が成り立つことは、表現方法から簡単にわかります。 また、$\roundp0=\poszero$ です。以下では $x\gt 0$ について考えるものとします。

何らかの形(文字列表現など)で与えられた実数(現実的には有理数)を浮動小数点数に正確に丸めることを考えます。 その実数を $x$ とし、丸めた結果 $\roundp x$ との誤差について考えます。

$\hfloor x\le x\lt 2\hfloor x$ が成り立つことに注意します。ある整数 $e$ が存在して、$\hfloor x = 2^e$ と表せます。overflow しない前提においては、$2^e$ は浮動小数点数として表せることに注意します。同様に $2\hfloor x = 2^{e+1}$ となります*18。 よって、丸めが起きた際でもそれらの範囲には収まるため、 $$\hfloor x \le \roundp x \le 2\hfloor x$$ となります。

このとき、上記の $e$ に対して、次のいずれかが成り立ちます。

  • ある整数 $m\in[2^{p-1}\lldot 2^p)\cap\Z$ が存在して $\roundp x = m\times2^{e-(p-1)}$
    • $\hfloor x\le \roundp x\lt 2\hfloor x$ のケース ... (1)
  • $m = 2^{p-1}$ に対して $\roundp x = m\times 2^{e+1-(p-1)}$
    • $\roundp x = 2\hfloor x$ のケース ... (2)

これらを整理して、ある整数 $m\in[2^{p-1}\lldot 2^p]\cap\Z$ が存在して $\roundp x = m\times 2^{e-(p-1)}$ と書けることになります。

(1) のケースにおいて、$x\lt\roundp x = m\times 2^{e-(p-1)}$ となる場合は $m\ge 2^{p-1}+1$ かつ $$ (m-1)\times 2^{e-(p-1)} \le x \lt m\times 2^{e-(p-1)} $$ が成り立ちます。$x\ge\roundp x = m\times 2^{e-(p-1)}$ となる場合は $m\le 2^{p-1}-1$ かつ $$ m\times 2^{e-(p-1)} \le x \le (m+1)\times 2^{e-(p-1)} $$ が成り立ちます。(2) のケースでは、$m = 2^{p-1}$ かつ $$ (m-1)\times 2^{e-(p-1)} \lt x \lt m\times 2^{e-(p-1)} $$ が成り立ちます。 以上より、 $$i\times 2^{e-(p-1)} \le x \le (i+1)\times 2^{e-(p-1)}$$ なる $i\in[2^{p-1}\lldot 2^p)\cap\Z$ が一意に定まるので、丸めた際の誤差はその $i$ を用いて $$ \begin{aligned} |x-\roundp x| &= \min{\{(i+1)\times 2^{e-(p-1)} - x, x-i\times 2^{e-(p-1)}\}} \\ &\le \tfrac12\times 2^{e-(p-1)} \\ &= 2^{e-p} = \hfloor x\cdot 2^{-p} \end{aligned} $$ とわかります。これは $x$ の 0.5 ULP 相当です。

あるいは、$|\varepsilon|\le 2^{-p}$ を満たす $\varepsilon$ を用いて $\roundp x = x+\varepsilon\hfloor x$ と書くこともできます。 $\hfloor x\le x$ であることや、$\hfloor x=x$ のときは誤差が生じないことを踏まえて、$|\varepsilon'|\lt 2^{-p}$ を用いて $\roundp x = (1+\varepsilon')x$ とすることもできます。

さらに、$|x-\roundp x|\le\hfloor x\cdot 2^{-p}\le \roundp x\cdot 2^{-p}$ を踏まえると、$|\varepsilon''|\le 2^{-p}$ を用いて $\roundp x = \tfrac x{1+\varepsilon''}$ とすることもできます[Theorems 2.2–2.3 in Hi02]

有理数の丸め

文字列からの変換やリテラルからの変換について考えます。 浮動小数点数とは限らない整数 $u$, $v$ を用いて $\roundp{\tfrac uv}$ と書けます。

$2^k = \hfloor{\tfrac uv}$ とすると、$2^{p-1} = \hfloor{\tfrac uv\cdot 2^{p-1-k}}$ です。 $p-1-k$ の正負に応じて、分子または分母の適切な方を $2^{|p-1-k|}$ 倍して得られる整数を $u'$, $v'$ とします(少なくとも片方は元と同じ整数)。

このとき、$r = u'\bmod v'$ として $r\lt v'/2$ であれば $$\biggroundp{\frac uv} = \Floor{\frac{u'}{v'}}\times 2^{k-(p-1)} = {\frac{u'-r}{v'}}\times 2^{k-(p-1)}$$ となります。また、$r\gt v'/2$ であれば $$\biggroundp{\frac uv} = \Ceil{\frac{u'}{v'}}\times 2^{k-(p-1)} = {\frac{u'+(v'-r)}{v'}}\times 2^{k-(p-1)}$$ となります。$r=v'/2$ のときは仮数部が偶数となるように丸められるので $$ \begin{aligned} \biggroundp{\frac uv} &= \left(2\cdot\Ceil{\frac{\floor{u'/v'}}2}\right)\times 2^{k-(p-1)} \\ &= \left(2\cdot\Floor{\frac{\floor{u'/v'}+1}2}\right)\times 2^{k-(p-1)} \\ &= \left(2\cdot\Floor{\frac{\tfrac{u'+v'}{v'}}2}\right)\times 2^{k-(p-1)} \\ &= \left(2\cdot\Floor{\frac{u'+v'}{2v'}}\right)\times 2^{k-(p-1)} \end{aligned} $$ などと表せます。

評価する文脈によっては、$r\ge v'/2$ の場合は指数部が $k+1$ に繰り上がる場合もあることに注意する必要があるでしょう。

例として、binary64 における $\roundp{0.1}$ について考えます。 $\hfloor{0.1} = 0.0625 = 2^{-4}$ であり、$(u', v') = (1\times 2^{(53-1)-(-4)}, 10) = (2^{56}, 10)$ です。 $2^{56}\bmod 10 = 6 \gt 5$ であることから、 $$ \begin{aligned} \roundp{0.1} &= \frac{2^{56} + (10-6)}{10}\times2^{-56} \\ &= \frac1{10} + \frac{10-6}{10} \times 2^{-56} \\ &= \tfrac1{10}(1 + 2^{-54}) \end{aligned} $$ となります。

文字列がとても長い場合

まず、$\mathit{Pmin}$ と呼ばれるパラメータが各フォーマットに対して定義されています*19

$\mathit{Pmin}$
binary16 $5$
binary32 $9$
binary64 $17$
binary128 $36$

また、拡張の型がある場合は、精度を $p$ bits として $\mathit{Pmin} = 1+\ceil{p\log_{10}(2)}$ とします。

そして、サポートしている型の $\mathit{Pmin}$ の最大値を $M$ とし、$H = M+3$ とします。 $H$ 桁よりも大きい文字列から変換する場合は、一旦 $H$ 桁に丸めてから該当の型に変換することが許容されているようです。

glibc においては、下記の記述があるので correctly rounded であること(すなわち丸めが二回行われないこと)を信頼してよさそうな気がします。

Except for certain functions such as [...], and conversions between strings and floating point, the GNU C Library does not aim for correctly rounded results for functions in the math library [...].

$\eod$

基本演算の丸め

四則演算については正確に丸められることが規定されています。上記と同様にして、$|\varepsilon|\le2^{-p}$ を用いて $x\oplus y = (1+\varepsilon)(x+y)$ などのように書くことができます。$\hfloor{x+y}$ を用いる評価などについても同様です。

また、$x = x'\times2^{e_x}$, $y = y'\times2^{e_y}$ とすると $$ \begin{aligned} x\pm y &= (x'\pm y'\times2^{e_y-e_x})\times2^{e_x}, \\ x\times y &= (x'\times y')\times2^{e_x+e_y}, \\ x\div y &= (x'\div y')\times2^{e_x-e_y} \end{aligned} $$ となります。overflow などを気にしなくてよい前提であれば、$2^k$ を掛ける操作で誤差は生じないため、指数部に関しては一般性を失わずに除去して考えられる場合が多いです。

ここでは、一般の実数 $x$, $y$ から丸めて計算した場合の誤差 $|(x+y)-(\roundp x\oplus \roundp y)|$ ではなく、浮動小数点数 $x$, $y$ 同士の計算の誤差 $|(x+y)-(x\oplus y)|$ について考えるものとします。一般の演算を一回行う場合については、上記のような $(1+\varepsilon)^{\pm1}$ 倍を考えるくらいしかないので、特殊な場合や複数回の場合について述べていきます。

差について

正確に差を表せる十分条件を考えます。すなわち、$x\ominus y=x-y$ となる条件です。 ただし、実数における差が $0$ になった場合については、$x\ominus y\in\{\poszero, \negzero\}$ であれば、正確に表せているということにします。

Claim 1: $\poszero\le\tfrac x2\le y\le 2x$ または $2x\le y\le \tfrac x2\le\poszero$ のとき、$x\ominus y=x-y$ が成り立つ。

Proof

$x=y$ の場合は、$x\ominus y=\poszero$ なので成り立つ。これは、$x=y=\poszero$ または $x=y=\negzero$ のケースも含む。 また、$(x, y) = (\poszero, \negzero)$ と $(x, y) = (\negzero, \poszero)$ については、それぞれ $x\ominus y = \poszero$ と $x\ominus y = \negzero$ であり成り立つ。

そうでないとき、すなわち $x\ne y$ かつ $x$ と $y$ の少なくとも一方の絶対値が $0$ でないとき、前提条件から他方の絶対値も $0$ でなく、$x$ と $y$ は同じ符号を持つことがわかる。


まず、$\poszero\lt\tfrac x2\le y\le 2x$ かつ $x\ne y$ の場合について考える。

整数 $m_x$, $m_y$, $e_x$, $e_y$ を用いて $x = m_x\cdot 2^{e_x}$ および $y = m_y\cdot 2^{e_y}$ ($m_x, m_y\in[2^{p-1}\lldot 2^p)$) とする。 $\tfrac x2\le y\le 2x$ より、$\tfrac12\le \tfrac{m_y}{m_x}\cdot 2^{e_y-e_x} \lt 2$ が成り立つ必要がある。 また、$\tfrac{m_y}{m_x}\in(\tfrac{2^{p-1}}{2^p}\lldot\tfrac{2^p}{2^{p-1}}) = (\tfrac12\lldot 2)$ が成り立つ。

$e_y\le e_x-2$ のとき $\tfrac yx=\tfrac{m_y}{m_x}\cdot 2^{e_y-e_x}\lt 2\cdot2^{-2} = \tfrac12$ となる。 同様に、$e_y\ge e_x+2$ のとき $\tfrac yx\ge 2$ となる。 よって、$e_y \in \{e_x-1, e_x, e_x+1\}$ のケースについて考えればよい。

Case 1 ($e_y = e_x$):

$x-y = (m_x-m_y)\cdot 2^{e_x}$ と表せる。 ある非負整数 $k\lt p-1$ が一意に存在し $2^k \le |{m_x-m_y}| \lt 2^{k+1}$ と書けるので、この $k$ を用いて $$ 2^{p-1} \le |{m_x-m_y}|\cdot 2^{(p-1)-k} \lt 2^p $$ と書ける。すなわち、 $$ x\ominus y = (m_x-m_y)\cdot 2^{(p-1)-k} \times 2^{e_x - ( (p-1)-k)} = x-y $$ となる。

Case 2 ($e_y = e_x-1$):

$x-y = (2m_x-m_y)\cdot 2^{e_x-1}$ と表せる。 ここで、$\tfrac x2\le y$ より $$ \tfrac12\cdot m_x\cdot 2^{e_x} \le m_y\cdot 2^{e_y} = m_y\cdot 2^{e_x-1} $$ であり、$m_x\le m_y$ が従う。すなわち、$2m_x-m_y\le m_x\lt 2^p$ となる。

また、$2m_x-m_y\ge 2\cdot 2^{p-1} - (2^p - 1) = 1$ より、$2m_x-m_y$ は $2^p$ 未満の正整数であることがわかる。 あとは Case 1 同様にして示せる。

Case 3 ($e_y = e_x+1$):

$e_x = e_y-1$ なので、Case 2 より $y\ominus x=y-x$ が成り立つ。また、$x\ne y$ なので $$ \begin{aligned} x\ominus y &= \roundp{x - y} \\ &= \roundp{-(y - x)} \\ &= -\roundp{y-x} \\ &= -(y\ominus x) \\ &= -(y - x) \\ &= x-y \end{aligned} $$ が成り立つ。


最後に、$2x\le y\lt \tfrac x2\lt\poszero$ かつ $x\ne y$ の場合について考える。 $\poszero\lt -\tfrac x2\le -y\le -2x$ であり、Cases 1–3 から $(-x)\ominus(-y) = (-x)-(-y) = -(x-y)$ が成り立つ。 あとは Case 3 同様にして $x\ominus y = x-y$ が示せる。$\qed$

Claim 2: $t$ が $x$ の mod $2^k$ での tail であるとき、$x\ominus t = x-t$ が成り立つ。ただし、$t$ が $x$ の mod $2^k$ での tail であるとは、 $$ x \equiv t\pmod{2^k} $$ かつ $|t|\le 2^{k-1}$ が成り立つことを言う*20。ここで、$k$ は任意の整数である。

なお、整数 $a$, $b$ および整数とは限らない実数 $m$ に対して、$a \equiv b \pmod m$ は、ある整数 $q$ が存在して $a-b = qm$ が成り立つことを意味する。

Proof

$|x| = 0$ のとき $t \equiv 0 \pmod{2^k}$ かつ $|t|\le 2^{k-1}$ が成り立ち、$|t| = 0$ となる。よって、このとき $x\ominus t = x-t\in\{\poszero, \negzero\}$ となる。

以下、一般性を失わず $x\gt 0$ かつ $\hfloor x=2^{p-1}$ とする。このとき $x$ は整数であり、$x\equiv 0\pmod 1$ である。

$k\le 0$ のとき、$x\equiv t\equiv 0\pmod {2^k}$ より($|t|$ の範囲に注意して)$t=0$ が従う。

$k\gt p$ のとき、$0\lt x\lt 2^p \lt 2^k$, $|t|\le 2^{k-1}$, $x-t\equiv 0\pmod{2^k}$ より $x=t$ が従う。

$0\lt k\le p$ のとき、整数 $q$ を用いて $x-t=q\cdot 2^k$ と書ける。$|t|\le 2^{p-1}$ なので $|x-t|\le 2^p$ となるが、$k$ の範囲から $q\le 2^{p-1}$ であるので、$q\cdot 2^k$ は浮動小数点数として表せる。すなわち、$x\ominus t=x-t$ となる。$\qed$

多少知っている人は桁落ちについて気になるかもしれませんが、これについては別で説明します。分けた理由についてもそこで述べます。

複数回の演算

複数回の演算を逐次的に行う際の誤差について考えます。

最初の例として、$a_1 + a_2 + a_3 + a_4$ を計算してみます。 ある $\delta_1, \delta_2, \delta_3\in[-2^{-p}\lldot2^{-p}]$ が存在して $$ ( (a_1 \oplus a_2) \oplus a_3) \oplus a_4 = ( ( (a_1 + a_2)(1+\delta_1) + a_3)(1+\delta_2) + a_4)(1+\delta_3) $$ と書けます。

ここで、 $$ \begin{aligned} 1+\theta_1 &= (1+\delta_3), \\ 1+\theta_2 &= (1+\delta_2)(1+\delta_3), \\ 1+\theta_3 &= (1+\delta_1)(1+\delta_2)(1+\delta_3) \end{aligned} $$ とおくと、 $$ a_1\oplus a_2\oplus a_3\oplus a_4 = (a_1+a_2)(1+\theta_3) + a_3(1+\theta_2) + a_4(1+\theta_1) $$ と書けます。ここで、$\theta$ の添字が $(1+\delta_i)$ の次数に対応していて、これが大きい方が誤差が大きくなりうることを意味します。

より一般に、$(1\pm \delta_i)^{\pm 1}$ を $k$ 個掛けたものの積を $1+\theta_k$ とおきます。$\delta_i$ は、その時々に応じて $|\delta_i|\le 2^{-p}$ の範囲で固定されます。これに対して常に $|\theta_k|\le\gamma_k$ と抑えられるような簡単な形の $\gamma_k$ を定義し、それを用いて解析する方法があります。次のように定義されます[Hi02]。 $$\gamma_k = \frac{2^{-p}\cdot k}{1-2^{-p}\cdot k}.$$

exercise: $\max_{\pm}{(1\pm 2^{-p})^{\pm k}} = (1-2^{-p})^{-k}\le 1+\gamma_k$ を示せ。ここで、$\max_{\pm}$ は、複号の任意の取り方について最大値を求めることを意味している。

exercise: 任意の $k\in [1\lldot 2^p)$ に対して $\gamma_k\lt\gamma_{k+1}$ を示せ。

これを用いると、 $$ \begin{aligned} &\phantom{{}={}} |(a_1\oplus a_2\oplus a_3\oplus a_4) - (a_1+a_2+a_3+a_4)| \\ &= |(a_1 + a_2)\theta_3 + a_3\theta_2 + a_4\theta_1| \\ &\le |a_1\theta_3|+|a_2\theta_3|+|a_3\theta_2|+|a_4\theta_1| \\ &\le |a_1|\gamma_3 + |a_2|\gamma_3 + |a_3|\gamma_2 + |a_4|\gamma_1 \\ &\le (|a_1|+|a_2|+|a_3|+|a_4|)\cdot\gamma_3 \end{aligned} $$ などが得られます。

同様にして、$n$ 個の要素を先頭から順に足していく方法では、誤差が $$ |(a_1\oplus\dots\oplus a_n)-(a_1+\dots+a_n)| = (|a_1|+\dots+|a_n|)\cdot \gamma_{n-1} $$ で抑えられます。相対誤差は $\tfrac{\sum|a_i|}{|{\sum a_i}|}\cdot \gamma_{n-1}$ で抑えられ、この $\tfrac{\sum|a_i|}{|{\sum a_i}|}$ は、この問題の 条件数 (condition number) と呼ばれる値です[Ri66, Hi02, OgRiOi05]。$\operatorname{cond}(\sum, a)$ と書くことにしておきます。

条件数に関して

条件数は、その問題において、入力のズレに対して出力がどの程度変化しうるかを表す値です。 一般の問題に対する具体的な計算方法に関してはここでは割愛します。

ここでは、総和 $\sum_{i=0}^{n-1} a_i$ を求める問題と入力 $a$ のペアとして $\operatorname{cond}(\sum, a)$ と書いていますが、単に $C(a)$ や $\kappa(a)$ のように書かれることもありそうです[Hi02]。$\eod$

また、上記と異なる括弧のつけ方での計算を考えます。$\delta_i$ や $\theta_i$ は適当に固定するものとします。 $$ \begin{aligned} (a_1\oplus a_2)\oplus(a_3\oplus a_4) &= (a_1+a_2)(1+\delta_1) \oplus (a_3+a_4)(1+\delta_2) \\ &= ( (a_1+a_2)(1+\delta_1) + (a_3+a_4)(1+\delta_2))(1+\delta_3) \\ &= (a_1+a_2)(1+\theta_2) + (a_3+a_4)(1+\theta'_2) \end{aligned} $$ より一般に、下記のように隣り合う二個ずつを足していくことで、$1+\delta_i$ の次数を抑えることができます(セグ木のようなイメージ)。 $$ ( (a_1\oplus a_2)\oplus(a_3\oplus a_4))\oplus( (a_5\oplus a_6)\oplus(a_7\oplus a_8)) $$ $k$ 個の総和を求めるときの相対誤差が、$\operatorname{cond}(\sum, a)\cdot\gamma_{\ceil{\log_2(k)}}$ で抑えられます。誤差の係数のオーダーが $\gamma_{\Theta(k)}$ から $\gamma_{\Theta(\log(k))}$ に改善されていますね。$p = 53$ の下では、$\gamma_{10^6} \approx 1.1\times 10^{-10}$、$\gamma_{20} \approx 2.2\times 10^{-15}$ 程度です。$\operatorname{cond}(\sum, a)$ が十分大きい状況では、気をつけた方がいいかもしれません?

加減算のみの場合は、構文木の深さが $\gamma$ の次数に対応しそうなので、あまり深くならないようにするとよさそうに見えます。乗除算が絡む場合は $\gamma$ が木の深さだけからは定まらなさそう(部分木内の乗除算の個数とかも関係しそう)なので一般にはちょっと難しくなるかもしれませんが、構文木を考えるアイディア自体はよさげな気がしています(未検証)。

区間和を求める際にも、累積和を用いるよりもセグ木を用いた方が精度が高くなる気がします。実際には結合的ではない(モノイドではない)ですが、この用途においては問題ない気がします。

なお、最近の研究によって、総和の計算においては $\gamma_k$ として $2^{-p}\cdot k$ を使えることがわかったみたいです[PiRu13]内積の計算ではまた違った値を使えるみたいです。

related keywords: pairwise summation, condition number

補題たち

ここでは正規化数についてのみ考えるものとします。

浮動小数点数 $x$ と $y$ に対して、$x\lt z\lt y$ となるような浮動小数点数 $z$ が存在しないことを示したくなるときはしばしばあります。 直感的に言えば、「$x$ と $y$ は等しいか、隣り合っている」「$x$ と $y$ はくっついている」ということです。 正式な表現があるか知らないので、下記では「くっついている」を用いることにします。

Lemma 3: $|y-x| \lt \hfloor{\min {\{x, y\}}}\cdot 2^{2-p}$ のとき、$x$, $y$ はくっついている。

Proof

一般性を失わず、$x\le y$ とする。このとき $x = m\times 2^{e-(p-1)}$ ($m\in[2^{p-1}\lldot2^p)$) とすると、$x$ より大きい最小の正規化数は $(m+1)\times 2^{e-(p-1)}$ となる。

$(m+1)\times 2^{e-(p-1)} \lt 2^{p-1}\times 2^{(e+1)-(p-1)}$ のとき、その次に大きく最小なのは $(m+2)\times 2^{e-(p-1)}$ である。 $(m+1)\times 2^{e-(p-1)} = 2^{p-1}\times 2^{(e+1)-(p-1)}$ のとき、その次に大きく最小なのは $$ \begin{aligned} (2^{p-1}+1)\times 2^{(e+1)-(p-1)} &= 2^{p-1}\times 2^{(e+1)-(p-1)} + 2^{(e+1)-(p-1)} \\ &= (m+1)\times 2^{e-(p-1)} + 2\cdot 2^{e-(p-1)} \\ &= (m+3)\times 2^{e-(p-1)} \end{aligned} $$ である。

よって、$y\lt (m+2)\times 2^{e-(p-1)} = x + \hfloor x\cdot 2^{2-p}$ を示せば十分とわかる。$\qed$

また、特定の浮動小数点数を固定し、その値に丸められる実数の範囲を考えたいこともしばしばあります。

Lemma 4: 実数 $x$ と浮動小数点数 $\hat x$ に対し、$\roundp x=\hat x$ となる十分条件として、 $$ \hat x-\hfloor{\hat x}\cdot 2^{-1-p} \le x \lt\hat x + \hfloor{\hat x}\cdot 2^{-p} $$ がある。$\hfloor{\hat x}\lt\hat x$ であれば $$ \hat x-\hfloor{\hat x}\cdot 2^{-p} \lt x \lt\hat x + \hfloor{\hat x}\cdot 2^{-p} $$ とできる。

Proof

$\hat x = m\times 2^{e-(p-1)}$ ($m\in[2^{p-1}\lldot 2^p)$) とし、$\hat x$ より小さい最大の浮動小数点数について考える。 $m = 2^{p-1}$ のときは、$(2^p-1)\times 2^{(e-1)-(p-1)}$ である。$m\gt 2^{p-1}$ のときは $(m-1)\times 2^{e-(p-1)}$ である。

それらと $\hat x$ との距離は、前者は $1\times2^{(e-1)-(p-1)}=2^{e-p}$、後者は $1\times2^{e-(p-1)}$ である。 丸めモードは roundTiesToEven を仮定しているので、距離の半分未満であれば、$\hat x$ に丸められるのに十分である。 また、$m = 2^{p-1}$ が偶数であることから、そのケースにおいては距離がちょうど半分でも $\hat x$ に丸められる。

前者の方が厳しい条件であることから、$\hat x-\hfloor{\hat x}\cdot 2^{-1-p} \le x$ であれば $\roundp x=\hat x$ となる。

上からの評価についても、上記同様の議論をすることで得られる。$\qed$

$\hfloor x\le x$ などと合わせて使うことも多いでしょうか。

桁落ちについて

よく「値が近い 2 数を引くと桁落ちによって誤差が出る」という説明がなされがちです。 しかし、この説明はやや不十分です。たとえば、binary64 において隣り合う 2 数 $(1+2^{-52})$ と $1$ の差 $2^{-52}$ は binary64 で正確に表せるため、誤差は生じません。

一方で $\roundp{1+2^{-53}}$ と $\roundp{1+2^{-53}+2^{-54}}$ の差も、上記の 2 数に丸めてから計算されるため、$2^{-54}$ ではなく $2^{-52}$ となってしまいます。これを踏まえると、引数たちが誤差を含む値かどうかが肝になっていそうです*21

これらの入力に関して

実際には ${1+2^{-53}+2^{-54}}$ などを直接得ることはできません。一方で、$2^{-53}+2^{-54}$ などは計算することができるため、文字列表現を得てから $1$ を足すなどの方法で対処することが可能です。

>>> x = f'{2**-53 + 2**-54:.54f}'.replace('0.', '1.')
>>> print(x)
1.000000000000000166533453693773481063544750213623046875
>>> print(float(x))
1.0000000000000002
>>> print(f'{float(x):.100g}')
1.0000000000000002220446049250313080847263336181640625

あるいは、$1\oplus(2^{-53}+2^{-54})$ として得る方法もあります。 $1\oslash49\otimes49 = 1-2^{-53}$ を踏まえると、$3\otimes(0.5\ominus(0.5\oslash49\otimes49)) = 2^{-53}+2^{-54}$ となります。

もちろん、Decimal などを使う方法もあります。$\eod$

そこで、実数 $x+\delta$ と $x$ の差 ($|\delta|\ll x$) を計算する状況を考えます。引数が丸められることに気をつけてください。ある $\varepsilon$, $\varepsilon'$ ($|\varepsilon|, |\varepsilon'|\le2^{-p}$) が存在して次のようになります。 $$ \begin{aligned} \roundp{x+\delta}-\roundp x &= (1+\varepsilon)(x+\delta) - (1+\varepsilon')x \\ &= (x+\delta) + \varepsilon(x+\delta) - x - \varepsilon'x \\ &= \delta + \varepsilon(x+\delta) - \varepsilon'x \\ &= (1+\varepsilon)\delta + (\varepsilon-\varepsilon')x. \end{aligned} $$ 実際には $\roundp{x+\delta}\ominus\roundp x$ なので $\roundp{(1+\varepsilon)\delta + (\varepsilon-\varepsilon')x}$ となるわけですが、一旦は上記の値に注目します。

$(1+\varepsilon)\delta$ の項については、$\delta$ に対して微小な誤差なのでよいでしょう。 一方、$(\varepsilon-\varepsilon')x$ については、本来消えていてほしいはずの $x$ についての項になっています。$\varepsilon-\varepsilon'$ は $2^{1-p}$ 程度に大きくなりうるので、$\delta$ の小ささによっては無視できない大きさになりそうです。

Example 5: $\sqrt{2^{52}+2} - \sqrt{2^{52}+1}$

たとえば、$x+\delta=\sqrt{2^{52}+2}$, $x=\sqrt{2^{52}+1}$ について考えてみましょう。 $$ \begin{aligned} \roundp{\sqrt{2^{52}+2}\,} &= 67108864.00000001490{\small 116119384765625}, \\ \roundp{\sqrt{2^{52}+1}\,} &= 67108864.0 \end{aligned} $$ とのことで、 $$ \begin{aligned} \roundp{\sqrt{2^{52}+2}\,} \ominus \roundp{\sqrt{2^{52}+1}\,} &= 1.490116119384765625\times 10^{-8} \\ &= 2^{-26} \end{aligned} $$ となります。 実際の値は $\delta\approx 7.450580596923826884\times 10^{-9} \approx 2^{-27}$ 程度であるため、$2\delta$ 程度の値が出てきてしまっていることがわかります。

この例においては、次のように式変形したものを計算することで、絶対値が近いかつ誤差を含む値同士の引き算を回避することができます。 $$ \begin{aligned} \sqrt{2^{52}+2} - \sqrt{2^{52}+1} &= \frac{(2^{52}+2) - (2^{52}+1)}{\sqrt{2^{52}+2}+\sqrt{2^{52}+1}} \\ &= \frac{1}{\sqrt{2^{52}+2}+\sqrt{2^{52}+1}}. \end{aligned} $$ 実際、 $$ \frac{1}{\roundp{\sqrt{2^{52}+2}}\oplus\roundp{\sqrt{2^{52}+1}}} = 7.450580596923828125\times10^{-9} $$ であり、$\delta$ に近い値になっています。あるいは、$\sqrt{2^{52}+1}+\sqrt{2^{52}+2} \approx 2\cdot\sqrt{2^{52}}$ を踏まえて、$\tfrac1{2\sqrt{2^{52}}} = 2^{-27}$ に近いことも確認できます。$\eod$

Example 6: $(3\cdot2^{51} + 2)^2 - (3\cdot2^{51} + 1)^2$

別の例として、$a\approx b$ を満たす浮動小数点数 $a$, $b$ に対して $a^2-b^2$ を計算してみることを考えます。 $a = 3\cdot2^{51} + 2$, $b = a - 1 = 3\cdot2^{51} + 1$ とします。これらは正確に表せる浮動小数点数であることに注意します。 $$ \begin{aligned} a\otimes a &= \roundp{9\cdot2^{102} + 3\cdot 2^{53} + 4} \\ &= \roundp{2^{105} + 2^{102} + 2^{54} + 2^{53} + 2^2} \\ &= 2^{105} + 2^{102} + 2^{54} + 2^{53}, \\ b\otimes b &= \roundp{9\cdot2^{102} + 3\cdot 2^{52} + 1} \\ &= \roundp{2^{105} + 2^{102} + 2^{53} + 2^{52} + 2^0} \\ &= 2^{105} + 2^{102} + 2^{54} \end{aligned} $$ より、$(a\otimes a)\ominus(b\otimes b) = 2^{53}$ とわかります*22。一方、真の値は $$ \begin{aligned} a^2-b^2 &= (9\cdot 2^{102} + 3\cdot 2^{53} + 4) - (9\cdot 2^{102} + 3\cdot 2^{52} + 1) \\ &= 3\cdot 2^{52} + 3 \end{aligned} $$ です。$2^{53} \approx 0.667 \cdot (3\cdot 2^{52})$ であり、33 % 程度の誤差が出ていることがわかります。 $a^2-b^2 = (a-b)(a+b)$ に基づいて計算すると、 $$ \begin{aligned} (a\ominus b)\otimes(a\oplus b) &= 1\otimes\roundp{3\cdot 2^{52}+3} \\ &= 3\cdot 2^{52}+4 \end{aligned} $$ であり、近い値が得られます。前者の方法では $a^2$ や $b^2$ を計算して誤差を生じさせた後、それらの差を計算しているため、過剰な誤差が出ています。 後者の方法では $a+b$ の計算による誤差程度しか生じず、高い精度で計算できています。

一般に、$(a+\delta_a)(b+\delta_b) \approx ab + a\delta_b + b\delta_a$ より、多少の誤差 ($\delta_a\ll a$, $\delta_b\ll b$) を含む値同士を掛け合わせても、誤差が主要項に大きく影響することはなさそうです。

なお、$(a, b) = (2^{53}, 1)$ のようなケースにおいては、$(a\otimes a)\ominus(b\otimes b)$ で計算した方が精度が高くなることに注意しましょう。$\eod$

最終的な計算結果として欲しいものが $\delta$ で、かつ絶対誤差が小さければよい状況では、$\delta$ 自体が小さければ桁落ちが起きても問題ないかもしれません*23。 主要項が $\delta$ に大きい影響を受けるような状況だと、桁落ちの影響を無視できず、基本的にうれしくないことになりそうです。

絶対値が近い 2 数の差に関して、引数が誤差を含まないものを benign cancellation、誤差を含んでいるものを catastrophic cancellation (severe cancellation) と呼ぶようです[Go91]

Re: Rump’s example

さて、先ほどの例に関して見てみます。$\tfrac a{2b}$ の項が十分に小さいことは見た通りなので、前半の項について確認してみます。

$$ f(a, b) = 333.75b^6 + a^2 (11a^2b^2 - b^6 - 121b^4 - 2) + 5.5b^8 + \frac a{2b} $$ 括弧内をまず計算します。 $$ \underbrace{a^2 (\underbrace{\underbrace{\underbrace{\underbrace{11a^2b^2}_{\tiny 72586759116001040064} - b^6}_{\tiny -1314174461784456350457997632} - 121b^4}_{\tiny -1314174606957974558362483008} - 2}_{\tiny -1314174606957974558362483010})}_{\tiny -7917111779274712207494296632228773890} $$ 全体としてはこうなります。 $$ \underbrace{\underbrace{\underbrace{333.75b^6\mathstrut}_{\tiny 438605750846393161930703831040} + \underbrace{a^2 (11a^2b^2 - b^6 - 121b^4 - 2)}_{\tiny -7917111779274712207494296632228773890}}_{\tiny -7917111340668961361101134701524942850} + \underbrace{5.5b^8\mathstrut}_{\tiny 7917111340668961361101134701524942848}}_{-2} $$ 各項は、計算途中で(絶対値が)$7.917\times10^{36}$ と大変大きくなり、$2^{122}$ を上回りますが、最終的には $-2$ まで小さくなっています。

なお、べき乗の部分を先に計算すると binary64 においては $-1.18\times 10^{21}$ 程度になります。

>>> a, b = 77617.0, 33096.0
>>> 333.75*b**6 + a**2*(11.0*a**2*b**2 - b**6 - 121.0*b**4 - 2.0) + 5.5*b**8 + a/(2.0*b)
-1.1805916207174113e+21

$|{-1.18\times10^{21}}| \approx 7.917\times10^{36}\times2^{-52}$ であり、別の精度で計算するとバラバラの値になります。 この莫大な誤差が生じることに触れて Rump’s example が紹介されることがしばしばあるようですが、元論文には下記のような記述があり、そういう解釈は意図と異なっていそうです。

The obtained values are the following:

     single precision    :   f  =  + 1.172603 ...
     double precision    :   f  =  + 1.1726039400531 ...
     extended precision  :   f  =  + 1.172603940053178 ...

All three values agree in the first 7 figures, whereas the
true value for f is

     exact value         :   f  =  - 0.827396059946821[3..4]

精度をある程度上げて計算しても答えがほぼ一致しているにも関わらず、真の値はどれとも異なっているという例として導入されていそうです。

後続の論文においては、下記の表現が提案されています[LW02]。 $$ f(a, b) = (333.75-a^2)b^6 + a^2 (11a^2b^2 - 121b^4 - 2) + 5.5b^8 + \frac a{2b} $$ これに従うと、べき乗を先に計算しても “所望” の結果を得られます。

>>> a, b = 77617.0, 33096.0
>>> (333.75-a**2)*(b**6) + a*a*(11.0*(a**2)*(b**2) - 121.0*(b**4) - 2.0) + 5.5*(b**8) + a/(2.0*b)
1.1726039400531787

なお、$77617^2 = 5.5\cdot33096^2 + 1 \,(= 6024398689)$ が成り立つので、$a^2 = 5.5b^2+1$ として整理することで次のものが得られます。 $$ f(a, b) = -5.5b^8 - 2 + 5.5b^8 + \tfrac a{2b} $$ この $5.5b^8\approx 7.917\times10^{36}$ によって $-2$ が無視された後に $5.5b^8$ が打ち消されることが本質で、$(-5.5b^8 \ominus 2) \oplus 5.5b^8$ を $-2$ ではなく $\poszero$ としてしまうことにより、$\tfrac a{2b} \approx 1.1726$ が全体的な出力となってしまっています。

>>> -5.5*b**8 - 2 + 5.5*b**8 + a/(2*b)
1.1726039400531787
>>> -2 + a/(2*b)
-0.8273960599468213

この $\pm5.5b^8$ たちがうまく打ち消し合ってくれないと、$5.5b^8\times 2^{-52} \approx 10^{21}$ 程度の誤差が生じてしまい、趣旨に反してしまいそうです(精度を変えると結果が変わって、バレてしまうため)。

精度を十分に上げる(113-bit 精度の long double などよりもさらに高くする)ことで、$(-5.5b^8\ominus 2)$ の部分で誤差が生じなくなり、正しい値が得られるようになります。計算途中の値の絶対値がどの程度の大きさになるかを意識する(あるいは大きくならないようにする)ことが重要です。

融合積和

別の例について考えます。

$(2 - 2^{-52})^2 = 4 - 2^{-50} + 2^{-104}$ であり、$(2-2^{-52})\otimes(2-2^{-52}) = 4 - 2^{-50}$ が成り立ちます。 そのため、$(2 - 2^{-52}) \otimes (2 - 2^{-52}) \ominus (4 - 2^{-50}) = \poszero$ となります。 すなわち、無限精度での計算と比べて $2^{-104}$ の情報が失われています。 $p$ 桁同士の乗算は $2p$ 桁程度の情報を持つのに対し、丸めの結果 $p$ 桁になってしまうため、後から下位 $p$ 桁を取り出すことはできないわけです。

そこで、規格で定義されている別の演算を紹介します。$a\otimes b\oplus c$ ではなく $\roundp{a\times b+c}$ を計算する命令です。 主に fma (fused multiply-add) と呼ばれることが多い気がしますが、fmac, fused mac (fused multiply-accumulate) などとも呼ばれることがあるようです。

これを用いることで、桁落ちによる誤差を回避して $$\roundp{(2-2^{-52}) \times (2-2^{-52}) + (-(4-2^{-50}))} = 2^{-104}$$ を得ることができるわけです。

とはいえ、残念ながら、fma が正しく実装されていない状況も多々あるようですので、困ります。

qiita.com

数学関数に関する式の計算

$f(x) = \tfrac1x(e^x-1)$ の計算を考えてみます。 $x\approx 0$ のとき $e^x\approx 1$ ですから、$e^x-1$ の部分で catastrophic cancellation が生じそうです。 真の値について、 $$ e^x = 1 + x + \frac{x^2}2 + O(x^3)\quad \text{as}\;x\to 0 $$ です*24から、 $$ \frac{e^x-1}x = 1 + \frac x2 + O(x^2)\quad \text{as}\;x\to 0 $$ で、$x\approx 0$ であれば $\tfrac1x(e^x-1) \approx 1+\tfrac x2$ であることがわかります。

$\gdef\codeexp{\mathtt{exp}}$ $\gdef\codelog{\mathtt{log}}$

実際に計算してみます。exp(x) は $\roundp{\exp(x)}$ を計算してくれるとは限らず、多少の誤差が出るかもしれないので、数式上では $\codeexp(x)$ と書いて区別します。$\exp(x)\approx1$ を確認しやすいように $\codeexp(x)$ の値も合わせて見てみます。

>>> from math import exp
>>> f = lambda x: (exp(x) - 1.0) / x
x exp(x) f(x)
1.0 2.718281828459045 1.718281828459045
1e-1 1.1051709180756477 1.0517091807564771
: : :
1e-7 1.000000100000005 1.0000000494336803
1e-8 1.00000001 0.999999993922529
1e-9 1.000000001 1.000000082740371
: : :
1e-14 1.00000000000001 0.9992007221626409
1e-15 1.000000000000001 1.1102230246251565
3.16e-16 1.0000000000000002 0.7026728003956687
1e-16 1.0 0.0

案の定、うれしくなさそうな値が得られていることがわかります。

ところで、$(\codeexp(x)\ominus1)\oslash(\codelog(\codeexp(x))$ に基づいて計算した方が、高い精度 (accuracy) の結果が得られることが知られています。ゼロ除算には注意しつつ、$\codelog(\codeexp(x))$ を計算した値も出力してみます。

>>> from math import exp, log
>>> f = lambda x: (exp(x) - 1.0) / x_ if (x_ := log(exp(x))) > 0.0 else 1.0
x exp(x) f(x)
1.0 2.718281828459045 1.718281828459045
1e-1 1.1051709180756477 1.0517091807564762
: : :
1e-7 1.000000100000005 1.0000000500000017
1e-8 1.00000001 1.000000005
1e-9 1.000000001 1.0000000005
: : :
1e-15 1.000000000000001 1.0000000000000004
3.16e-16 1.0000000000000002 1.0000000000000002
1e-16 1.0 1.0

先ほどのアルゴリズムのときよりも $1 + \tfrac x2$ に近い値が得られていそうです。

理由について説明します。$g(y) = (y-1)/\log(y)$ を考えると、$f(x) = g(e^x)$ です。 $\hat y=\codeexp(x)$ とすると、$\hat y\approx e^x$, $g(\hat y)\approx f(x)$ および $\roundp{\hat y}=\hat y$ が成り立ちます。

誤差のない値に対しては桁落ちが発生しても問題ないため、$\hat y\ominus 1$ を計算しても大丈夫です。そのため、$g(\hat y)$ は十分によい精度で計算することができます。

誤差の削減に関して

総和計算

Kahan algorithm や Kahan–Babuška algorithm など多数の亜種があります。 精度(誤差のオーダー)や計算量などにそれぞれ違いがあります。ここでは、Neumaier による改良版の Kahan–Babuška algorithm を挙げます。

function $\textsc{sum-compensated}(a)\colon$
$s \gets 0$
$c \gets 0$
for $a_i$ in $a$ do
    $s' \gets s\oplus a_i$
    if $|a_i| \le |s|$ then
        $c \xgets{\oplus} (s \ominus s')\oplus a_i$
    else
        $c \xgets{\oplus} (a_i \ominus s')\oplus s$
    end-if
    $s \gets s'$
end-for
return $s\oplus c$

ここで、$c$ は補正 (compensation) のための項で、$s\oplus a_i$ の誤差項の和を持っている項です。これによって、相対誤差を $\operatorname{cond}(\sum, a)\cdot\gamma_{|a|-1}^2$ 以下に抑えることができます。補正を複数段で行うことで、さらに精度を上げることも可能なようです[OgRuOi05]。 証明については次回以降の記事で改めて取り上げることにします。

Python 3.12 では sum() がこのアルゴリズムが採用され、誤差の上界が小さくなりました (cpython/Python /bltinmodule.c)。 float でない値が混じっている場合は期待したようにならないので、float での総和を計算したいときは各値を float だけにしておくのがよいかもしれません。 なお、correctly-rounded な値を返すようになったというわけではないことに注意してください (cf. math.fsum(iterable))。

% () { python3.11 <<< $1; echo; python3.12 <<< $1 } 'import sys
print(sys.version)
print(sum([0.1] * 10))
a, b = 77617., 33096.
print(sum([-5.5*b**8, -2.0, 5.5*b**8, a/(2*b)]))
print(sum([1e20, 1e40, 1.0, -1e40, -1e20]))
print(sum([1e40, 1e20, 1.0, -1e20, -1.0, -1e40]))'

出力は下記の通りです。

3.11.5 (main, Aug 24 2023, 15:09:45) [Clang 14.0.3 (clang-1403.0.22.14.1)]
0.9999999999999999
1.1726039400531787
-1e+20
0.0

3.12.0 (tags/v3.12.0:0fb18b02c8, Nov 12 2023, 15:16:03) [Clang 14.0.0 (clang-1400.0.29.202)]
1.0
-0.8273960599468213
0.0
1.0

最後の行のように、問題によっては v3.11 の方だけが(たまたま)正確な値になっていることもありますが、値自体は精度保証をしていないことを思い出しましょう(すなわち、v3.12 でたまたま値が正確でなくなったケースを指して改悪だというのは、やや短絡的そうということ)。とはいえ、sum などは広く使われていそうなのに、関数名やオプションなどではなくデフォルトの挙動でこうした修正をするのは、なかなか攻めたアップデートだなという気持ちもあります。

related topics: error-free transformations, Dekker’s algorithm, compensated summation

数学関数の精度について

$\text{\texttt{sqrt}}(x) = \roundp{\sqrt{x}}$ や $\text{\texttt{fma}}(x, y, z) = \roundp{x\times y+z}$ などは correctly rounded であることが規定されていますが、それ以外の数学関数 (sin, log, atan, cbrt, etc.) については、そうなるとは限りません。 IEEE 754 が定める関数は、required functions と recommended functions に分かれており、四則演算などは前者、三角関数などは後者に含まれます。 前者は「それらの演算の correctly-rounded な値を用意する必要がある」であるのに対し、後者は「それらの数学関数の correctly-rounded な値を得る関数を提供することが好ましい」というもので、提供していなくても IEEE 754 に準拠していないことにはなりません。

たとえば glibc においては、それらを正確に返すことを目標にしていないという記述があり、どの程度の誤差が出るかの表を公開しています。 また、無限精度では「この関数は(この区間で)単調増加である」などの性質がある場合でも、それが満たされるとは限らないようです。

related topics: table maker’s dilemma, Lindemann–Weierstrass theorem

実際の例

wandbox 上で gcc 8.4.0 と gcc 9.3.0 で、std::norm(std::complex<long double>) の挙動が変わる例を最近見たので、紹介しておきます。

#include <cfloat>
#include <complex>
#include <cstdio>

int main() {
    fprintf(stderr, "LDBL_MANT_DIG: %d\n", LDBL_MANT_DIG);
    fprintf(stderr, "%.100Lg\n", std::norm(std::complex<long double>(49921.0L, 50077.0L)));
    fprintf(stderr, "gcc %d.%d.%d\n", __GNUC__, __GNUC_MINOR__, __GNUC_PATCHLEVEL__);
}

各バージョンの出力はそれぞれ次の通りです。

LDBL_MANT_DIG: 64
4999812169.9999999995343387126922607421875
gcc 8.4.0
LDBL_MANT_DIG: 64
4999812170
gcc 9.3.0

ICPC 2023 Asia Yokohama Regional K 問題の AOJ ミラー においては、これに起因して(おそらくは)AOJ の gcc が 8.3.1 である関係で、

answer 50002 50077 49921

が正答となるテストケースにおいて

query 81 0 81 100000

0.0000432 が返ってくるようです(真の値が $0$ で絶対誤差が $10^{-6}$ より大きいので不正)。Juding Details 的には gcc 11.4.0 のようなので、本番のジャッジ環境では(少なくともこのケースでは)問題なさそうです。$\eod$

なお、core-mathcr-libm のように、correctly-rounded な数学関数を実装しているライブラリもあるようです。後者はもうメンテされていないようです、かなしい。

比較演算について

$\gdef\nan{\bot}$ $\gdef\imcomp{\mathrel{\,\xcancel{\phantom:}\xcancel{\phantom:}\,}}$

浮動小数点型で表せる値全体からなる集合 $F$ を考えたり、それらの演算を定義する上では、$=$ を形式的に定義する方が無難そうです。 たとえば、$\poszero\ne\negzero$ や $\nan=\nan$ などです。こうしておくことで、$\poszero\oslash\poszero = \nan$ のようにして記述することもできます($\nan$ で NaN を表しています)。

一方、浮動小数点型の演算においては +0.0 == -0.0true だったり、nan == nanfalse だったりします。 これに関しては、$F$ 上の二項関係 $\le_F$ を別途定義して、それに基づく比較を == などで表していると見るのがよい気がしています。

全ての要素間での比較が定義されているわけではなく、less than ($\lt_F$), equal to ($=_F$), greater than ($\gt_F$) に加えて unordered ($\imcomp_F$) となる場合もあります。 たとえば $\negzero =_F \poszero$ (-0.0 == 0.0) や、$\nan\imcomp_F\nan$ (nan != nan) や、$1\imcomp_F\nan$ (1.0 != nan) となります。 unordered である要素同士を <, <=, == などで比較したときは false が返ってくるため、たとえば x <= y!(x > y) が等価とは限らないことに気をつけましょう。 両辺がどちらも $\nan$ でない場合は unordered にはならないので安心です。

演算子について、どの条件が満たされるときに true が返るかの表を示します。✅ が true、❌ が false を意味します。

$\lt_F$ $=_F$ $\gt_F$ $\imcomp_F$
<
<=
==
>=
>
!=

たとえば、x <= y は $(x\lt_F y)\vee(x=_F y)$ と同値、x != y は $(x\lt_F y)\vee(x\gt_F y)\vee(x\imcomp_F y)$ と同値です。$\imcomp_F$ に注意しましょう。

こうした局面においては bool ではなく Option<bool> などが返ってきてほしい気もします。 Rust の .partial_cmp() においては Some(Less), Some(Equal), Some(Greater), None を返すようになっており、よさがあります。

雑多トピック

Decimal を使いましょう」について

Python で「Decimal であれば精度が上がる」という主張をしている人はよく見かけます。この精度というのは precision なのか accuracy なのかどっちの意図で言っているのでしょう。おそらくは特に区別せず言っていると思うのですが。

binary64 での 0.1 では $0.1$ を正確に表せないが Decimal('0.1') は正確に $0.1$ なので Decimal の方が精度が高い、という主張も(残念ながら)しばしば見ます。$\tfrac1{10}$ が表せているだけで、たとえば $\tfrac13$ などはやはり正確には表せないことは理解されているのでしょうか。

十進の小数表示に対して(適切な桁数設定の下で)正確に計算したいという状況では問題ないでしょうが、「いつでも Decimal を使っておけば常にいい感じになる」と解釈されかねない主張は、教育としても有害のような気がします。 実際、先の Rump’s example においてもデフォルトの精度 (getcontext().prec == 28) において正しくない値を返すことがわかりました。

long double を使ってみよう」について

double だと WA になったのでとりあえず long double に変えてみよう」という 愚かな ことをしたことがある競プロ er は多数いることと思います。

たとえば、double では区別できなかった 2 数を区別するために long double を使いたい(たとえば Sqrt Inequality の正しい long double 解法など)のであれば正当そうな感じはします。

丸め誤差によって真の値より大きくなったり小さくなったりしていることを解決したい場合は、(precision を上げれば正確に表せる場合を除いて)不当そうな感じがします。 ここでは double (53-bit precision) による丸めを $\roundp{x}_{53}$、long double (64-bit precision) による丸めを $\roundp{x}_{64}$ と表すことにしてみると、たとえば $$ \gdef\sroundp#1{(\hspace{-.15em}[#1]\hspace{-.15em})} \begin{aligned} % \tfrac1{100}(1+12\cdot2^{-59}) = \roundp{0.01}_{53} &\gt 0.01 \gt \roundp{0.01}_{64} \gt \tfrac1{100}(1-24\cdot 2^{-70}) \\ % \tfrac1{10^6}(1-213696\cdot 2^{-72}) = \roundp{0.1^6}_{53} & \lt 0.1^6 \lt \roundp{0.1^6}_{64} \lt \tfrac1{10^6}(1+350592\cdot 2^{-83}) \underbrace{\tfrac1{100}(1+12\cdot2^{-59})}_{\sroundp{0.01}_{53}} &\gt 0.01 \gt \underbrace{\tfrac1{100}(1-24\cdot 2^{-70})}_{\sroundp{0.01}_{64}} \\ \underbrace{\tfrac1{10^6}(1-213696\cdot 2^{-72})}_{\sroundp{0.1^6}_{53}} & \lt 0.1^6 \lt \underbrace{\tfrac1{10^6}(1+350592\cdot 2^{-83})}_{\sroundp{0.1^6}_{64}} \end{aligned} $$ となり、大小関係が逆転しています。すなわち、precision を上げた方が値が大きくなる場合や小さくなる場合が(当然)あり、何らかのよい性質を期待するのは妥当ではないように思えます。

double だと期待した値より小さくなってしまったので long double を使ったところ、期待した値以上になり、うまくいった」というのは、long double なら別の反例で落ちるであろうということです。

しかし、「53-bit precision を落とすケースは用意していても 64-bit precision や 113-bit precision を落とすケースは用意していないだろう」という期待であれば実る可能性が(残念ながら)そこそこ高いような気もします。ただ、多くの場合は嘘解法でしょうから、コンテスト後に撃墜ケースを自分で作ってみることをおすすめします。

もちろん、「Decimal を使ってみよう」に対しても上記のことが当てはまります($0.1$ の例に関しては $1/3^k$ などを考えてください)。

なお、

long double eps = 1.0e-6;

のように書くと、double 型のリテラルであるところの 1.0e-6 を計算したものが long double にキャストされて eps に入ることに注意してください。すなわち、$\roundp{10^{-6}}_{64}$ ではなく $\roundp{10^{-6}}_{53}$ が入っているということです。long doubleリテラル1.0e-6L のような形式です。

printf に関して

保持している値の確認などの際に、出力の挙動を勘違いしているとすべてを勘違いしそうなので、気をつける必要がありそうです。

Java や Kotlin での挙動

%f による指定では、Double.toString(double) で返ってきた文字列を、桁数に応じて round half-up で丸めたりゼロ埋めしたりすると書かれています。 Double.toString(double) では、double において隣接する値と区別するのに足りる長さしか返ってこなさそうなので、double が表している値を正確に得ることは期待できなさそうです。

先の $\mathit{Pmin}$ に基づいて計算した $H$ に関して、$H$ 桁を超える部分に関してはゼロ埋めしてよいことになっていますが、この方針だと $H$ に満たない桁数しか出てこなさそうです。これは大丈夫なのでしょうか。

Kotlin の REPL、kotlinc で試してみます。

>>> "%.60f".format(0.1)
res0: kotlin.String = 0.100000000000000000000000000000000000000000000000000000000000

まるで 0.1 を正確に表せているかのような出力が出てきますが、かなり misleading な気がします。

>>> "%.60f".format(0.3)
res0: kotlin.String = 0.300000000000000000000000000000000000000000000000000000000000
>>> "%.60f".format(0.1 + 0.2) //       ~v~
res1: kotlin.String = 0.300000000000000040000000000000000000000000000000000000000000

浮動小数点数に怪しいイメージを持っている人が見れば、「0.1, 0.2, 0.3 のような値は正確に表せるが、+ によって誤差が生じた」のような誤解をしてもおかしくなさそうです。

>>> "%a".format(0.3)
res0: kotlin.String = 0x1.3333333333333p-2
>>> "%a".format(0.1 + 0.2)
res1: kotlin.String = 0x1.3333333333334p-2

16 進法の出力もサポートしているものの、「10 進法で正確に保持しているが、16 進法として出力する際に(打ち切り)誤差が生じた」という解釈をされそうな気もします。観測できる挙動からは反駁できない気もします。Double.MAX_VALUE1.7976931348623157E308 が返ってきますが、「10 進法で持っているが、たまたまそういう最大値だった」と言われると困りそうです。

C での仕様

fprintf の変換指定子の f, F の項目には下記の記載しかないように見えます。

The value is rounded to the appropriate number of digits.

IEEE 754 で言うところの convertToDecimalCharacter に従っていると解釈していいと思います。5.12 Details of conversion between floating-point data and external character sequences を参考にするとよいようです。

Conversions to or from decimal formats are correctly rounded.

と書かれています。実際、たとえば fesetround(FE_DOWNWARD) のように丸めモードを変えると printf("%.0f\n", 1.5) の挙動に差異が生じることが確認できます。

roundTiesToEven の説明において「仮数部が偶数の方(状況によってはこれが一意とは限らず、その場合は絶対値が大きい方)」と述べましたが、その状況の一例がここになります。 たとえば 95.0"%.0e" で丸める場合、9e+011e+02 が考えられます。これは仮数部がどちらも奇数なので、roundTiesToEven の最初の tiebreak では定まりません。次の tiebreak によって、絶対値が大きい方である 1e+02 が選ばれます。

四捨五入の話

無限精度においては、正の数 $x$ の四捨五入 $\rounded x$ は $\floor{x+0.5}$ のように書くことができるのは比較的有名です。 では、浮動小数点数 $x$ に対して $\floor{x\oplus\roundp{0.5}}$ (floor(x + 0.5)) と書くのは正当でしょうか。

note: 実際には $\roundp{0.5} = 0.5$ ですが、暗黙にごまかす癖がつくのを避けるために明示しました。

$x_1 = \tfrac12(1-2^{-53})$ や $x_2 = 2^{52}+1$ などが反例として知られています。それぞれ $x_1\oplus 0.5 = 1$、$x_2\oplus 0.5 = 2^{52}+2$ となってしまいます。 自力では思いつけなかった人も、「なんでそんなことになるの?」とはならずに理由を説明できるようになっていたらうれしいです。

指数部を見ながら適切に切り捨てることで $\floor{\bullet}$ の部分は無限精度でできるため、そこの演算では誤差は生じないものとします。 別の方法として、$x \ominus \floor x\ge 0.5$ のとき $\floor x\oplus 1$ を、そうでないとき $\floor x$ を返すのは正当でしょうか。

$x\ge 1$ であれば $\tfrac x2\le\floor x$ ですから、Claim 1 より $x\ominus\floor x = x-\floor x$ が成り立ちます。また、$\floor x=0$ のときも $x\ominus\floor x = x-\floor x$ より、左辺は誤差なしで計算できます。 $x - \floor x$ が整数でないとき、$\floor x\oplus1 = \floor x+1$ となることは示せる($x$ が整数でないことから $x\lt 2^p$ であり、$\floor x\oplus 1$ は誤差なく表せる)ので、$\floor x\oplus 1$ の部分も誤差なしで計算できます。よって、この方法は正当であることがわかります。

また、別の例として、10 進法で書かれた値を小数第ナントカ位まで四捨五入したいようなケースもあります。 たとえば 1.15 を小数第 1 位まで四捨五入しましょう。 $$\roundp{1.15} = {\small 1.149999999999999911}{\scriptsize 18215802998747676610}{\tiny 9466552734375}$$ であるため、これを正確に四捨五入しようとすると、切り捨て方向を選びたくなる気がします。 もちろん 1.1 も正確に表せないため、$1.1$ ではなく $$ \roundp{1.1} = {\small 1.100000000000000088}{\scriptsize 81784197001252323389}{\tiny 0533447265625} $$ が得られていることにも注意が必要です。

たとえば Pythonround(1.15, 1) のように書くと、1.1 が得られます。 一方、Ruby1.15.round 1 と書くと 1.2 が得られます。もちろん 1.2 についても実態は $$ \roundp{1.2} = {\small 1.199999999999999955}{\scriptsize 59107901499373838305}{\tiny 47332763671875} $$ であることに注意しましょう。

Python においては、正確な丸めを用いた変換をできる場合はそれを用いているようです (cpython/Objects/floatobject.c)。 Ruby では $\rounded{\roundp{1.15}\otimes10}\oslash 10$ に基づいているようです (ruby/numeric.c)。

そもそも、10 進法で表したときの桁に基づく操作をしたいときに、2 進法を経由して行おうとするのがよくなさそうです。

下記に諸々のことが載っています。

hnw.hatenablog.com

note: 1.15.round 11.2 になることをもって「Ruby の方が直感的でよい」のような認識をしてしまっている人は、冷静に考え直した方がよいかもしれません。

また、桁数が増えた場合は Ruby (irb 1.11.0 (2023-12-19)) でも次のような挙動を観測できます。特に最後の例では、.round を適用した方が端数が出ているように見えます。

irb(main):001:0> 1.1234567890123455
=> 1.1234567890123455
irb(main):002:0> 1.1234567890123455.round 15
=> 1.123456789012345
irb(main):003:0> 123456789012345.5
=> 123456789012345.5
irb(main):004:0> 123456789012345.5.round 2
=> 123456789012345.52

人間が入力してプログラムに読ませるときは「10 進法から 2 進法へ変換し、有限桁で打ち切る際の誤差」が生じ、コンピュータが人間に読ませる出力をするときは「(読みやすさのために)10 進法の適当な桁で打ち切る差異の誤差」が生じます。 実際に四捨五入の対象となるのは出力時の誤差が生じる前の値なので、混乱したりコンピュータに逆ギレしたりするわけですね。読みやすさのために適当な桁数に丸めるのは、混乱を招くのに大いに貢献している気がします。

二重丸めの話

double での correctly-rounded な値を計算したいときは、long double などの精度の高い型で計算したものを丸めればいい」という迷信はしばしばあります。

四捨五入を複数回行うとまずいという日常的な例 (e.g. 4445 → 4450 → 4500 → 5000 ???) から自然にわかるように、これは当然嘘です。

64-bit 仮数部の long double を用いて考えます。 $x = 2^{64} + 2^{12}$ と $y = 2^{11} - 2^{0}$ とすると、$\roundp{x}_{64} = \roundp{x}_{53} = x$ かつ $\roundp{y}_{64} = \roundp{y}_{53} = y$ であることはわかります。

$p$-bit での $\oplus$ を $\stackrel{p}{\oplus}$ と書くことにすると、 $$ \begin{aligned} x\stackrel{64}{\oplus} y &= x+y+1, \\ x\stackrel{53}{\oplus} y &= x, \\ \roundp{x\stackrel{64}{\oplus} y}_{53} &= x+2^{12} \end{aligned} $$ となり、correctly-rounded な値にはなっていないことがわかります。競プロの文脈で correctly-rounded を厳密に要求する場面は少ない気もしますが、誤差が余分に増えることもあることには注意するとよいかもしれません。

FMA の節で挙げた記事についても見てみるとよいでしょう。

不等式評価や eps の話

何らかの実数 $x$, $y$ に対して $x\lt y$ かどうかを判定したいことはしばしばあります。 それらを浮動小数点数で計算した結果を $\hat x$, $\hat y$ としたとき、$\hat x\lt \hat y$ とすると意図しない結果が得られることがあります。

たとえば、$x = 1+2^{-54}$ と $y = 1+2^{-53}$ のとき、$x\lt y$ ですが、$\roundp{x}=\roundp{y} = 1$ となってしまうため、$\hat x=\hat y$ となります。

>>> 1.0 + 2.0**-54 < 1.0 + 2.0**-53
False

簡単のため、実数 $d$ に対して $d\lt 0$ の判定について考えます。

$\gdef\lomax#1{{#1}_{\text{lo-max}}}$ $\gdef\himin#1{{#1}_{\text{hi-min}}}$

さて、$d\lt 0$ となるケースにおける $\hat d$ の最大値を $\lomax{\hat d}$、$d\ge 0$ となるケースにおける $\hat d$ の最小値を $\himin{\hat d}$ と書くことにします。

このとき、$\lomax{\hat d} \lt \poszero\le \himin{\hat d}$ であれば、$\hat d\lt \poszero$ で判定して問題ないでしょう。 また、$\lomax{\hat d}\ge\himin{\hat d}$ となると、どうしようもなさそうです。アルゴリズムを改善して精度を高めるか、浮動小数点数を介さずにできる方法を考える必要がありそうです。わかりやすい例として、$\lomax{\hat d}=\himin{\hat d}=\poszero$ のときに $\hat d=\poszero$ をどう扱うか考えるとよいかもしれません。

残った $\lomax{\hat d}\lt \himin{\hat d}$ のケースについて考えます。 $\varepsilon\in[\lomax{\hat d}\lldot\himin{\hat d})$ を選んで $\hat d\le \varepsilon$ を判定すればよさそうです。$\hat d$ との判定の等号の有無次第で、$\varepsilon$ を選んでくる区間の開閉は変わります。 また、最初の $\lomax{\hat d}\lt\poszero\le\himin{\hat d}$ は、これの特殊なケースです。

なお、$\hat x\ominus\hat y \lt \varepsilon$ と $\hat x \lt\hat y\oplus\varepsilon$ の真偽は一致するとは限らないことに注意しましょう。$\varepsilon$ を足した際に誤差が生じうるためです。たとえば、$\hat x = 1+2^{-52}$、$\hat y = 1$、$\varepsilon = 5\cdot 2^{-54}$ としてみます。

>>> x = 1.0 + 2.0**-52
>>> y = 1.0
>>> eps = 5.0 * 2.0**-54
>>> x - y < eps
True
>>> x < y + eps
False

結局、$\lomax{\hat d}$ や $\himin{\hat d}$ を見積もることが重要になります。この見積もりの際に、問題ごとの性質を使って ad hoc な考察が必要になることもしばしばあります。「何も考察をせずに $\roundp{10^{-6}}$ などを使って祈る」「テンプレートに eps = 1e-8 と書いておく」というような態度はあまり賢くないように思えます。

INF = 1e9 などについてもそうですが、「たまたまうまくいくことが多いが、完全ではないもの」を何も考えずに使うのは、うまくいかない場面で泥沼にはまる原因になるので、あまりおすすめしないです。使いたい人は勝手に使えばいいと思っているので、別に反論はいらないです。

サイトや文献

  • The Herbie Project, Herbie web demo
    • 数式を入力すると改善してくれるらしい
  • Python 3.12.0 documentation, 15. Floating Point Arithmetic: Issues and Limitations
  • freebsd-src/lib/msun/src
    • 数学関数たちのソフトウェア実装がある。コメントも多数あり。
  • The GNU C Library, 19.7 Known Maximum Errors in Math Functions
  • [Ri66] Rice, John R. “A theory of condition.” SIAM Journal on Numerical Analysis 3, no. 2 (1966): 287-310.
  • [De71] Dekker, Theodorus Jozef. “A floating-point technique for extending the available precision.” Numerische Mathematik 18, no. 3 (1971): 224–242.
  • [Ru88] Siegfried M. Rump “Algorithms for verified inclusions: Theory and practice.” In Reliability in computing, pp. 109–126. Academic Press, 1988.
  • [Go91] Goldberg, David. “What every computer scientist should know about floating-point arithmetic.” ACM computing surveys (CSUR) 23, no. 1 (1991): 5–48.
  • [GKP94] Ronald L. Graham, Donald E. Knuth, Oren Patashnik. Concrete Mathematics: A Foundation for Computer Science, 2nd Edition. Addison-Wesley Professional, 1994.
  • [KD98] William Kahan, Joseph D. Darcy. “How Java’s floating-point hurts everyone everywhere.“ In ACM 1998 workshop on java for high-performance network computing, p. 81. Stanford University, 1998.
  • [Hi02] Higham, Nicholas J. Accuracy and stability of numerical algorithms. Society for industrial and applied mathematics, 2002.
  • [LW02] Loh, Eugene, and G. William Walster. “Rump's example revisited.” Reliable Computing 8, no. 3 (2002): 245–248.
  • [BC04] Bender, Michael A., and Martin Farach-Colton. “The level ancestor problem simplified.” Theoretical Computer Science 321, no. 1 (2004): 5–12.
    • これは浮動小数点関連ではないですが、$\hfloor{\bullet}$ の記法が載っていたものです。
  • [OgRuOi05] Ogita, Takeshi, Siegfried M. Rump, and Shin'ichi Oishi. “Accurate sum and dot product.” SIAM Journal on Scientific Computing 26, no. 6 (2005): 1955–1988.
  • [PiRu13] Jeannerod, Claude-Pierre, and Siegfried M. Rump. “Improved error bounds for inner products in floating-point arithmetic.” SIAM Journal on Matrix Analysis and Applications 34, no. 2 (2013): 338–344.
  • [Ar14] Arnold, Jeff. "Techniques for Floating-Point Arithmetic." (2014).
  • [Kn14] Donald E. Knuth. The Art of Computer Programming, Volume 2: Seminumerical Algorithms, 3rd Edition. Addison-Wesley Professional, 2014.
  • [Std754] “IEEE Standard for Floating-Point Arithmetic,” in IEEE Std 754-2019 (Revision of IEEE 754-2008), vol., no., pp. 1–84, 22 July 2019, doi: 10.1109/IEEESTD.2019.8766229.

関連記事

これまで書いたものたち (👉 浮動小数点数):

浮動小数点数の性質に関して、いくつかの証明をしました。挙動が規格で定められていて、性質を数学的に証明できることがわかってくると、だいぶ仲よしになれるのではないでしょうか。 「(無限精度の演算と比較すると)非自明な挙動をするので、いちいち証明が必要で面倒」といった感覚は否めないかもしれません。

感覚としては、「整数に対して、特定の桁数に丸めてくれる四則演算をする型」のようなものに近く、整数関連の証明が得意な人は力を発揮できるような気がしています。一般性を失わずに指数部を固定して、整数のケースに帰着させるようなことはしばしばあります。

謝辞

今回は何名かの方に review をおねがいしました。ぽかいこちゃん、mod_poppo 先生、hidesugar2 さん、Y.Y. さん、ありがとうございました。

ぽかいこちゃん (MMNMM) は AtCoder の解説 (e.g. ABC 263 Ex) でしばしば浮動小数点数の誤差評価などについて書いており、えびちゃんのここ最近の浮動小数点数シリーズの動機づけの大きい部分です。最近あまり直接お話できていないですが、えびちゃんのことを見てくれているみたいでうれしいです (#3, #5)。

mod_poppo 先生は浮動小数点数オタクの方で、諸々の記事を書いておられます。数年前に書いておられた記事を読んで「浮動小数点数についてちゃんと勉強したいな」と思った記憶があります。具体的にどの記事だったかは記憶が定かではありませんが、本文中で言及した FMA の記事も当時読んでいたと思います。また、えびちゃんはあまりこういったレビュー依頼などが得意ではないのですが、最近書いておられた記事 を読んで、依頼しようと思ったという経緯が実はありました。その結果、人々からためになるコメントがいくつもいただけたのでよかったなぁと思っています。

hidesugar2 さんは Project Euler をたくさんやっていたりして、数学に明るい方です(とえびちゃん的には思っています)。証明でえびちゃんが見落としている部分や、油断して書いていた記法などについても鋭く指摘してくださり、大変ためになりました。この記事に限らず、詰めの甘いツイートなどについてもしばしば指摘してくださったりして、えびちゃんもそうなりたいなぁと思うばかりです。

Y.Y. さんは数値計算というより最適化界隈の方だと思っているのですが、レビューを依頼する流れを(半ば無理やり)作ってしまいました。お忙しい中お読みくださって、わかりにくい箇所の改善の提案や細かいミスの指摘をしてくださいました。以前えびちゃんが 計算量やオーダーに関する記事 を書いたときも Y.Y. さんがお褒めくださったのを覚えています。うれしかったです。

依頼にあたって皆さまが快く引き受けてくださり、とてもうれしかったです。ありがとうございました。

あとがき

めちゃくちゃ長くなってびっくりです。(数式パートも含めて)7 万文字を超えています。すみませんでした。 「下書きを更新する」を押してから更新完了するまでに 15 秒くらいかかっていて、そのたびに長さを実感しました。

たしか 2023 年の 11 月くらいから書いていた記憶はあるのですが、何回か丸っと没にしたりした記憶もあります。

ここ数ヶ月の浮動小数点数シリーズを書きながら、それ以前の自分が知らなかったこと・なんとなくやっていたことをちゃんと考えたり調べたりできてよかったと思います。 皆が皆えびちゃんのように浮動小数点型と仲よくなりたい人ばかりではないと思っているものの、仲よくなりたい人が仲よくなるための一助になれたらうれしいかなという気持ちです。

今回は主に基本的な算術まわりについて調べましたが、基本的な数学関数の correctly-rounded な値を求めるところにも興味が(当然)出てきてしまったので、また調べていきたいです。とはいえ直近は別のことに興味が出てきてしまったので間は空くかもしれません。

また、IEEE 754 に則っていないような処理系が現在どの程度あるのかとかの方面の調査は今回ほぼしておらず、あくまで則っている処理系を前提にしています。AtCoder など競プロ文脈においては問題ないんじゃないかなと思っているのですが、組み込みとかゲームとかの分野ではそうもいかないのかもしれません? 少なくとも、「(IEEE 754 にちゃんと準拠している処理系においてさえ)予測・解析のできない謎挙動をする」という思い込みから脱却できたらいいかなという気持ちです。

ある程度数式の読める競プロ er を想定読者層に(暗に)している気はしており、そうでない人にとっては障壁があるかもしれません。内容が内容でもあるので、そうした人々にリーチする必要があるかは怪しいところがあるかもしれません? 数式なしだとふわっとした理解しかしにくい分野な気がしていて、そういう記事は現状でもたくさんネットにあるような気がします。

さて、えびちゃんが最近やりたいことは「多くの競プロ er がちゃんと知らないまま過ごしている(という偏見がある、えびちゃんのあまり知らない)分野」をちゃんとつぶすことで、浮動小数点数はその手始めという感じでした。次のトピックになりそうと思っているのは簡潔データ構造や組合せゲーム(の特に不偏でないもの)などです。他にはなにがあるでしょう。正直なところ、ABC に出てくる出てこないみたいな基準はあまり興味がありません(数年前は FPS を想定した高度な問題なんて ABC に全然出ませんでしたし)。競プロ er 向け自作平衡 $b$-分木みたいなのも、もしかしたら需要があるかもしれません? いや、需要のあるなしも実はあまり興味がないかもしれません。

それはそれとして、ABC-F くらいで解けていないものも「ちゃんとわかっていない分野」な気がします。はい反省しています。 とはいえこういうのは興味深さ優先探索なので、まぁ好きに生きていけばいいかなとふんわり思っています。こういうところは学生の頃からあまり変わっていない気がします。変わるべきかもしれません?

あとがきに書きたい内容が尽きてきたのでこのくらいで締めにします。次のトピックでお会いいたしましょう。

おわり

おわりです。お疲れさまでした。

*1:ところで浮動小数点数のことを「浮動小数」と呼ぶ人がしばしばおられますが、“浮動” “小数点” “数” なので、あまり妥当な略という感じがしません。Neodym をネオジウムと呼んでしまうのと同様の無頓着さを感じます。小数点を小と略すタイプの人であれば納得します。中文では浮点数と呼ぶらしいです。

*2:automaton の複数形が automata であることとは関係なく DFA (deterministic finite automaton) の複数形が DFAs と呼ばれるの同様、略称の複数形はそういう感じになりそうです。

*3:必ずしも桁数で測る必要はないと思います。

*4:四則演算については可能という話で、一般の数学関数についてはここでは言っていません。

*5:ここを暗黙にやると混乱の元になる気がしています。

*6:UnicodeUTF-8 などが別のレイヤーの話であるのと同様な気がします。

*7:求めた値がどの程度 accurate であるかは、得られた値や precision からはわからず、使ったアルゴリズムを評価する必要があります。

*8:あるいは、有効数字 2 桁で何らかの計算をした結果が 3.0 になったような状況でも、“3.0” というの自体が 2.95 以上 3.05 未満の区間を意味しているわけではないのと同様です。

*9:浮動小数点数はわけのわからない挙動をするという話ではなくて、実際の挙動を正しく認識しましょうという意味です。

*10:仮数部を $1.m_1m_2\lldot m_{p-1}$ のような小数と見るのが冗長になりそうだったので、整数と見るように導入しました。また mantissa の複数形として、mantissas でなく mantissæ, mantissae もあるらしいです。

*11:桁数が収まらない場合は適宜丸められます。$\text{\texttt{0x}}$ の直後の部分は $[1\lldot2)$ に収める必要はありませんが、$\text{\texttt p}$ 以降の部分を省略することはできません。

*12:IEEE 754 が規定される前は、様々なプロセッサなどが思い思いの挙動をしたりして、再現性・可搬性がめちゃくちゃだったらしいです[Go91]

*13:細かい話ではありますが、binary64 でできなかっただけで、コンピュータにとって不可能と言うのはよくないと思います。コンピュータは巡回セールスマン問題を愚直に解くことしかできない、という決めつけの詭弁と似ています。

*14:言語によってはこうした手法がうまくいかないこともありそうです。一旦はうまくいく言語で試すことで我慢します。Kotlin で試したところだと、誤差があるにも関わらず、誤差がないかのような出力をされて困りました。

*15:浮動小数点数は誤差が出るものだから」のような曖昧な解答ではなく、数式での説明を期待しています。

*16:もちろん、$\roundp{0.5} = 0.5$ などが示すように、小数の場合でも誤差が出るとは限りません。

*17:べき乗の部分を先に計算すると想定の挙動と変わるので、おそらくは $5.5b^8$ は $5.5\otimes(b\otimes\dots\otimes b)$ ではなく $5.5\otimes b\otimes\dots\otimes b$ の解釈で合っているのだと思います。

*18:実際には、$\hfloor x = 2^{\emax}$ の際にはこれは該当の型では表せませんが、overflow しない前提なので大きな問題は起きないでしょう。気をつける必要がある状況においては気をつけてください。

*19:$P_{\text{min}}$ などと書きたい気持ちがありますが、規格の記載をそのまま採用しました。

*20:直感的には、下位桁が一致しているというようなことですね。

*21:あまり詳細に踏み込まない一般論で言うなら、そもそも浮動小数点数は誤差を含みがちなので、そうした前提で話されているのかもしれません? 実際にはあまり考えずに言っている人も多そうです。

*22:$\text{\texttt{b*b}}$ ではなく $\text{\texttt{b**2}}$ で計算したところ、$2^{105}+2^{102}+2^{53}$ になりました。環境などによるかもしれません。

*23:$1\ll\delta\ll x$ のようなケースでは困りそうです。

*24:多くの競プロ er が見慣れている $O(f(n) )$ は暗に $n\to\infty$ を前提としていますが、ここでは $x\to0$ なのでやや違和感を抱く人もいるかもしれません。どちらの場合でも、$O(f(\bullet) )$ は「$f(\bullet)$ よりも寄与が小さい項は無視できる」のようなイメージで捉えるとよいと思います。

今年のえびちゃん 2023

今年もやってみます

rsk0315.hatenablog.com

今年書いたもの

カスの記事

カスの記事を上げることでしか消化できない欲求があります。

いま確認したところ、期間限定公開の方は web.archive.orgアーカイブが残っていましたが、公開終了した後のものでした。アーカイブされたのを確認してから公開終了という流れにしてもよかったかなあ(インターネットしぐさ)。

はてなブログ関連

気まぐれアルゴリズム

かわいそうシリーズ

振り返り

DP 関連

メンタルモデルを形成したり、愚直に考える(というか、計算量を一旦忘れて、入出力がなんなのかを数学的に意識する)ことをしたりするのは大事だよな〜と思ったシリーズでした。

数式がちゃがちゃ

このあたりで数学がちゃがちゃをしていたことが、後に浮動小数点数シリーズに足を踏み入れるきっかけになったような気もします。

ブラックボックスのお勉強

assembly を眺めて考えて、むむーとなって楽しかったです。コンパイラ側のコードを読めというのは正しいと思っているものの、二通りある楽しみ方の片方を選んじゃったかなという感じです。

浮動小数点数

最近はここに興味を持っていますね。 入門用記事を書いている途中なのですが、年内はちょっと厳しそうです。 1 月中には上げられたらうれしいですね。あるいは 2024 年内には上がらないかもしれません。

Sqrt Inequality に関しては、ブラックボックスのお勉強シリーズにも含まれるかもしれません。公式解説で言及されていた eps = 1.0E-14 がよくわからなかったので、ちゃんとわかりたい気持ちがあったんですよね。記事を書いた当時、出題から 3.5 年も経っていたことに驚きました。最近始めた人は知らない問題かもしれません。

今年買ったもの

浮動小数点数のお勉強のときに、

などを買ったりしました。前者二つは、8 億円くらいするようなイメージだったのですが、割引コードによって計 1 万円弱で買えて、大変お得でした。最後のはもう少し高かったです。Concrete Math は浮動小数点数関連書籍ではないのですが、前々から気になっていた本だったのと、割引適用のために 2 冊以上買う必要があったので、ちょうどいいと思って買いました。

あと ダウナー系お姉さんに毎日カスの嘘を流し込まれる音声【CV:餅梨あむ】 も書いました。こういうのってもっと高いのかなと思っていました。“““カスの嘘 100 連発””” というタイトルから感じる狂気が好きですね。100 連発を 100 ファイルに分けてシャッフルして、カスの嘘イントロクイズをやりたい気持ちになります。ループすると、#100 の嘘と #1 の嘘でテンションがやや違って面白いです。

他にも、MacBook Pro を買ったり、めがね(かわいい)を買ったりしました。

衝動買いには Mac も含まれています。実はね。

そういえば、【天才博士bot】博士らしい格言31語 日めくりカレンダー も買いました。適切にめくる習慣をつけない場合、特定の日付だけを見る可能性がグッと高まります。不要な月極サービスを解約しやすくなりました*1

来年に向けて

浮動小数点数まわりの話がある程度落ちついたら、FPS や ABC-G 典型をわかっていきたい気持ちがあります。いつまでも解けないままでいるのはつまらないので。 あと ARC 寄りの頭になれるとうれしい気持ちもあるものの、いまいち続かないので困ります。こういうのは(少なくとも、特定の期日までに成長する必要がない限りは)興味深さ優先探索でいい気がしているので、一旦は ABC の方をつぶしたいかなという気持ちです。

  • FPS 全般
    • 話としてはある程度わかっているはずだけど、いまいち自分の道具になっていない感がある
    • SPS ってなに?1, 2
  • stack を使って高速化できる系のやつ
  • フローに帰着できるやつ
    • というよりは最小カットとかを学ぶべきそう
  • lowlink
    • 5 年前くらいから知ってはいるのに、まだお友だちじゃない
  • 高速何々変換シリーズ
    • これもわかってなさの方向性としては FPS に近いか
  • 非不偏ゲーム
    • e.g. ABC 229 H
    • ABC 卒業的な意味で優先度は高くない気がするものの、こういうのに惹かれちゃう

ABC-{G, H, Ex} くらいに関しては、解説をざっと一通り斜め読みしてジャンル分けして、知りたいところをつまめるように索引っぽくしておくとよさそうな気もしています。それはそれとして ABC-F も埋めなきゃだと思うので、G 以降のお勉強が飽きたら交互に、って感じにしようかなあ。

結局のところ、「レートを上げるためのお勉強」という感じではなくて、自分が知らないところをつぶしたいという感じに近そうです。自分が複数人いて複数のトピックを並行してつぶせたらうれしいのになあという気持ちです。 そういう意味合いだと、学校での授業は、複数の科目を時間ごとに区切って並行して学べるのでうれしいというような気もします。

労働している場合ではなくて、隠居して勉強だけしてたくない?という身勝手な欲求も生えそうです。 「社会人になると時間が取れないというのは甘えで、捻出できる時間はいくらでもあるでしょ」のような指摘は正しい側面もあるとは思うんですが、的外れでもある気はします。そもそもいくらでもあるというのが嘘なので。

競プロ以外

うれしい

えびちゃんが Twitter を始めてすぐ(11 年前くらい)からのフォロヮと初めてオフで会いました。 会いたい人とは会えるうちに会っておいた方がいいですね。 一部のフォロヮやインターネットストーキングが得意な人は写真を見たかもしれません。

え、11 年も Twitter やってたんですか、こわすぎる。

おいしい

何年か前によく行っていたところ(生活圏内から遠め)と同じ系列のお店が、実は生活圏内に存在していたことがわかったので、最近行きました。

かわいい

環お姉ちゃんに沼りそうです。水篠みずしの たまき ちゃんです。一人称が「環ちゃん」なの、えびちゃんと親近感がありますね。

M.L.V.G は環お姉ちゃんの曲ではないですね。環お姉ちゃんの曲は出るのでしょうか。

M.L.V.G がテトラの子たちのペンライトカラーの頭文字 (midnight blue, lemon chiffon, vermilion, gainsboro) になっていることと、譜面が文字押し(押し?)になっていることに最近気づきました。文字押しと言えば、jubeat の Qubellic Prism ラストの TYFP (thank you for playing) を最近見て、いいなーってなりました(オタクはそういうのが好きなので)。

かっこいい

3b1b/manim のアニメーション、かっこいいです。 あと 3B1B の日本語翻訳版のナレーションの人の声や話し方がめちゃ好きです。

たのしい

wafflegame を続けています。5/5 を狙うようになってから別ゲーと化したような気がします。royale や deluxe はむずかしいです。

「J が入っていたら ninja, emoji になりがち。たまに eject もある」などの典型があります。ところで、🥷 は ninja でもあり emoji でもあって面白いです。

5/5 でクリアするのは、(答えがわかっている前提なら)下記を $1\le A_i\le26$ で解く(復元つき)のに帰着できる気がします。$N$ の上限は $16$ (daily) だったり $26$ (deluxe) だったりです。royale はよくわかりませんが $21$ 前後くらいでしょうか。

最適行動をしたときのムーブ数が既知情報なので、それも使えるのかもしれません?

atcoder.jp

おわり

来年もよろしくおねがいします。

*1:買った人にだけしか伝わらない気がしますが、まぁよいでしょう。

ポインタ木なしで償却 O(log(n)²) 時間の predecessor query

お風呂で思いついたシリーズです。

できる操作は次の通りです。

  • $\texttt{new}()$
    • $∅$ で初期化する
  • $S.\texttt{insert}(x)$
    • $S ← S ∪ \{x\}$ で更新する
  • $S.\texttt{remove}(x)$
    • $S ← S \smallsetminus \{x\}$ で更新する
  • $S.\texttt{less}(a)$
    • $\max {\{x∈S \mid x\lt a\}}$ を返す
  • $S.\texttt{lesseq}(a)$
    • $\max {\{x∈S \mid x\le a\}}$ を返す

これらについて、償却 $O(\log(n)^2)$ 時間で処理します。$n$ はその時点のクエリ数だったり、$|S|$ に改善できたりします。

基本方針

静的データ構造を $\log(n)$ 個持つことで動的にするテクです。

区間最小値を求めるセグ木を $\log(n)$ 個持ちます。 $i$ 番目のセグ木は、長さが $0$ または $2^i$ のいずれかです。

$S.\texttt{insert}(x)$ の際は、長さが $0$ となるセグ木のうち番号が最小のものが $i$ 番だとして、$[0\ldots i)$ 番目のセグ木の要素と $x$ を昇順に連結させたセグ木を $i$ 番目に置き、$[0\ldots i)$ 番目のセグ木は空にします。

$S.\texttt{remove}(x)$ の際は、セグ木上の二分探索を行うことで $x$ の検索を行うことができるため、見つかった場合は $∞$ で更新します。各セグ木について検索を行います。

$S.\texttt{less}(a)$ や $S.\texttt{lesseq}(a)$ の際は、各セグ木に対して二分探索を行い、条件を満たす最大の元(空の場合は無視)を集め、全体の最大値を求めます。

細部

有限値だけ見ると昇順に並んでいますが、$∞$ が入ることもあるので、気をつけて場合分けをします。

挿入の際、重複して要素を入れないように注意しましょう。入れたい場合も注意します。

$[0\ldots i)$ 番目のセグ木と $\{x\}$ を組み合わせるパートでは、マージソートで使うサブルーチンの要領で、ソート済みの列二つをソート済みになるように並べます。等比数列の和の形になるので、$O(2^i)$ 時間で可能です。 $i$ 番目のセグ木を作るためには $2^i$ 回の挿入が必要なため、セグ木の部分に関しては償却 $O(1)$ 時間です。 挿入操作の前に重複の確認が入り、$O(\log(n)^2)$ 時間です。

削除の際、$i$ 番目のセグ木に $∞$ が $2^{i-1}$ 個以上できた場合、$i-1$ 番目のセグ木の長さが $0$ であれば降格させるか、そうでない場合は統合させる($i-1$ 番目のセグ木の要素とくっつけて $i$ 番目のセグ木を再構築する)こともできます。これにより、全体のうちの $∞$ が半分未満になるようにでき、管理する要素が $2|S|$ 程度で抑えられます。

successor query に対応させる場合、同じものを逆向きの比較関数で構築してもよいですし、隣り合う要素の添字を管理するなどで対応してもよいです。

実装

judge.yosupo.jp

すいません、Rust です。 本当は Python のことを考えて思いついたのですが、Python で実装して速くなる気がしなかったので実装していません。有志の方にお願いします。 案の定 Rust でも十分遅かったです。

普通に TLE すると思っていたんですが、9170 ms / 10 s で間に合ったのでうれしかったです。

その他

値の削除と predecessor query と要素の列挙が高速にできるデータ構造であれば、セグ木の部分は置き換えられそうです。

また、今回は各セグ木での答えの max を考えましたが、このようにいくつかに分割したものでの答えを組み合わせられるタイプのクエリであれば、別の種類のクエリにも答えられそうです。

例:Aho–Corasick のやつ

rsk0315.hatenablog.com

数年前に書いた記事ですが、手法はこれと同じです。他にも、binomial heap も似たような感じで $2^i$ サイズごとに分けたりしますね。

あとがき

ふとしたときに思いついた log 倍悪い系のデータ構造でした。他には、区間を管理するデータ構造をビットごとに持つことで RUQ するやつが(自分の中では)印象深いです。

rsk0315.hatenablog.com

これはお皿洗いしているときに思いついたシリーズでした。

浮動小数点数の記事はまだ時間がかかりそうです、すみません。

おわり

こういうふわっとした記事もいいですね。

xx.yy を 100 倍したら xxyy にできる?

小数点以下が高々 2 桁まである値が十進表記で与えられたとき、「それを浮動小数点数型に parse してから 100 倍する」ことで、元の値の小数点を 2 つぶん動かしたものにできますか?という話です。

>>> print(f'{float("10.00")*100:.100g}')
1000
>>> print(f'{float("10.01")*100:.100g}')
1001
>>> print(f'{float("10.1")*100:.100g}')
1010

お、いけそう?

>>> print(f'{float("10.03")*100:.100g}')
1002.9999999999998863131622783839702606201171875
>>> print(f'{float("10.05")*100:.100g}')
1005.0000000000001136868377216160297393798828125

すみません、失敗してしまいました。ということで不可能ということがわかりました。いかがでしたか?


もちろん、これで終わってはつまらないので、もう少し有益な話をします。たとえば、次のような方針たちが自然に思いつくと思います。

  • 100 倍してから、ある微小値を足し、小数部分を切り捨てる
  • ある微小値を足してから 100 倍し、小数部分を切り捨てる
  • 100 倍した後、最も近い整数に丸める

これらの正当性について考えていきたいです。100 倍に限らず、実際に出題された例に応じた桁数を考えます。

導入

整数型以外を好まない皆さまにおかれましては、

>>> (lambda x, y: 100 * x + y)(*map(int, input().split('.')))  # 10.03
1003

int ipart, fpart;
scanf("%d.%d", &ipart, &fpart); // 10.03
int x = ipart * 100 + fpart; // 1003

のように、「小数点で分割してから整数として parse し、ナントカ倍して足す」という方針を取ることが多いでしょう。 Python においては int('10.05'.replace('.', '')) のように小数点を消す方針でやることもありますが、桁数が固定でない場合は厄介です。

もちろんこの方法が悪いとは思いませんし、競技プログラミングにおいて初心者ができるようになるべき方針は、このような整数型を用いる方法だと思います。

藍染惣右介「整数型とは 整数型に縋らねば生きて行けぬ者の為にあるのだ」

BLEACH #47, p. 50 (407. DEICIDE9)

浮動小数点数型でやった場合でも正当性を示せる人の方がかっこいいと思うので、そうなりたいです。

前提・記法

下記を先に読むとよいかもしれません。証明などは詳しく追わなくても大丈夫だと思います。

rsk0315.hatenablog.com

ここで定義されている $\roundp x$, $\hfloor x$, $x\otimes y$ に関して、同様のものを用います。また、$x\oplus y \triangleq \roundp{x+y}$ も用います。

例題

両方向の反例(少し小さくなるケース・少し大きくなるケース)もそれぞれ挙げておきます。

# 1.2
>>> print(f'{1.13*100:.100g}')
112.9999999999999857891452847979962825775146484375
>>> print(f'{1.09*100:.100g}')
109.0000000000000142108547152020037174224853515625

# 2.3
>>> print(f'{16.002*1000:.100g}')
16001.999999999998181010596454143524169921875
>>> print(f'{16.001*1000:.100g}')
16001.000000000001818989403545856475830078125

# 5.4
>>> print(f'{10000.0009*10000:.100g}')
100000008.99999998509883880615234375
>>> print(f'{10000.0004*10000:.100g}')
100000004.00000001490116119384765625

# 4.9
>>> print(f'{1024.000000006*1e9:.100g}')
1024000000005.9998779296875
>>> print(f'{1024.000000011*1e9:.100g}')
1024000000011.0001220703125

これらの例は binary64 の形式で roundTiesToEven による丸めを前提としています(いわゆる C++ で言う double のデフォルトの挙動)が、long double__float128 などを使ったり、別の丸めモードにした場合でも反例は存在します*1

なお、x.ye9 のようなリテラルの形だと、丸めが $\roundp{x.y}\otimes 10^9$ ではなく $\roundp{x.y\times 10^9}$ のように行われるため、結果が異なる場合があることに注意してください。 今回の問題設定では「文字列を一度 parse してから数倍する」であり、一旦 parse した時点で丸めが生じています。

>>> print(f'{1024.000000006*1e9:.100g}')
1024000000005.9998779296875
>>> print(f'{1024.000000006e9:.100g}')
1024000000006

考察

上で述べた方針たちについて考えていきます。

方針 1

Claim 1: $p=53$ のとき、$\delta = \roundp{0.003}$ とすると、任意の $0\le m\le 10^{13}$ に対して $$\floor{\roundp{m\times 10^{-k}} \otimes 10^k\oplus \delta} = m$$ が成り立つ。ただし、$k$ は overflow・underflow が起きない範囲で任意の非負整数とする。

Proof

$\roundp{m\times 10^{-k}} = m\times 10^{-k}+\varepsilon$ と書ける。ここで、$|\varepsilon|\le \hfloor{m\times 10^{-k}}\cdot 2^{-p}$ を満たす。

$$ \begin{aligned} \roundp{m\times 10^{-k}}\times 10^k &= (m\times 10^{-k}+\varepsilon)\times 10^k \\ &= m + 10^k\varepsilon \end{aligned} $$ であり、 $$ \begin{aligned} |10^k\varepsilon| &\le 10^k\cdot\hfloor{m\times 10^{-k}}\cdot 2^{-p} \\ &\le 10^k\cdot(m\times 10^{-k})\cdot 2^{-p} \\ &= m\cdot 2^{-p}. \end{aligned} $$ すなわち、 $$ m(1-2^{-p})\le\roundp{m\times10^{-k}}\times10^k\le m(1+2^{-p}) $$ であり、 $$ \begin{aligned} \roundp{m\times10^{-k}}\otimes10^k &\ge m(1-2^{-p})-\hfloor{m(1-2^{-p})}\cdot 2^{-p} \\ &\ge m(1-2^{-p})-m(1-2^{-p})\cdot 2^{-p} \\ &= m(1-2^{-p})^2 \\ %&= m(1-2^{1-p}+2^{-2p}) \\ &= m - m\cdot(2^{1-p}-2^{-p}). \end{aligned} $$

$m\le10^{13}$, $p=53$ とすると、 $$ \begin{aligned} \roundp{m\times10^{-k}}\otimes10^k &\ge m-10^{13}\cdot (2^{-52}-2^{-106}) \\ &\gt m-10^{13}\cdot 2^{-52} \\ &\gt m - 0.0022205 \end{aligned} $$ となる。上方向からも、同様に $$ \roundp{m\times10^{-k}}\otimes10^k \lt m + 10^{13}\cdot 2^{-52} $$ と抑えられることがわかる。

そこで、簡単のため $\delta = \roundp{0.003}$ とし、 $$\roundp{m\times10^{-k}}\otimes10^k\oplus\roundp{0.003}\lt m+1$$ を示せばよい。丸めの範囲を考えて $$\roundp{m\times10^{-k}}\otimes10^k+\roundp{0.003}\lt (m+1)\cdot (1-2^{-53})$$ を示す。 $$ \begin{aligned} &\phantom{{}={}} (m+1)\cdot(1-2^{-53}) - \left(\roundp{m\times10^{-k}}\otimes10^k\right) \\ &\gt (m+1)\cdot(1-2^{-53}) - (m + 10^{13}\cdot 2^{-52}) \\ %&= m\cdot(1-2^{-53}) + (1-2^{-53}) - m - 10^{13}\cdot 2^{-52} \\ %&= -2^{-53}m + (1-2^{-53}) - 10^{13}\cdot 2^{-52} \\ &\ge -10^{13}\cdot2^{-53} + (1-2^{-53}) - 10^{13}\cdot 2^{-52} \\ &\gt 0.9 \\ &\gt \roundp{0.003}. \quad\qed \end{aligned} $$

note: $\roundp{0.003} = \tfrac3{1000}\cdot(2^{61}+48)\times 2^{-61} \gt 0.003$ です。

方針 2

Claim 2: $p=53$ のとき、$\delta = \hfloor{\roundp{m\times 10^{-k}}}\times 2^{1-p}$ とすると、任意の $0\le m\le 10^{13}$ に対して $$\floor{(\roundp{m\times 10^{-k}} + \delta) \otimes 10^k} = m$$ が成り立つ。ただし、$k$ は overflow・underflow が起きない範囲で任意の非負整数とする。

なお、$\roundp{m\times 10^{-k}} + \delta$ は $\roundp{m\times 10^{-k}}$ より大きい最小の浮動小数点数であり、nextafter()next_up() などの関数を用いることで得ることができる。

Proof

$\roundp{m\times 10^{-k}} = m\times 10^{-k}+\varepsilon$ と書ける。ここで、$|\varepsilon|\le \hfloor{m\times 10^{-k}}\cdot 2^{-p}$ を満たす。 よって、 $$ \begin{aligned} (\roundp{m\times 10^{-k}}+\delta)\times 10^k &\le (m\times 10^{-k}+\delta+\varepsilon)\times 10^k \\ &= m + (\delta+\varepsilon)\times 10^k \end{aligned} $$ となる。ここで、 $$ \begin{aligned} \delta+\varepsilon &\ge \hfloor{\roundp{m\times10^{-k}}}\times 2^{1-p} - \hfloor{m\times 10^{-k}}\times 2^{-p} \\ &\ge \hfloor{m\times10^{-k}}\times 2^{-p} - \hfloor{m\times 10^{-k}}\times 2^{-p} \\ &\ge 0 \end{aligned} $$ より、$(\roundp{m\times 10^{-k}} + \delta) \otimes 10^k\ge m$ は成り立つ。また、 $$ \begin{aligned} \delta+\varepsilon &\le \hfloor{\roundp{m\times10^{-k}}}\times 2^{1-p} + \hfloor{m\times 10^{-k}}\times 2^{-p} \\ &\le \hfloor{m\times10^{-k}}\times 2^{2-p} + \hfloor{m\times 10^{-k}}\times 2^{-p} \\ &\le (m\times10^{-k})\times 2^{2-p} + (m\times 10^{-k})\times 2^{-p} \\ &\le 5(m\times 10^{-k})\times 2^{-p} \\ \end{aligned} $$ より、 $$ \begin{aligned} m + (\delta+\varepsilon)\times 10^k &\le m + 5m\cdot 2^{-p} \\ &\le m + 5\cdot 10^{13}\cdot 2^{-53} \\ &\lt m + 0.005552 \\ \end{aligned} $$ となる。方針 1 同様、これは $(m+1)\cdot(1-2^{-53})$ より小さいことが示せるので、$(\roundp{m\times 10^{-k}} + \delta) \otimes 10^k\lt m+1$ となる。$\qed$

方針 3

Claim 3: $p=53$ のとき、任意の $0\le m\le 10^{13}$ に対して $$\rounded{\roundp{m\times 10^{-k}}\otimes 10^k} = m$$ が成り立つ。ただし、$k$ は overflow・underflow が起きない範囲で任意の非負整数とする。

なお、実数 $x$ に対して、$x$ に最も近い整数を $\rounded x$ によって表す。round() などの関数で得られる。 ここでは、最も近い整数が複数あるケースについては考えなくてよいことが示せる。

Proof

Claim 1 の証明の中で、該当の範囲において $$ m - 0.0022205 \lt \roundp{m\times10^{-k}}\otimes10^k\lt m + 0.0022205 $$ が成り立つことを示した。これで十分である。$\qed$

方針 4

$10^k$ を掛けた後に nextafter() を行い、切り捨てる方針ではどうでしょう。証明・反証は読者への課題とします。 気が向いたら追記するかもしれません。

関連文献

qiita.com

あとがき

競プロの文脈では $10^{13}$ 個の全探索は絶望的ですが、こうした文脈で多少の時間を掛けてもよい状況であれば $10^{13}$ 個は全然いける寄りの個数なので、一応ちゃんと全てテストしています。特に、浮動小数点数型に関する簡単な演算であれば比較的軽い気がします。

方針 1–3 それぞれについて、1 秒あたり $10^9$ 回程度のテストができています。下記がテスト全体・ケースあたりにかかった時間です (Rust)。

  • 方針 1 → 2:33:16 (0.9196 ns per case)
  • 方針 2 → 4:00:15 (1.4415 ns per case)
  • 方針 3 → 2:29:34 (0.8974 ns per case)

Rust に限らず C++ でも試してみましたが、.next_up() std::nextafter() は重いのかな?という感じでした。

AtCoder 環境でも $2\times 10^9$ ケース回したところ、方針 1 が 2359 ms、方針 2 が 8771 ms、方針 3 が 4760 ms (std::round()), 6893 ms (std::lround()) でした。プロセッサだったり最適化のかかり方などにもよるのかもしれません。そこの調査に関してはこの記事の範囲外という気がしています。

ともかく、こうした見積もり方がわかってくると、競プロ er あるあるであるところの

eps ってよくわかんないけど、強い人が 1e-6 とかにしてるし、とりあえず決め打ちでテンプレに #define eps 1e-6 って書いてる

のようなものが愚かに感じられるのではないでしょうか(この記事を読んで「どうやら eps3e-3 で決め打ちすればいいらしい」と解釈する人がいたら困ります。考えて読んでほしい)。 にぶたんの上限を常に 1e9 とか 1e18 とかで決め打ちするのが愚かだと思わないタイプの人とは思想が合わないと思うので、いちいち反論してくれなくても大丈夫です。

また、冒頭ではいくつかの成功する例も挙げましたが、浮動小数点数あるあるの「いくつかのケースではなんかうまくいく」に対する注意喚起でもあります。 貪欲法と並び、浮動小数点数型の演算も「未証明の方法は大抵誤っている」側だと思うので、証明の方針めいたものすらない状態で「いくつかテストが通ったからよさそう」というのは控えた方がよいです。

さて、一旦書こうと思っていた問題提起系の紹介記事は済ませたので、次は基礎の話を書こうかなと思います。少し期間が空くかもしれません。

これまで書いたものたち (👉 浮動小数点数):

rsk0315.hatenablog.com

rsk0315.hatenablog.com

rsk0315.hatenablog.com

rsk0315.hatenablog.com

おわり

お疲れさまでした。そろそろ Advent Calendar の時期ですね。架空 Advent Calendar 2023 の -18 日目の記事でした。

*1:入力が非負として、roundTowardPositive を使いつつ切り捨てればよいのでは?と思いつつ、それができる言語ばかりではないですね。

u64 の平方根を f64 で計算するやつに関して

符号なしの 64-bit 整数型の値に関して平方根を計算したい局面はそこそこあります。 しかし、多くの言語ではそのための関数が用意されておらず、浮動小数点数型用のメソッドを使う人が多くいます。

今回は、そうした方針が妥当か?(誤差が出うることを踏まえつつ、正しい値を復元するのは簡単か?)という話です。

記法・前提

浮動小数点数型に対する一般的な知識はあるとします。ない人向けの記事はそのうち書きます。

実数 $x$ を該当の浮動小数点数型の値に丸めたものを $\roundp x$ と書くことにします。

IEEE 754 の規格において、sqrt(x) は $\roundp{\sqrt x\,}$ を返すことを規定しているので、ここではそれを仮定します。

note: 他の関数の誤差に関して

GCC において、sin(x), log(x), cbrt(x) などは、$\roundp{\sin(x)}$, $\roundp{\ln(x)}$, $\roundp{\sqrt[3]{x}\,}$ を返すとは限りません*1sqrt などは IEEE 754 の規格で required の枠に含まれるのに対し、sin などは recommended の枠に含まれるので、建前上は「recommended の方は提供していませんよ〜」ということなのでしょうか。今後、正確に丸めた (correctly-rounded) $\roundp{\sin(x)}$ を返す crsin(x) のような関数を提供するかも?とのことです。

www.gnu.org

任意の正規化数 $x$ に対し、$2^e \le x\lt 2^{e+1}$ が成り立つような $e$ が一意に存在します。この $e$ を用いて、$\hfloor x=2^e$ と定義します。 $\hfloor x\le x\lt 2\hfloor x$ や $\tfrac x2\lt\hfloor x\le x$ が成り立ちます。 また、$x$ が正整数の下では $\hfloor x\le x\le 2\hfloor x-1$ や $\tfrac{x+1}2\le\hfloor x\le x$ が成り立ちます。

下記では、次のような関数を作ることを想定します。

/// $\floor{\sqrt n}$ を返す。
fn uint_floor_sqrt(u64 n) -> u64 {
    f64 tmp = (n as f64).sqrt();
    // tmp が持ちうる誤差に対して何かしらの補正をする
    todo!()
}

どのように補正を行うかについてはまだ与えないことにします(適当に考えた方針が「実はよかったよね」とするのではなく、考察ベースで自然なアルゴリズムを考えたいため)。考察次第では、.sqrt() をする前に何かしらをするかもしれません。

??「仮数部が 64 bits 以上の long double がある環境なら long double でやればよくないか?」

それはそう... そうした型のない言語を想定しているとかの建前でやり過ごします。

考察

note: $n \le 2^{53}$ では $\roundp n=n$ ですが、それより大きい場合はそうとは限らないため、気をつける必要があります。

$\gdef\roundsqrt#1{\roundp{\sqrt{#1}\,}}$

Observation 1:

u32 に関してはすぐに全探索できるため、その範囲を f32仮数部 24 bits)で計算してみます。

fn main() {
    for i in 0..=u16::MAX as u32 {
        for j in i * i..=i * i + 2 * i {
            let r = (j as f32).sqrt() as u32;
            assert_eq!(r, i, "failed on {}", j);
        }
    }
    println!("passed all tests.");
}
assertion `left == right` failed: failed on 16785407
  left: 4097
 right: 4096

やはり、これはどうやらうまくいかないようです。 $$ \floor{\roundsqrt{(2^{12}+1)^2-2}} = 2^{12}+1. $$

note: f64 においても、$\floor{\sqrt n}$ を固定して二分探索を使ったり、固定した際の最大値を試してから降順に見ていくなどで、同様の探索ができます。 最小の反例は 4503599761588224 で、下記のようになっています。 $$ \floor{\roundsqrt{(2^{26}+1)^2-1}} = 2^{26}+1. $$ $(2^{\floor{p/2}}+1)^2$ より少し小さいあたりでだめになってしまう感じでしょうか。

Observation 2:

では、$+1$ の誤差は許容してみましょう。

fn main() {
    for i in 0..=u16::MAX as u32 {
        for j in i * i..=i * i + 2 * i {
            let r = (j as f32).sqrt() as u32;
            assert!([i, i + 1].contains(&r), "failed on {}", j);
        }
    }
    println!("passed all tests.");
}
passed all tests.

これはよいらしいです。

Observation 3:

ということで、u64 にして次のような関数を考えてみます。

fn u64_floor_sqrt(n: u64) -> u64 {
    let tmp = (n as f64).sqrt() as u64;
    let tmp_m1 = tmp.saturating.sub(1); // if tmp >= 1 { tmp - 1 } else { tmp };

    // ((tmp - 1) + 1)**2 <= n, in an overflow-safe way
    if tmp_m1 * (tmp_m1 + 2) < n { tmp } else { tmp_m1 }
}

$n = 0$ の場合や、$n = 2^{64}-1$ で $2^{32}$ が返ってきた場合などにオーバーフローしないように気をつけています。 そうしたケースに気を払う必要がない状況では、もう少し雑にやっても大丈夫そうです。

fn main() {
    for i in 0..=u32::MAX as u64 {
        let actual_left = u64_floor_sqrt(i * i);
        let actual_right = u64_floor_sqrt(i * i + 2 * i);
        let expected = i;
        assert_eq!(actual_left, expected);
        assert_eq!(actual_right, expected);
    }
    println!("passed all tests.");
}
passed all tests.

これもよいらしいです。$\floor{\sqrt n}-1$ や $\floor{\sqrt n}+2$ になることはなさそうです。

$\roundp x$ や $\sqrt x$ の単調性から、u64_floor_sqrt は正当そうです。 これで「u64 の範囲ではこれでいけます、めでたし」で終わってもいいにはいいんですが、たとえば C++ において unsigned __int128 の範囲を __float128 で計算するのを同じ方針で調べることはできないため、もう少し数学的に考えます。

もしかすると 128-bit では反例が見つかるかもしれないし、そうではないかもしれませんから。

ここまでの話をまとめると、次のようなことを示したいです。

Conjecture 4: それなりに大きい $n$ に対しても $\floor{\sqrt n} \le \roundsqrt{\roundp n} \lt \floor{\sqrt n}+2$ が成り立つ。

まず、浮動小数点数 $\hat x$ と実数 $x$ に対して、$\roundp x = \hat x$ となる十分条件について触れておきます。 $$ \roundp x=\hat x \impliedby \begin{cases} \hat x-\hat x\cdot 2^{-1-p} \le x\le \hat x+\hat x\cdot 2^{-p}, & \text{if }\hfloor{\hat x}=\hat x; \\ \hat x-\hfloor{\hat x}\cdot 2^{-p}\lt x\lt\hat x+\hfloor{\hat x}\cdot 2^{-p}, & \text{if }\hfloor{\hat x}\gt\hat x. \end{cases} $$ $\roundp x\ge \hat x$ などを示したいときは片側の不等式が成り立てばよいです。

$0\le k\le 2^{p-1}$ なる整数 $k$ に対して、整数 $n\in[k^2\ldots(k+1)^2)$ が conjecture 4 を満たすことを示します。

整数性から $n\in[k^2\ldots k(k+2)]$ であることと、単調性から $$\roundsqrt{\roundp{k^2}}\le\roundsqrt{\roundp n}\le \roundsqrt{\roundp{k(k+2)}}$$ であることに注意します。そのため、 $$ k\le \roundsqrt{\roundp{k^2}} \wedge k+2\gt \roundsqrt{\roundp{k(k+2)}} $$ を示せば十分です。考えている状況に合わせ、$\roundp{k^2}\lt\infty$ である前提とします。 また、この範囲で $\roundp k=k$ が成り立つことに注意しておきます。 加えて、$p\ge 3$ も仮定しておくことにします。

Lemma 5: $0\le k\le 2^p$ なる整数 $k$ に対して、$k\le \roundsqrt{\roundp{k^2}}$ が成り立つ。

disclaimer: もともと $k\le 2^p$ で示す予定だったのですが、Lemma 6 が $k\le 2^{p-1}$ でしか成り立たなさそうということが Lemma 5 を示した後にわかりました。Lemma 5 を直すのが面倒だったのでそのままになっています。

Proof

ある $e$ が存在して $k = 2^e$ のとき、$2^{2e} = 2^{p-1}\times 2^{2e+1-p}$ と表せるため、$\roundp{2^{2e}}=2^{2e}$ となる。 $\roundsqrt{\roundp{k^2}} = \roundsqrt{k^2} = \roundp k=k$ が成り立つ。すなわち、$k\le \roundsqrt{\roundp{k^2}}$ となる。

$k = 0$ のときも同様にして成り立つことが示せる。

次に、$k = 2^e$ なる $e$ が存在しないかつ $k\ge1$ とし、 $$ k-\hfloor{k}\cdot 2^{-p}\lt \sqrt{\roundp{k^2}} $$ を示す。これらを評価すると、 $$ \begin{aligned} k-\hfloor k\cdot 2^{-p} &\le k-\tfrac{k+1}2\cdot 2^{-p} \\ &= (1-2^{-1-p})k - 2^{-1-p}; \\ \sqrt{\roundp{k^2}} &\ge \sqrt{k^2-\hfloor{k^2}\cdot 2^{-p}} \\ &\ge \sqrt{k^2-k^2\cdot 2^{-p}} \\ &= k\cdot \sqrt{1-2^{-p}} \end{aligned} $$ となるので、$(1-2^{-1-p})k-2^{-1-p}\lt k\cdot\sqrt{1-2^{-p}}$ を示せばよい。非負であることは明らかであるから二乗の差を考え、 $$ \begin{aligned} &\phantom{{}={}} (k\cdot\sqrt{1-2^{-p}})^2 - ( (1-2^{-1-p})k-2^{-1-p})^2 \\ &= (1-2^{-p})k^2 - ( (1-2^{-1-p})k-2^{-1-p})^2 \\ %&= (1-2^{-p})k^2 - ( (1-2^{-1-p})^2k^2-2^{-p}\cdot(1-2^{-1-p})k+2^{-2-2p}) \\ %&= (1-2^{-p})k^2 - ( (1-2^{-p}+2^{-2-2p})k^2-2^{-p}\cdot (1-2^{-1-p})k+2^{-2-2p}) \\ %&= -(2^{-2-2p}\cdot k^2-2^{-p}\cdot (1-2^{-1-p})k+2^{-2-2p}) \\ %&= -2^{-2-2p}\cdot k^2+2^{-p}\cdot (1-2^{-1-p})k-2^{-2-2p} \\ %&= 2^{-p}\cdot (-2^{-2-p}\cdot k^2 + (1-2^{-1-p})k - 2^{-2-p}) \\ %&= 2^{-2p}\cdot (-2^{-2}\cdot k^2 + (2^p-2^{-1})k - 2^{-2}) \\ &= 2^{-2-2p}\cdot (-k^2 + (2^{2+p}-2)k - 1) \\ \end{aligned} $$ となる。

$-k^2+(2^{2+p}-2)k-1$ の部分が、主要項が負の $2$ 次関数になっていることに注意すると、これが $k=1$ および $k=2^p$ で正になっていることを示せば、$1\le k\le 2^p$ の範囲で正であることがわかる。 $$ \begin{aligned} -1^2 + (2^{2+p}-2)\cdot 1-1 &= 4\cdot(2^p-1); \\ -(2^p)^2 + (2^{2+p}-2)\cdot 2^p-1 %&= -2^{2p} + 2^{2+2p} - 2^{1+p} - 1 &= 2^p\cdot(3\cdot 2^p-2)-1 \end{aligned} $$ より、$p\ge 1$ でこれらは正となることがわかる。

以上より、$0\le k\le 2^p$ なる任意の整数 $k$ に対して $k\le\roundsqrt{\roundp{k^2}}$ が成り立つ。$\qed$

Lemma 6: $0\le k\le 2^{p-1}$ なる整数 $k$ に対して、$k+2\gt \roundsqrt{\roundp{k(k+2)}}$ が成り立つ。

Proof

任意の $0\le k\le 2^{p-1}$ に対して、 $$ \sqrt{\roundp{k(k+2)}} \lt (k+2)-\hfloor{k+2}\cdot 2^{-p} $$ を示せばよい。両辺を評価して $$ \begin{aligned} \sqrt{\roundp{k(k+2)}} &\le \sqrt{k(k+2)+\hfloor{k(k+2)}\cdot 2^{-p}} \\ &\le \sqrt{k(k+2)+k(k+2)\cdot 2^{-p}} \\ &= \sqrt{k(k+2)}\cdot \sqrt{1+2^{-p}}; \\ (k+2)-\hfloor{k+2}\cdot 2^{-p} &\ge (k+2)-(k+2)\cdot 2^{-p} \\ &= (k+2)\cdot (1-2^{-p}) \end{aligned} $$ となるので、 $$ \sqrt{k(k+2)}\cdot\sqrt{1+2^{-p}} \lt (k+2)\cdot(1-2^{-p}) $$ を示す。Lemma 5 同様にして、 $$ \begin{aligned} &\phantom{{}={}} ( (k+2)\cdot(1-2^{-p}))^2 - (\sqrt{k(k+2)}\cdot\sqrt{1+2^{-p}})^2 \\ &= (k+2)^2\cdot(1-2^{1-p}+2^{-2p}) - k(k+2)\cdot(1+2^{-p}) \\ &= 2^{-2p}(k+2)\cdot\left( (2^{2p}-2^{1+p}+1)(k+2)-(2^{2p}+2^p)k\right) \\ %&= 2^{-2p}(k+2)\cdot\left( (2^{2p}-2^{1+p}+1)k+2(2^{2p}-2^{1+p}+1)-(2^{2p}+2^p)k\right) \\ %&= 2^{-2p}(k+2)\cdot\left( (-2^{1+p}+1)k+2(2^{2p}-2^{1+p}+1)-2^p\cdot k\right) \\ &= 2^{-2p}(k+2)\cdot\left( (-2^{1+p}-2^p+1)k+2(2^{2p}-2^{1+p}+1)\right) \\ \end{aligned} $$ となる。$(-2^{1+p}-2^p+1)k+2(2^{2p}-2^{1+p}+1)$ の部分が、主要項が負の $1$ 次関数であることに注意すると、これが $k=0$ および $k=2^{p-1}$ で正であることを示せば、$0\le k\le 2^{p-1}$ の範囲でこれが正であることがわかる。$p\ge 3$ の前提で、 $$ \begin{aligned} &\phantom{{}={}} (-2^{1+p}-2^p+1)\cdot 0+2(2^{2p}-2^{1+p}+1) \\ &= 2^{1+2p}-2^{2+p}+1; \\ &\phantom{{}={}} (-2^{1+p}-2^p+1)\cdot 2^{p-1}+2(2^{2p}-2^{1+p}+1) \\ &= 2^{1+2p}+2^{p-1} - (2^{2p}+2^{2p-1}+2^{2+p}) + 2 \\ &\ge 2^{1+2p}+2^{p-1} - 2^{1+2p} + 2 \\ &\gt 2^{p-1}+2 \end{aligned} $$ より、これらは正であることがわかる。

すなわち、$0\le k\le 2^{p-1}$ なる整数 $k$ に対して、$k+2\gt \roundsqrt{\roundp{k(k+2)}}$ が成り立つ。$\qed$

よって、Lemmata 5–6 より、$0\le n\le 2^{2(p-1)}$ なる整数 $n$ に対して $$\floor{\sqrt{n}} \le \floor{\roundsqrt{\roundp n}} \le \floor{\sqrt{n}}+1$$ が成り立ちます。すなわち、上記の関数の正当性が示されたことになります。

note: (n as f64).sqrt() as u64 で計算しているのが $\floor{\roundsqrt{\roundp n}}$ です。n as f64 が $\roundp n$ に、x.sqrt() が $\roundsqrt x$ に、r as u64 が $\floor r$ に対応します。

実装例

Rust

fn u64_floor_sqrt(n: u64) -> u64 {
    let tmp = (n as f64).sqrt() as u64;
    let tmp_m1 = tmp.saturating.sub(1);
    if tmp_m1 * (tmp_m1 + 2) < n { tmp } else { tmp_m1 }
}

C++

uint64_t uint64_floor_sqrt(uint64_t n) {
    if (n == 0) return 0; 
    uint64_t tmp_m1 = std::sqrt(n) - 1;
    return tmp_m1 * (tmp_m1 + 2) < n ? tmp_m1 + 1 : tmp_m1;
}

Python

def int_floor_sqrt(n: int) -> int:
    tmp = int(n**0.5)
    return tmp if tmp**2 <= n else tmp - 1

下記のケースでテストしています。

$n$ $\floor{\sqrt n}$
$0$ $0$
$1$ $1$
$(2^{26}+1)^2-1$ $2^{26}$
$(2^{26}+1)^2$ $2^{26}+1$
$2^{64}-1$ $2^{32}-1$

あとがき

つかれました。おなかがすきました。

証明は一旦 iPad で下書きをしているのですが、組版する段階になって符号が間違っていることに気づいたり、指数部の符号も間違っていたり、示したい命題が違っていたりで散々でした。 下書きの証明をそのまま写すのではなく、組版しながらまた別途考えているので、「あれ、こんな式書いた記憶なくない?」というところで気づいたりしています。

さておき、ある程度は浮動小数点数とも仲よくなってきたんじゃないかなという気持ちがあります。 初心者の頃、「平方根sqrt(n) で求めるのは誤差が出ちゃう」と聞いて、返り値の ±30 くらい前後を探すようなコードを書いていたのが懐かしいです。

また、今回の話を踏まえると、素数判定のときにループの上限を i <= sqrt(n) とするのは正当そうということもわかりますね。 約数列挙についてはどうでしょう。あまり考えていませんが、少なくとも u64f64 で全探索した限りでは、$n=i\cdot(i+1)$ のとき $\floor{\roundsqrt{\roundp n}}=i$ となるようだったので、重複して列挙されることはなさそうな気がします。

いかがでしたでしょうか。 今回の話は、ここ数回の中でも競プロ文脈に比較的近そうな気がします。 とはいえ、この記事を見て「なるほど!」となって、rated のコンテストでも f64 ベースの実装で平方根を求める人がいたらなかなか勇気があるなという気もします。 そんなこともないかも。 証明が嘘だったときに責任を取るのは嫌なので、各自で証明し直してほしい気はします。

当初はこの記事は「書く予定リスト」に含まれていなかったため、予告していた「浮動小数点数初心者向け」の記事は次の次になりそうです。ごめんね。

おわり

おわりです。

*1:「有限桁なので、数学的な厳密な値からは誤差が出て当然」という話ではなく、「無限精度で計算した値を正確に丸めた値が欲しいが、それとは異なる値が返ってきうる」という話です。

(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}$ と書くことにします。 型としては doubleIEEE 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$ です。

floor(n / 10.0) と floor(n * 0.1) は等価?

これを double で計算したときに常に等しくなりますか?というクイズです。 ndouble で表せる整数であることにします。

答え

>>> 6755399441055749 / 10
675539944105574.9
>>> 6755399441055749 * 0.1
675539944105575.0

より、前者の floor と後者の floor は異なることになります。 すなわち、等価ではありません。

解説など

真であるなら証明をする必要があるんですが、偽だったので反例を挙げて終了です。 それはそれとして、「これの範囲では成り立つ」という話をしても損はしないので、してみましょう。

浮動小数点数型についてある程度の知識は持っているものとします。そうでない人のための記事は近々(次か次の次くらいに)書く予定です。

$\gdef\fround#1{(\hspace{-.2em}[#1]\hspace{-.2em})}$

まず、0.1 と言っているのは、double(かつ、多くの処理系でデフォルトになっている丸め方向)の文脈においては $$ \begin{aligned} &\phantom{{}={}} {\footnotesize 0.100000000000000005}{\scriptsize 55111512312578270211}{\tiny 81583404541015625} \\ &= \texttt{ccccccccccccd}_{(16)} \times 2^{-55} \\ &= (\tfrac45(2^{52}-1)+1)\cdot 2^{-55} \\ % &= (\tfrac45(2^{-3}-2^{-55})+2^{-55}) \\ &= \tfrac1{10}+\tfrac15\cdot 2^{-55} \\ \end{aligned} $$ と等しい値になります*1

実数 $x$ を浮動小数点数型(ここでは double)で表せる値に丸めたものを $\fround{x}$ と書くことにします*2。 たとえば、$\fround{0.1} = 0.1 + \tfrac15\cdot 2^{-55}$ や $\fround{10} = 10$ です。

また、浮動小数点数 $x$, $y$ に対して $\fround{x\times y}$ および $\fround{x\div y}$ を、それぞれ $x\otimes y$ および $x\oslash y$ と書くことにします。 浮動小数点数型の値同士の演算のみしか考えていないことに注意してください。

さて、今回やりたいのは $\floor{n\oslash 10}$ と $\floor{n\otimes\fround{0.1}}$ の比較です。

note: double仮数部は $p = 53$ bits なので、$2^{53}$ 以上の(浮動小数点数型の)値 $x$ は整数となるため、$\floor{x}=x$ となります。

実数 $x\ge 0$ と整数 $k$ に対して $\floor{\fround{x}} = k$ となる条件は、$\fround{x}$ の整数部分が $k$ となることです。丸め誤差の関係で、$x$ の整数部分が $k$ であることとは同値ではありません。

さて、浮動小数点数で表せる整数 $n = 10q+r$ ($q\in\N, 0\le r\le9$) を考えます。$2^e\le 10q+r \lt 2^{e+1}$ とします。$0\le e\lt p$ としておきます。

$\fround{0.1} = \tfrac1{10} + \tfrac15\cdot 2^{-55} = \tfrac1{10}\cdot(1+2^{-54})$ なので、 $$ \begin{aligned} (10q+r)\cdot\fround{0.1} &= (q+\tfrac r{10})\cdot(1+2^{-54}) \\ &= q+(\tfrac r{10} + (q+\tfrac r{10})\cdot 2^{-54}) \end{aligned} $$ となります。一旦 $r$ が最大の $9$ である場合について考えます。 $$ \begin{aligned} (10q+r)\cdot\fround{0.1} &\le q+(\tfrac9{10} + \tfrac1{10}(10q+r)\cdot 2^{-54}) \\ &\lt q+\tfrac9{10} + \tfrac1{10}\cdot 2^{e+1-54} \\ &= q+\tfrac9{10} + \tfrac1{10}\cdot 2^{e-53}. \end{aligned} $$ $\tfrac1{10}(10q+r)\lt\tfrac18(10q+r)\lt 2^{e+1-3}$ より、$\tfrac1{10}(10q+r)$ の 0.5 ULP は $2^{e-56}$ 以下です。 小数部分を考えて、$\tfrac9{10}+\tfrac1{10}\cdot 2^{e-53} \lt 1-2^{e-56}$ の範囲においては、$(10q+r)\cdot\fround{0.1}$ は $q$ の方に丸められます。 これにより、$e\le51$ の範囲では $\floor{n\oslash 10} = \floor{n\otimes\fround{0.1}}$ となります。 また、同様の考察により、$r\le 8$ では $e=52$ でも $\floor{n\oslash 10} = \floor{n\otimes\fround{0.1}}$ となることが示せます。

よって、$\floor{n\oslash 10} \ne \floor{n\otimes\fround{0.1}}$ となる最小の $n$ は、$e = 52$ かつ $r = 9$ においてこれを満たす最小の $q$ によって作られることがわかります(この時点では存在性は示していませんが、これから作るので大丈夫です)。 $(10q+r)\cdot\fround{0.1}$ の小数部分を考えて*3、 $$ \begin{aligned} \tfrac9{10} + (q+\tfrac9{10})\cdot 2^{-54} &\ge 1-2^{-4} = \tfrac{15}{16} \\ q+\tfrac9{10} &\ge \tfrac3{80} \cdot 2^{54} \\ &= \tfrac35\cdot 2^{50}. \end{aligned} $$ よって、条件を満たす最小の $q$ は $\ceil{\tfrac35\cdot 2^{50}-\tfrac9{10}} = \tfrac15\cdot(3\cdot2^{50}-2)$ となります。 これにより、最小の $n$ は $10\cdot\tfrac15\cdot(3\cdot 2^{50}-2)+9 = 3\cdot2^{51}+5$ となります。これは、冒頭に挙げた例そのものです。

disclaimer: ここでは、該当の範囲で $(10q+r)\oslash 10 = q$ となることをちゃんと示していません。よしなに。

一般のケース

一般に、$\floor{n\oslash d}$ と $\floor{n\otimes\fround{\tfrac1d}}$ を考えてみます。 $d = (2k+1)\cdot 2^s$ とします。二冪のケースは指数部が変わるだけなので $k\ge 1$ とします。 $$ \tfrac1d = \tfrac{2^t}{2k+1}\times 2^{-(s+t)} $$ と表せます。ここで、$t$ は $2^{p-1}\le\tfrac{2^t}{2k+1}\lt 2^p$ となるように定めます。 すなわち、$t = (p-1)+\ceil{\log_2(2k+1)}$ となります。 $$ \tfrac1d = \underbrace{\tfrac{2^{(p-1)+\ceil{\log_2(2k+1)}}}{2k+1}}_{\mu}\times 2^{-(s+t)}. $$ この $\mu$ の部分を $\floor{\mu}$ または $\ceil{\mu}$ に適切に丸めたものを $m$ として、 $$ \fround{\tfrac1d} = m\times 2^{-(s+t)} $$ となります。 $$ (2^{(p-1)+\ceil{\log_2(2k+1)}})\bmod (2k+1) \le k $$ であれば $m=\floor{\mu}$ に、そうでなければ $m=\ceil{\mu}$ になります。 よって、$0\lt u\le k$ とおき、丸め方向に応じて適切な符号を用いることで $$ \fround{\tfrac1d} = \tfrac{2^t\pm u}{2k+1}\times 2^{-(s+t)} $$ と書くことができます。$2^t$ の持つ素因数が $2$ のみであることと、$2k+1\ge3$ が $2$ 以外の素因数のみを持つことから、$u=0$ にはなりません。

符号ごとに分けて考えます。

まず、$\fround{\tfrac1d} = \tfrac{2^t+u}{2k+1}\times2^{-(s+t)}$ のケースです。$\fround{\tfrac1d}\gt\tfrac1d$ に注意します。 $n = dq + r$ とし、$2^e\le d q+r\lt 2^{e+1}$ とします。

$$ \begin{aligned} (dq+r)\cdot\fround{\tfrac1d} &= (dq+r) \cdot \left(\tfrac{2^t+u}{2k+1}\times 2^{-(s+t)}\right) \\ &= (dq+r) \cdot \tfrac{2^t+u}d\cdot 2^{-t} \\ &= (q+\tfrac rd) \cdot (1+u\cdot 2^{-t}) \\ &= q + \tfrac rd + (q+\tfrac rd)\cdot u\cdot 2^{-t} \\ &= q + \tfrac rd + (q+\tfrac rd)\cdot u\cdot 2^{-(p-1)-\ceil{\log_2(2k+1)}}. \end{aligned} $$ あとは、先ほど同様に小数部分を考えていけばよいでしょう。疲れてきたので割愛です。

次に、$\fround{\tfrac1d} = \tfrac{2^t-u}{2k+1}\times2^{-(s+t)}$ のケースです。$\fround{\tfrac1d}\lt\tfrac1d$ に注意します。

いきなりですが、実はめちゃくちゃ反例があります。下記は Haskell の対話型ツール GHCi での出力です。

ghci> take 10 . map round $ filter (\x -> x * (1.0 / x) /= x / x) [1..]
[49,98,103,107,161,187,196,197,206,214]

$49\otimes\fround{\tfrac1{49}} \lt 1$ や $197\otimes\fround{\tfrac1{197}} \lt 1$ などが成り立ちます。先の $\fround{\tfrac1{10}}$ のケースでは身近な範囲(たとえば int に収まる程度)では大丈夫だったりしましたが、こちらに関しては身近な範囲でも全然事故りそうですね。 注意として、浮動小数点数の精度によってこの反例は異なります。下記は Rust の対話型ツール evcxr での出力です*4

>> (1..).filter(|&x| 1.0 / x as f32 * x as f32 != 1.0).take(10).collect::<Vec<_>>()
[41, 47, 55, 61, 82, 83, 94, 97, 107, 109]
>> (1..).filter(|&x| 1.0 / x as f64 * x as f64 != 1.0).take(10).collect::<Vec<_>>()
[49, 98, 103, 107, 161, 187, 196, 197, 206, 214]

この $49\otimes\fround{\tfrac1{49}}$ のようなものについては、書きたいことがいろいろあるので、別の記事で書くことにします。

また、興味のある読者は floor で満足せず、ceil について考えてみてもよいでしょう。ところで $\ceil{x} = -\floor{-x}$ や $\fround{-x} = \fround{x}$ などが成り立ちます。 f32 においては、仮数部が 24 bits しかないため、それなりに十分なテストを行うこともできます。C++ では floating-point promotion によって float が勝手に double になったりするため、Rust などを用いる方が楽そうです。

速度に関して

整数の除算が遅いというのは比較的有名な話かと思います。

そこで、$\sum_{i=1}^{10^{10}} \floor{\tfrac i{10}}$ を、整数型 (i64) と浮動小数点数型 (f64) の $i\oslash 10$ 版と $i\otimes \fround{\tfrac1{10}}$ 版で計測してみました。 また、定数除算最適化の影響を考慮するため、除数を定数として用意するものと実行時に与えるものを比較しました。

もちろん値自体は $O(1)$ で(というか紙とペンで)求めることができて、 $$ \begin{aligned} \sum_{i=1}^{10^{10}} \floor{\tfrac i{10}} &= 9\sum_{i=0}^{10^9-1}i + \sum_{i=1}^{10^9}i \\ &= 10\sum_{i=1}^{10^9-1}i + 10^9 \\ &= 5\cdot (10^9-1)\cdot 10^9 + 10^9 \\ &= 5\cdot 10^{18} - 4\cdot10^9 \\ &= 4999999996000000000 \end{aligned} $$ です。これは $2^{63}$ よりも小さいことを指摘しておきます。

足し合わせる変数には i64 を用い、除算パートのみで上記の違いがある状態にします。 なお、$10^{10} \lt 3\cdot 2^{51}+5$ なので、この範囲で正しい値を計算できることに注意してください。 計測結果は下記の通りでした。

i64 f64 $\oslash\, {10}$ f64 $\otimes\,\fround{0.1}$
定数 2.26 s 2.86 s 2.24 s
変数 5.75 s 2.90 s 2.20 s

どうやら、いずれも f64 で $\otimes$ を用いるものが最も速かったようです。あらあら。 あと、やはり変数での除算は重いようですね。

今回は floor の方のみ実験しましたが、ceil の方においても適切に丸めた方を用いたりすることで高速に正しい値を得ることが可能そうな気がします。 浮動小数点数型と仲よくなると、適切に証明した上で高速化をすることができそうでうれしいですね。

おきもち

コンテスタント目線では、この手の誤差評価などは “非本質” などと扱われがちな部分な気はしています。 一方で、先ほど言及したように高速化に寄与する場合があり、知っていると役に立つ場面もあるのかなとも思います。

あるいは、テストケースを作成する側に回った際にも、「浮動小数点数型を用いたこうした解法は実は正当である」とか、「誤差がこうなることがわかるので、こういうケースを置いておけば落とせそう」とか、そうした考え方ができるとうれしそうです。AtCoder では(特に double を想定したケースは)考慮されていがちな気はしますよね(long double は見落とされがちな気もしますが)。

ソフトウェア開発の文脈だと、「数学的な証明をしました」で単体テストを書かないのは不安な感じがある気がします。 一方、この手の浮動小数点数が絡むものにおいては、単体テストで適当にいくつかのケースを選んで「大丈夫そうですね」と言うのは不安な気持ちになります。テストしていないケースにおいて誤差でめちゃくちゃになることが容易に起こりうるからです。多少詳しい人でないと、先の例の $\fround{\tfrac1{49}}$ などの反例を挙げることは難しそうですよね。

とはいえ、浮動小数点数に明るくない人(世の中のマジョリティ)にそういう事情を説明してわかってもらうのは大変そうな気がしていて、うまいことやっていく必要がありそうだなという気持ちになったりしています。 この辺の話は正規表現などのトピックでも当てはまりそうですね。

特になにかあったわけではないのですが、競プロ寄りの人々には興味のない話かなと思って薄めの文字で書きました。

あとがき

$\fround{x}$ の記法は、今回思いつきで導入したものになります。 下記の前回の記事では $\hat x$ とか $\operatorname{round}(x)$ とか書いていたものですが、長めの式に対してもよい感じに書けるものがいいなと思ったので考えてみました。

rsk0315.hatenablog.com

経緯としては、最近(文脈自体はあまり関係ないものの)$[\![\texttt{a+b}]\!]$ のような括弧が 使われているの を見て、「括弧をそうやって作るのも “アリ” なんだな」と思ったのがきっかけになります。操作としては round を表すので丸い括弧がいいなと思って、$(\hspace{-.2em}|x|\hspace{-.2em})$ や $\fround{x}$ などを試して決めたものです。 $(\hspace{-.2em}|x|\hspace{-.2em})$ もいいなと思ったのですが、別の文脈*5での使われを思い出したことと、$[\ldots]$ の方が Gauß 記号や $\floor{\ldots}$ $\ceil{\ldots}$ を連想させて、丸めに近いかな?という気持ちがあったことによります。定義は下記です*6

$\gdef\fround#1{(\hspace{-.2em}[#1]\hspace{-.2em})}$

もしかしたらこの括弧に別の意味があったり、該当の概念に別の記法があったりするのかもしれませんが、あまりサーベイしていません。

おわり

最近は需要のなさそうな記事ばかり書いていますが、何らかの局面で誰かしらの役に立ったらいいなと思っています。あるいは他の人の役に立っていなくても、書いているときの自分が気持ちよくなっているので十分です。

*1:「誤差でぐちゃぐちゃになった意味不明な値」のようなイメージがある人もいるかもしれませんが、数学的に書ける値であることを意識するとよいかもです。

*2:仕様に応じて適切な丸め方向を採用するものとします。

*3:仮数部の偶奇に応じて等号を含むかが変わりますが、一旦等号を含むものから考えてみて、だめならその次のものを採用することにします。

*4:いろいろな対話型ツールがあってたのしいですね。

*5:The SATySFibook の record の部分です。こっちも別の出典があるのかも。分野がやや違いそうなので、同じ記号でも気にしてなくていいと思ってはいます。

*6:論文とかで何らかの独自の記法を使うときは、$\LaTeX$ での書き方を footnote にでも書いておいてくれと思っているタイプです。