「decimal は演算が正確である」という主張をちゃんと書く

## 概要

よくある

* (二進) 浮動小数点数は正確な演算結果を返さないことがある

という話題に

* Python の `decimal` (等の十進浮動小数点数) なら正確な演算が可能

といったコメントが付いたりしてそれに対して

* `decimal` でも正確な演算ができないことはある

という指摘があったりしてじゃあ何故「`decimal` は演算が正確である」という印象を持たれるのか考えるために問題をきちんと定式化して論じようと思った.  

この記事は「主張をちゃんと書こう」という試みに過ぎず, その主張が妥当かどうかはまた別の問題であることに注意されたい.  

## 準備

まずは「リテラル」「データ」「演算」を定義して「正確な演算結果を返さない」という言葉が厳密に意味するところを明らかにしておく.  

### リテラル

(実数の) **リテラル**は十進法による有限文字列として表現できる実数であり, $10^dn$ ($n, d \in \mathbb{Z}$) の形の有限小数として定式化される[^literal].  
この集合を,  

[^literal]: 実際はリテラルは表記方法に依存し, `0.01` のような有限小数表記や `1e-2` のようないわゆる E 表記はリテラルとなる (ことが多い) が, ここでの定義に従う `1.0*10**(-2)` のような指数表記や `1/100` のような有理数表記はリテラルではなくリテラルと演算子を組み合わせた式であるとされる (ことが多い).  

$$ 
\mathcal{L} := \\{10^dn \ \vert \ n, d \in \mathbb{Z}\\}, 
$$

