Pythonの計算機イプシロン


## 計算機イプシロンとは

計算機イプシロンとは，計算機の浮動小数点における「1より大きい最小の値」と1との差のことです。
有限の桁しか持たない浮動小数点での2つの値の一致判定などに使用されます。

### よくある誤解

計算機イプシロンは，しばしば「浮動小数点で表せる0より大きい最小の数」と誤解されます。
もしこのような定義だと，指数部分を小さくすることで仮数部分の精度に関係なく小さな値を表現できてしまうことになります。
そのため，「1より大きい」という少しややこしそうな定義になっています。
(この記事を書こうと思ったきっかけも，「計算機イプシロン Python」でググった時に出てきた記事で誤った計算方法を使用しているのを発見したからです)

## 計算方法

### 定義通りの実装

計算機イプシロンの定義通りに計算する方法を実装すると，以下のようになります。

```python
EPS = 0.0
tmp = 1.0
while 1.0 + tmp > 1.0:
  EPS = tmp
  tmp /= 2
print(EPS)
```

### ループを使用しない実装

コメントで教えていただいた，ループを使用しない方法です。

```python
a = 4 / 3
b = a - 1
c = b + b + b
EPS = 1 - c
print(EPS)
```

<details><summary>解説</summary><div>

この方法はコードを眺めただけでは何が起きているのか分かりづらいため解説を載せておきます。
(倍精度浮動小数点数の詳細は割愛します)

計算機イプシロンとは，定義から

```
                                                         仮数の一番下が1
                                                         ↓
   0b10000000000000000000000000000000000000000000000000001 × 2⁻⁵² … 1より大きい最小の値
-) 0b10000000000000000000000000000000000000000000000000000 × 2⁻⁵² … 1
------------------------------------------------------------------
   0b00000000000000000000000000000000000000000000000000001 × 2⁻⁵²
```
であるので，上記のコードでこれが計算できていることを確認します。


`4 / 3`の計算結果は`0x3ff5555555555555`であり，

- 符号`0x0`
- 指数部`0x3ff=1023`
- 指数バイアス`0x3ff=1023` … 定義値
- 仮数部`0x5555555555555`

であるので，値は
$$
a = 2^{1023-1023}\times(\texttt{0x15555555555555}\times2^{-52})
= 1.3333333333333333
$$
となります。
(正確には1.3333333333333332593184650249895639717578887939453125)

`a`の値から1を引いてみると

```
   0b10101010101010101010101010101010101010101010101010101 × 2⁻⁵² … a
-) 0b10000000000000000000000000000000000000000000000000000 × 2⁻⁵² … 1
------------------------------------------------------------------
   0b00101010101010101010101010101010101010101010101010101 × 2⁻⁵²
 = 0b10101010101010101010101010101010101010101010101010100 × 2⁻⁵⁴
(= 0x3fd5555555555554)                                   ↑
                                                         精度の都合で切られる
```
となるので，`b = 0.33333333333333326`となります!
(正確には0.3333333333333332593184650249895639717578887939453125)
お気づきの通り，`0x3fd5555555555555`とすると(人間が期待する)0.3333333333333333となりますが，1.3333333333333333-1とすることによって誤差が生じます。
(`1 / 3`の結果は`0x3fd5555555555555`になるので0.3333333333333333です。)

この1/3より微妙に小さい'b'を2つ足してあげると

```
   0b10101010101010101010101010101010101010101010101010100 × 2⁻⁵⁴
+) 0b10101010101010101010101010101010101010101010101010100 × 2⁻⁵⁴
------------------------------------------------------------------
  0b101010101010101010101010101010101010101010101010101000 × 2⁻⁵³
 = 0b10101010101010101010101010101010101010101010101010100 × 2⁻⁵³
(= 0x3fe5555555555554
 = 0.6666666666666665)
```
となり，これにさらに`b`を足すと

```
   0b10101010101010101010101010101010101010101010101010100 × 2⁻⁵³
+) 0b10101010101010101010101010101010101010101010101010100 × 2⁻⁵⁴ ← 1ビットずれている
------------------------------------------------------------------
   0b11111111111111111111111111111111111111111111111111110 × 2⁻⁵³
(= 0x3feffffffffffffe)
```
という一番下のビットだけが0となっている値が得られます。
(計算すると0.9999999999999997779553950749686919152736663818359375)
これより，`c = 0.9999999999999998`となります。

最後に，1から'c'を引いてみると

```
   0b10000000000000000000000000000000000000000000000000000 × 2⁻⁵² … 1
-) 0b11111111111111111111111111111111111111111111111111110 × 2⁻⁵³ … c
------------------------------------------------------------------
   0b00000000000000000000000000000000000000000000000000001 × 2⁻⁵²
(= 0b10000000000000000000000000000000000000000000000000000 × 2⁻¹⁰⁴ … 正規化
 = 0x3cb0000000000000
 = 2.220446049250313080847263336181640625 × 10⁻¹⁶)
```
となり，確かに計算機イプシロンを計算できていることがわかりました。

</div></details>

## 実行結果

Win10のPython3.6.0で実行したところ，「2.220446049250313e-16」が表示されました。
これは，IEEE754のbinary64(倍精度)の計算機イプシロンです。

### ライブラリとの比較

Pythonでは，`sys.float_info.epsilon`なる値が計算機イプシロンを表しています。(Cで言う`FLT_EPSILON`的なやつ)
先ほどと同じ環境で確認してみると，

```py
>>> import sys
>>> sys.float_info.epsilon
2.220446049250313e-16
```

と表示され，自分で計算した値と一致しています。

## 利用例

冒頭で書いたように，計算機イプシロンは浮動小数点の一致判定などで役に立ちます。
この記事をここまで読んでくださった方はご存知の方も多いかと思いますが，コンピュータの浮動小数点では実数を正確に表現することはできません。
例えば，実数では0.1の10倍と0.1を10回足したものは一致しますが，これをPythonで表現すると

```py
>>> a =  0.1*10
>>> b = sum(0.1 for _ in range(10))
>>> a, b
(1.0, 0.9999999999999999)
>>> a == b
False
```

となってしまいます。
コンテキストにもよりますが，一般的にこれは望まれざる挙動でしょう。

そこで，計算機イプシロンを使用して一致判定を行います。
具体的には，相対誤差が十分に小さいとき，2つの値が等しいとみなします。
これをPythonで表現すると

```py
>>> import sys
>>> a =  0.1*10
>>> b = sum(0.1 for _ in range(10))
>>> a, b
(1.0, 0.9999999999999999)
>>> abs(a-b) < sys.float_info.epsilon * max(abs(a), abs(b))
True
```

となります。
相対誤差を使用する理由などの厳密な議論に興味がある方は，参考文献に挙げたページをご覧ください。

## 追加・変更点など

### 追加 (2019/08/30)

計算機イプシロンによる浮動小数点の一致判定の実例を追加しました。

### 追加 (2020/02/21)

コメントで教えていただいた，ループを使用しない計算方法とその解説を追加しました。

## 参考文献

- [桁の落ちない話](https://qiita.com/azapen6/items/fd03c9a8a29a00bb4626)
- [浮動小数点の比較①](http://berobemin2.hatenablog.com/entry/2016/02/27/231856)
- [Pythonで浮動小数点数floatと16進数表現の文字列を相互に変換](https://note.nkmk.me/python-float-hex/)