えびちゃんの日記

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

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

浮動小数点数 ボリューム大きめ

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

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

まえがき

浮動小数点型といえば、競プロ 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]。

精度という語について

精度という言葉は、適切でない使われ方をすることが多そうな気がします。 英語では precision と accuracy という異なる 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.0 も 3.0 + 1e-324 も 2.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 double も Decimal も信用ならないですね。

こうした誤差が起きる理由は、先の $\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-math や cr-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.0 が true だったり、nan == nan が false だったりします。 これに関しては、$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_VALUE で 1.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+01 と 1e+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} $$ が得られていることにも注意が必要です。

たとえば Python で round(1.15, 1) のように書くと、1.1 が得られます。 一方、Ruby で 1.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 1 が 1.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:Unicode と UTF-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)$ よりも寄与が小さい項は無視できる」のようなイメージで捉えるとよいと思います。