と書く.  
$n' := 10n, d' := d-1$ に対して $10^{d'}n' = 10^dn$ となる等, リテラルのパラメータ $n, d$ は一意ではないことに注意しておこう.  

### データ

**データ**はコンピュータ上で実際に扱われる数であり, ここでは整数でパラメータ付けられた一般形を持つ数と定義しておく[^data].  
例えば `decimal` による十進浮動小数点数や固定小数点数はリテラルと同じく $10^dn$ と記述され, 二進浮動小数点数は $2^dn$ という一般形を持つためいずれもデータであると言える[^float].  

[^data]: ただしどのような関数を使ってパラメータ付けられるものかは定めていないので, 可算無限集合がデータであるということになる.  

[^float]: もっと一般的には, 浮動小数点数はこれらを正整数の定数 $r$ で割った数に拡張される. しかし標準的には $p$ 進浮動小数点数において $r$ は $p$ の累乗で定められるので今回の話の枠組みでは無視できる.  

通常はリテラルのパラメータは任意の値を取り得る[^parameter]のに対して, データのパラメータは有限範囲に制限されるためデータの集合は有限集合となる[^finiteset].  
議論の都合によりパラメータの制限を考えないデータを扱いたい場合もあり, その時データは**任意精度である**ということにしよう.  

[^parameter]: メモリや入力する人間の体力等の環境的な制約を除く. この辺の差別化は面倒なのでリテラルのパラメータは任意の値を取れるとしておく.  

[^finiteset]: [pure-Python の `decimal` は $n$ として任意の整数を取れる](https://stackoverflow.com/questions/28081091/what-is-the-largest-number-the-decimal-class-can-handle)ようなので有限集合にならないが, 本質的ではないところで厳密に語るのは煩雑になるのでデータ集合は有限集合であるとしておく.  

データは実数をコンピュータ上で扱うため制限を加えたものであり, **丸め関数**と呼ばれる実数集合 $\mathbb{R}$ からデータ集合 $\mathcal{X}$ への写像 $\rho: \mathbb{R} \rightarrow \mathcal{X}$ が用いられる.  
これにより実数 $x \in \mathbb{R}$ を扱いたい時は $\rho(x) \in \mathcal{X}$ を代わりに扱い, この $\rho(x)$ を得る操作を $x$ を**丸める**という.  
以降は実数を丸める操作が頻出するため, $x_\rho := \rho(x)$ という記法を導入しておこう.  

丸め関数に求められる性質としてはとりあえず,  

* データ集合上で恒等: $x \in \mathcal{X} \Rightarrow x_\rho = x$, 
* 広義単調増加: $x \leq y \Rightarrow x_\rho \leq y_\rho$, 

の二つを仮定しておく.  
これらの性質を持つ丸め関数は任意精度のデータに対しては一般的に定めることができず, 実際にプログラミング言語において任意精度のデータを扱う場合は, 任意のパラメータの制限付きデータ集合とそれに適した丸め関数を都度選択することで任意精度を実現する.  

### 演算

**演算**は写像とほぼ同義であり, 特に一つ以上のいくつかの集合の直積を定義域とするものを考える.  
$n$ 個の集合の直積を定義域とする演算を $n$ 項演算 ($n = 1$ なら単項演算) 等と呼ぶこともあり, 以降は簡単のため, 特に言及しない場合は二つの集合の直積からの写像である**二項演算**を対象にして話を進める.  
ここではリテラルやデータを入出力とするものを考えたいので, ($n$ 項) 演算は定義域が $\mathbb{R}^n$ である関数 (つまり値域は $\mathbb{R}$) としておこう.  

データ集合 $\mathcal{X} \subset \mathbb{R}$ 及び丸め関数 $\rho: \mathbb{R} \rightarrow \mathcal{X}$ を用いて二項演算 $f: \mathbb{R}^2 \rightarrow \mathbb{R}$ による演算結果を得たい場合,  

$$ 
f_\rho(x, y) := f(x_\rho, y_\rho), 
$$

で定まる演算 $f_\rho$ の結果をさらに $\rho$ で丸めた $(f_\rho(x, y))_\rho$ を $f(x, y)$ の代用とするのが通例になる[^surrogate].  
この $(f _\rho(x, y)) _\rho$ を返す関数の記号を定めないのは, 後に丸める前の $f _\rho(x, y)$ 自体を評価する必要が出て来ることの他, 後述の $= _\rho$ や演算の合成において丸めた値を再度丸める手間を省くためである.  

[^surrogate]: 計算効率の理由で厳密にこの手続きを踏まない演算も考えられるがここでは無視する.  

### 正確な演算結果を返す

データ集合 $\mathcal{X}$ と丸め関数 $\rho: \mathbb{R} \rightarrow \mathcal{X}$ を選んだ時,  

$$ 
\alpha =_\rho \beta \ \overset{\text{def}}{\iff} \ \alpha _\rho = \beta _\rho, 
$$

という記号 $=_\rho$ を用いて, リテラル $x, y \in \mathcal{L}$ で $f(x, y) \in \mathcal{L}$ となる物に対して,  

$$ 
f_\rho(x, y) =_\rho f(x, y), 
$$

が成り立つ, というのが「$f(x, y)$ が ($\mathcal{X}, \rho$ に関して) **正確な演算結果を返す**」という言葉の定義である.  
これは単項演算及び三項以上の演算に対しても同様に定義される.  
前提となる $x, y \in \mathcal{L}$ が $f(x, y) \in \mathcal{L}$ となるという条件は, **正確な演算結果を返すか判定可能である**と言うことにしよう.  
勝手な定義に思えるかもしれないが, この定義は「リテラル上で $x+y = z$ が成り立つ時に `x+y==z` が `True` を返すか」といった議論から構成したもので妥当な定義だと言えるだろう (と思ってもらいたい).  

注意として, ここで定義した「$f(x, y)$ が正確な演算結果を返さない」ことと「$f(x, y)$ が丸めによって直感と異なる演算結果を返す」こととは別の現象である.  
例えば Python の `float` による以下の例では $0.1+_\rho10^{-20} = _\rho 0.1$ となっている (対応するコードが `True` を返す) ため $0.1+10^{-20}$ の `float` における結果は直感に反している (はず) だが, $0.1+ _\rho10^{-20} = _\rho 0.1+10^{-20}$ であるため, $0.1+10^{-20}$ は `float` 及びその丸めに関して正確な演算結果を返していると言える.  
ただし, $+ _\rho$ は二項演算 $\mathrm{add}: (x, y) \mapsto x+y$ に対して $x+ _\rho y := \mathrm{add} _\rho(x, y)$ を返す演算子である.  

```python
print(0.1+0.00000000000000000001==0.1)
print(0.1+0.00000000000000000001==0.10000000000000000001)
```

    True
    True

### 演算の合成

本題に入る前に演算が複数の演算の合成として実装されている場合について言及しておく.  

対象の演算 $f$ がプログラム上で複数の演算 $f_0, f_1, \dots, f_{r-1}$ の合成写像として実装されている時, 個々の演算に丸めを加えた $(f_0)_\rho, (f_1) _\rho, \dots, (f _{r-1}) _\rho$ の合成を考える必要がある.  
$f$ が複数の演算の合成で実装されていない時 (演算の合成への分解を考えない時), $f$ は**単純な演算である**と言うことにしよう.  

単純でない $f$ が $x, y \in \mathcal{L}$ について正確な演算結果を返すか判定可能であることは, $f$ を構成する各演算 $f_i$ が正確な演算結果を返すか判定可能であることを意味しない.  
[参考資料二つ目](https://qiita.com/mod_poppo/items/dd2ac0af80993751ee5a)での,  

$$ 
\begin{gather*} 
f(x, y) := \mathrm{mul}(\mathrm{div}(x, y), y), \\\\ 
\mathrm{div}(x, y) := \frac{x}{y}, \quad \mathrm{mul}(x, y) := xy, \\\\ 
\end{gather*} 
$$

がその例で, $f(x, y) = x$ と整理できるため $f$ 自体は ($y = 0$ を除いて) 常に正確な演算結果を返すか判定可能だが, $\mathrm{div}$ については $x, y \in \mathcal{L}$ で $\frac{x}{y} \notin \mathcal{L}$ となる組が存在するためその時は正確な演算結果を返すか判定可能ではない.  

## 本題

### `decimal` は演算が正確だと言われる要因

まず簡単な事実として, リテラル $x, y \in \mathcal{L}$ について単純な演算 $f$ が正確な演算結果を返すか判定可能である時, $x, y$ を共に含むデータ集合 $\mathcal{X}$ と丸め関数 $\rho: \mathbb{R} \rightarrow \mathcal{X}$ に関して $f$ は $f(x, y) \notin \mathcal{X}$ であったとしても正確な演算結果を返す.  
何故なら $\rho$ は $x, y \in \mathcal{X}$ において恒等であり,  

$$ 
f_\rho(x, y) = f(x_\rho, y_\rho) = f(x, y), 
$$

より, 両辺に $\rho$ を作用させても等式は保たれて $f_\rho(x, y) =_\rho f(x, y)$ が成り立つためである.  

ここで, 十進浮動小数点数や固定小数点数は $10^dn$ という**リテラルと同じ一般形を持つ**ということが肝になることが分かる.  
つまり適当に選んだリテラルに対して, パラメータの範囲を十分に広く取れば選んだリテラルはデータ集合に属し, その時, 正確な演算結果を返すか判定可能な単純な演算は常に正確な演算結果を返すことになる.  

例えば具体的な演算といくつかのリテラルを選んで実際に正確な演算結果を返すか検証するとしよう.  

* 加算: $\mathrm{add}(x, y) := x+y$, 
* 減算: $\mathrm{sub}(x, y) := x-y$, 
* 乗算: $\mathrm{mul}(x, y) := xy$, 

といった基本的な演算を検証に用いる単純な演算として選んだ時, これらは任意のリテラルに対して正確な演算結果を返すか判定可能であり, さらに検証に用いるリテラルとして短い文字列で記述されて一般の十進浮動小数点数に収まるものを取ると, 必ず「十進浮動小数点数についての演算は常に正確な演算結果を返す」という結論が得られてしまう[^concluded].  
これが「`decimal` は演算が正確である」という主張が実際に意味するところだと考えられる.  

[^concluded]: 一応言及しておくと, 「十進浮動小数点数でも常に正確な演算結果を返すとは限らない」という主張と「十進浮動小数点数では特定の (しかし一般的な) 条件下で正確な演算結果を返すことが保証される」という主張は両立する. 要はデータや演算の性質を正確に把握して適切に運用することである.  

### 演算の合成による反例

先に触れた[参考資料二つ目](https://qiita.com/mod_poppo/items/dd2ac0af80993751ee5a)はこれを踏まえた反例で,  

* 除算: $\mathrm{div}(x, y) := \frac{x}{y}$, 

を用いた $f(x, y) = \mathrm{mul}(\mathrm{div}(x, y), y)$ という演算で検証を行っている.  
この場合, $\mathrm{div}$ は常に正確な演算結果を返すか判定可能とは限らず, $\mathrm{div}(x, y) \notin \mathcal{L}$ の時に続く $\mathrm{mul}$ にデータ集合に属さない値が渡されることになる.  
これにより前節の議論を適用できず, $f_\rho(x, y) =_\rho f(x, y)$ が成り立たない例が得られる.  
演算の合成によらず $f(x, y) = \mathrm{div}(x, y)$ とした場合は, $\mathrm{div}(x, y) \notin \mathcal{L}$ ならばリテラルを用いて $f _\rho(x, y) = _\rho f(x, y)$ が成り立つか確認することができず (つまり正確な演算結果を返すか判定可能でない), $\mathrm{div}(x, y) \in \mathcal{L}$ の時は $x, y \in \mathcal{L}$ が共にデータであれば正確な演算結果が返る.  

なお, じゃあ常に正確な演算結果を返すか判定可能な演算の合成なら大丈夫かというとそういうわけではなく, 全引数をデータ集合から取ったとしても前段の演算結果がデータ集合から外れていれば反例を得るのに十分である.  
次は $f(x, y, z) = \mathrm{add}(\mathrm{add}(x, y), z)$ として $\mathrm{add}(x, y)$ がデータ集合から外れるように $x, y, z$ を与えた時の例である.  
この例では `getcontext().prec = 1` により $-9 \leq n \leq 9$ に制限されている.  
$d$ についてはこの書き方だと $\pm999999999999999999$ が上下限とのこと[^range].  

[^range]: https://stackoverflow.com/questions/28081091/what-is-the-largest-number-the-decimal-class-can-handle, 脚注[^finiteset]も参照.  

```python
from decimal import *

getcontext().prec = 1
x = Decimal(6)
y = Decimal(6)
z = Decimal(3)
# インスタンス化時には getcontext().prec の精度で丸められない仕様なので 0 を足しておく
w = Decimal(6+6+3)+Decimal(0)

print(x+y+z==w)
```

    False

### 二進浮動小数点数との比較

前節と同様の議論で, 二進浮動小数点数に関して正確な演算結果を返さない場合について考えてみる.  

二進浮動小数点数の場合, 任意精度であっても表現できないリテラルが存在する.  
特に, 任意の $\alpha, \beta \in \mathbb{R}$ ($\alpha \lt \beta$) に対して $\alpha \lt x \lt \beta$ であるリテラル $x \in \mathcal{L}$ で任意精度の二進浮動小数点数で表現できないものが存在し (演習), どのように範囲を決めても正確な演算結果を返さない例が存在し得ることを示唆している[^suggest].  

[^suggest]: 「任意の二項演算 $f$ と二進浮動小数点数で表現できないリテラル $x \in \mathcal{L}$ に対して $y \in \mathcal{L}$ が存在して $f$ は $x, y$ について正確な演算結果を返さない」みたいなことを言えれば綺麗にまとまったが, $y$ を無視して $x$ を返すような演算を考えれば常に正確な演算結果を返すのでそれも無理だった.  

このように,  

* 十進浮動小数点数では正確な演算結果が返る条件を比較的容易に構成・実現できる,  
* 二進浮動小数点数では正確な演算結果が返るケースが強く限定される,  

という差から「`decimal` は二進浮動小数点数より演算が正確である」という認識が発生していると考えられる.  
また, $2$ は $10$ の約数であり, $d \lt 0$ の時 $2^dn = 10^d(5^{-d}n)$ と変形できて, 任意精度二進浮動小数点数は常に任意精度十進浮動小数点数によって記述できる.  
つまりパラメータの範囲を十分大きく取ると二進浮動小数点数に関して正確な演算結果を返すならば十進浮動小数点数に関しても正確な演算結果を返し, このこともこの認識の後押しとなっているだろう.  

リテラルとして二進数を文字列表記したもの, つまり $2^dn$ で定まる数を取ってもこれは常に $10^{d'}n'$ の形式でも書けるため, 上の意味の優劣関係は同等にはなっても逆転はしない (有効桁数の確認が面倒になるぐらい).  
これを破るには任意の $d$ に対して $10^d$ の約数でない整数 $p$ による $p$ 進数リテラルが必要だが, そうすると今度は「$p$ は $2^d$ の約数でない」も成り立つので, やはり二進浮動小数点数の方が有利とはならない.  

### その他の議論

$p$ 進浮動小数点数以外のデータに $p/q$ でパラメータ付けられる有理数がある.  
有理数は「任意進浮動小数点数」とも言えるデータで, 十進浮動小数点数の場合と同様にリテラルは常に有理数の集合に含まれるため $\mathrm{add}, \mathrm{sub}, \mathrm{mul}$ はパラメータの範囲が十分大きい有理数に関して常に正確な演算結果を返す.  
さらに $\mathrm{div}$ に関しても, 入力が共にリテラルならば $0$ 除算でない返り値はリテラルでないとしても常に有理数となっており, 先の $f(x, y) = \mathrm{mul}(\mathrm{div}(x, y), y)$ の例は $\mathrm{div}$ が正確な演算結果を返すか判断可能かどうかに関わらず, $0$ 除算が発生する場合を除いて $f$ は正確な演算結果を返すことが示される.  
当然, パラメータの制限から外れれば $f(x, y, z) = \mathrm{add}(\mathrm{add}(x, y), z)$ の例のようなことが発生する.  
これは有理数以外のより広いデータを採用した場合も同様の現象が起こる.  

演算について三角関数等のより一般的な関数を考えることは, ここでの「正確な演算結果を返すかどうか」という問題においてはあまり有用ではない.  
まず演算を厳密に計算できるアルゴリズムが存在するかという別の問題があり, 近似計算する実装方法を取ったとしても, 正確な演算結果を返すか判断可能なリテラルを探してここまでと同じ議論を適用することになるためである.  

## 参考資料

Qiita での関連記事達の一部を挙げる.  

#### 「0.1+0.2≠0.3」を説明できないエンジニアがいるらしい #Python - Qiita

二進浮動小数点数において $0.1+0.2$ が正確な演算結果を返さない詳細な解説とコメント欄で言語や型による差異について濃密な議論が行われている.  

https://qiita.com/higashi_nc/items/9a5ea00415a008f06843

#### decimal型（十進小数）に夢を見ている輩が多すぎる #浮動小数点数 - Qiita

「十進浮動小数点数なら $0.1+0.2$ を正確に計算できる」という主張に対して十進浮動小数点数でも正確な演算結果を返さない例を挙げて有理数やより一般的な実数に話を広げて議論している.  

https://qiita.com/mod_poppo/items/dd2ac0af80993751ee5a

#### 浮動小数点演算で数学的に正確な演算結果を期待したとして(いやしちゃダメなんだけど)裏切られる確率は型によって違うのかな #C++ - Qiita

タイトル通り各数値型で正確な演算結果を返す確率を検証している.  
当初の予定ではこの「浮動小数点演算で数学的に正確な演算結果を期待したとして(いやしちゃダメなんだけど)裏切られる確率」を計算して示すつもりだったが考えてるうちにリテラルと decimal が同じパラメータ付けを持ってるからと述べる方が簡単だと気付いたので没にした.  
記事内及びコメントの実験は小数点以下桁数を増やすとたぶん確率 $1/4$ に収束するんじゃないかと思っているが二つの引数がちょうど同じになる部分に引っ張られて少し小さめの確率になっているんじゃないかと予想している.  

https://qiita.com/fujitanozomu/items/b0f93067f49216ccf71a