3.7 数値型の誤差と範囲

(1/1)
円周率
3.3 比較演算子とブール演算子」で触れたように、コンピュータは2進数で計算するために、循環小数になる場合がある。これを回避する方法はないのだろうか。また、Pythonの計算は、電卓のように10進数8桁や12桁といった限界はないのだろうか。

今回は、数値型の範囲と誤差について学ぶ。最後に、巨大素数や円周率を求めるプログラムを紹介する。

目次

サンプル・プログラム

10進数の循環小数

まず、\( 1 \div 3 \) の計算結果を見てみよう。
# 10進数の循環小数
a = 1 / 3
b = a * 3
print("a = " + str(a))
予想通り、\( 0.3333... \) となる。これは10進数の循環小数だからどうしようもない。だが、この商に \( 3 \) を掛けると、 \( 1 \) に戻る。さすがに安物の電卓とは違う。
では、\( 1 \div 3 \) の商と、\( 7 \div 21 \) の商を比較したらどうなるだろうか。
# 小数の比較
a = 1 / 3
b = 7 / 21
print(a == b)
当然等しくなる(True)。

2進数の循環小数

# 小数の比較
a = 0.3
b = 0.1 + 0.2
print(a == b)
プログラム "float3.py" は、\( 0.3 \) と \( 0.3 \) を比較するだけなのだが‥‥等しくならない(False)。なぜか――。

コンピュータは、数値を内部的に2進数で処理していると説明した。整数は2のべき乗 \( \displaystyle 2^0,\ 2^1,\ 2^2... \) を使って表すことは「3.4 ビット演算子、シフト演算子」で説明したとおりだ。そして、小数点以下も2のべき乗 \( \displaystyle 2^{-1},\ 2^{-2},\ 2^{-2}... \) を使って表す。
わかりやすいように、2のべき乗(マイナス)を分数にしてみると、\( \displaystyle \frac{1}{2},\ \frac{1}{4},\ \frac{1}{8}... \) となる。

\( 0.3 \) を2進数の分数で表そうとすると、\[ \displaystyle 0.3 = \frac{0}{2} + \frac{1}{4} + \frac{0}{8} + \frac{0}{16} + \frac{1}{32} + \frac{1}{64} + \frac{0}{128} +... \] と、無限に足し算を続けることになってしまう。つまり、\( 0.3 \) は2進数の循環小数である。
だが、循環小数になるといっても、前述の \( 1 \div 3 \) と同じで、同じ \( 0.3 \) と \( 0.3 \) を比較するのだから、結果は True になるはずではないだろうか。そこで、プログラム "float4.py" を実行してほしい。
# 小数の比較
a = 0.3
b = 0.1 + 0.2
print(f"{a:.54f}")
print(f"{b:.54f}")
print文の引数は f"{a:.54f}" という見慣れない形式になっているが、これは、フォーマット文字列(「3.8 文字列」で詳しく説明する)と呼ばれ、変数 a に格納されている値を10進数で小数点以下54桁――.54f;fはFloatのf――で表示する、という意味になる。
実行環境によって若干異なるかもしれないが、\[
\begin{eqnarray*}
0.3 &=& 0.299999999999999988897769753748434595763683319091796875 \\
0.1+0.2 &=& 0.300000000000000044408920985006261616945266723632812500
\end{eqnarray*}
\] のように、内部の値が異なっていることがわかる。だから、比較演算で等しくならない。

では、では、\( 1 \div 3 \) の商と、\( 7 \div 21 \) の商はどうだろうか。プログラム "float5.py" を実行してほしい。
# 小数の比較
a = 1 / 3
b = 7 / 21
print(f"{a:.54f}")
print(f"{b:.54f}")
両者とも全く同じ値になる。

2進数の循環小数

それにしても、小数第一位で誤差が発生してしまうとなると、計算機として使い物にならない――どうしたものか。

