19
6

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

Pythonの計算機イプシロン

Last updated at Posted at 2018-11-25

計算機イプシロンとは

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

よくある誤解

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

計算方法

定義通りの実装

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

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

ループを使用しない実装

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

a = 4 / 3
b = a - 1
c = b + b + b
EPS = 1 - c
print(EPS)
解説

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

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

                                                         仮数の一番下が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⁻¹⁶)

となり,確かに計算機イプシロンを計算できていることがわかりました。

実行結果

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

ライブラリとの比較

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

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

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

利用例

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

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

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

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

>>> 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)

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

参考文献

19
6
2

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
19
6

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?