3.3 比較演算子とブール演算子」で紹介したが、Python には、許容誤差範囲内で引数が等しいかどうかを比較する math.isclose関数 が用意されている。比較演算があるなら、math.isclose関数を使うのが1つの解決策である。デフォルトの許容誤差は \( 1 \times 10^{-9} \) だ。
# 小数の比較
import math
a = 0.3
b = 0.1 + 0.2
print("a == b : " + str(a == b))
print("isclose(a, b) = " + str(math.isclose(a, b)))
また、比較演算を伴わない場合は、組み込み関数の round関数を使って、小数点以下の適当な桁で丸めるのも解決策になる。プログラム "round1.py" では、小数第10位で丸めている。
# 小数の比較
import math
a = 0.3
b = 0.1 + 0.2
print("a = " + str(round(a, 10)))
print("b = " + str(round(b, 10)))

decimalモジュール

Python には、小数演算を10進数で行うことができる decimalモジュールが用意されている。これを使うことで、2進数の循環小数と縁を切ることができる。プログラム "decimal1.py" を実行してみてほしい。
# decimalモジュール
import decimal
a = decimal.Decimal("0.3")
b = decimal.Decimal("0.1") + decimal.Decimal("0.2")
print(f"{a:.54f}")
print(f"{b:.54f}")
2進数の循環小数が発生していないことが分かるだろう。
decimalモジュールは、decimal.Decimal("数字") のようにして使う。ここで、数値はかならす数字(文字列)として渡すこと。数値で渡してしまうと、decimalモジュールに渡る前に2進数の循環小数になってしまうためだ。
算術演算子や算術関数は、そのまま使うことができる。

floatが扱える数値の範囲

次に、Python で扱える数値の範囲に制約があるかどうかを見ていくことにしよう。

Python の小数(float)は、IEEE (アイ・トリプル・イー) という国際会議が定めた IEEE 754 規格(倍精度浮動小数点数)に則って小数を扱う。
浮動小数点数
浮動小数点数
IEEE 754 浮動小数点数
IEEE 754 浮動小数点数
IEEE 754 では、小数を符合部、仮数部、指数部の3つにわけ、符合部1ビット、仮数部52ビット、指数部11ビットの合計64ビットになる。このように符合部、仮数部、指数部に分けた小数を浮動小数点数と呼ぶ。つまり、Python で扱える小数の範囲は ±2.2×10-308~±1.8×10308 となる。
# float型の範囲
import sys
print(sys.float_info)
プログラム "float_info1.py" は sys.float_infoメソッドを使い、float が扱える範囲の情報を表示する。
# float型のオーバーフロー
import sys
a = sys.float_info.max
print(a * 10)

b = 0 - sys.float_info.max
print(b * 10)
プログラム "floatOverflow1.py" は、float が扱える最大値 sys.float_info.max を超えて計算したらどうなるかをみる。計算結果は inf(無限大)となり、オーバーフローしてしまう。負の最大値に対しても同様に計算しようとすると、-inf(負の無限大)となり、やはりオーバーフローしてしまう。
# float型のアンダーフロー
import sys
a = sys.float_info.min
print(a - 1)

b = 0 - sys.float_info.min
print(b + 1)
プログラム "floatUnderflow1.py" は、float が扱える最小値 sys.float_info.min を下回るとどうなるかをみる。計算結果は -1.0となり、sys.float_info.min が0とみなされアンダーフローしてしまう。負の最小値に対しても同様に計算しようとすると、+1.0 となり、やはりアンダーフローしてしまう。

decimalが扱える数値の範囲

decimalモジュール はどうだろうか――。
# decimalモジュールの有効桁数
import decimal
from decimal import getcontext
print(getcontext())
a = decimal.Decimal("1.0")
b = decimal.Decimal("3.0")
c = a / b
プログラム "decimal2.py" は、decimal が扱える最大桁数 getcontext().prec を表示する。デフォルトでは28桁だ(実行環境によって異なるかもしれない)。
10進数の循環小数 \( 1 \div 3 \) に対して print(f'{c:.30f}') を実行すると、小数点以下29桁目からゼロになることがわかる。
# decimalモジュールの有効桁数
import decimal
from decimal import getcontext
getcontext().prec = 48
a = decimal.Decimal("1.0")
b = decimal.Decimal("3.0")
c = a / b
プログラム "decimal3.py" を実行してみてほしい。
decimal の最大桁数 getcontext().prec は変更が可能で、このプログラムでは48桁に増やしている。print(f"{c:.50f}") を実行すると、小数点以下49桁目からゼロになる。

intが扱える数値の範囲

整数(int )はどうだろうか――Python では、メモリが許すかぎり(正負ともに)巨大な整数を扱うことができる。

巨大素数を求める

このように Python では扱える整数に上限がないことから、巨大素数を求めるプログラムを作ってみた。"getBigPrime.py" を実行してみてほしい。画面に200桁以下の素数が表示される。最大桁数は変数 MAX_DIGITS に代入しており、これを変更してみてほしい。桁数に制約はないが、計算量が幾何級数的に増えるので、計算結果が出るまでに時間がかかるようになる。
巨大素数を求める
巨大素数は暗号化の分野で無くてはならない数字だ。ここでは、扱える整数に上限がないという感触が得られれば十分である。
プログラムの解説は省略するが、興味がある方は「PHPとPythonで巨大素数を扱う」をご覧いただきたい。素数を求める手続きを説明している。

円周率を求める

Python では扱える小数の桁数に制約がないことから、円周率を求めるプログラム "pi.py" を作ってみた。画面に小数点以下第1000位までの円周率を表示する。
小数点以下の桁数は変数 MAX_DIGITS に代入しており、これを変更してみてほしい。桁数に制約はないが、計算量が幾何級数的に増えるので、計算結果が出るまでに時間がかかるようになる。
円周率を求める
ここではチュドノフスキーの公式を使って円周率を求めている。巨大素数と同じくプログラムの解説は省略するが、興味がある方は「C++ で円周率を計算する」をご覧いただきたい。

まとめ

今回学んだように、Python では、decimalを使えば2進数の誤差に惑わされることなく任意の有効桁数の小数演算ができ、巨大な整数(int)も扱うことができる。Python が科学技術計算に向いているというのは、こうした仕様によるところが大きい。
ただし、小数の有効桁数を大きくすればするほどメモリを消費し、演算速度は低下する。整数についても同じことが言える。残念ながら、Python他のインタプリタ型言語より処理速度が遅い
計算速度をとるか、有効桁数をとるかは、プログラム設計時に決めておく必要がある。

練習問題

次回予告

これまで見てきたように、Python ではダブルクォーテーション "..." またはシングルクォーテーション '...' で囲んだデータを文字列型(str)として扱う。次回は、基本的な文字列に加え、水平タブや改行といった表示制御文字、Unicode文字などを表すことができるエスケープシーケンス、整数にカンマ区切りを付けたり小数点以下の桁数の指定を行うなどの書式指定ができるフォーマット文字列書式演算子について学ぶ。

ループカウンタはなぜ「i」なのか

ジョン・バッカス
ジョン・バッカス
1954年(昭和29年)にIBMのジョン・バッカスが開発した世界初の高級言語「FORTRAN」は、変数名として6文字まで利用でき、データの型宣言が不要だった。とはいっても、動的型付けだったわけではない。変数名の頭文字がI~Nのいずれかだった場合、整数型として扱うというルールだった。

FORTRAN を考案したジョン・バッカスは、のちに、「Iから始まる6文字に決まったのは、みんな添え字にI、J、Kをつかっていたので、(L、M、Nを)気前よく増やした」と語っている。Nまで増やしたのには、INtegerというウィットを効かせたからとも言われている。

以来、多くのプログラミング言語で、(英小文字にはなったものの)ループカウンタにはi, j, kを用いるのが慣習となった。
(この項おわり)
header