はじめに
自分用メモです。
python3でシミュレーションをしていて問題に当たりました。問題が起きた部分を、簡単にすると以下のコードのようになります。
def f(a):
return 1 - (1 - a)**0.5
f(0.00000000001)
# 5.000000413701855e-12
f(0.000000000001)
# 5.000444502911705e-13
f(0.0000000000001)
# 5.007105841059456e-14
f(0.00000000000001)
# 4.9960036108132044e-15
f(0.000000000000001)
# 5.551115123125783e-16
f(0.0000000000000001)
# 0.0
f(0.00000000000000001)
# 0.0
数学的に言えば、関数fは引数aの値が0に向かって小さくなるにつれて、返り値が0に近くはずですが、aが0でない限り、結果は0にならないはずです。数値計算の限界として、ある程度aが小さくなった時に返り値が0になるのは仕方がないとして、aの値がたかだがe-17程度で結果が0になってしまっています。もうちょっと頑張ってもらえないか、というのが今回のテーマです。
decimalモジュールを使う
公式ドキュメントを見て解決しました。decimal モジュールを使えば良いようです。
decimalモジュールをimportした上で、有効桁数を指定します。
from decimal import *
getcontext().prec = 100
数値をDecimalオブジェクトに変換してから計算します。Decimalオブジェクトは、足し算引き算にもちいることはできますが、0.5乗とかはできず、0.5をDecimalオブジェクトに変換してから計算する必要があるようです。
def g(a):
return 1 - (1 - Decimal(a))**Decimal('0.5')
#実行
g(0.00000000000000000000000000000000001)
#Decimal('5.0000000000000000392877259729119015321129309725540314195553793115E-36')
ここで返り値はDicimalクラスのオブジェクトです。
Decimalからfloatに戻すには、float()に渡します。
def h(a):
x = 1 - (1 - Decimal(a))**Decimal('0.5')
return float(x)
h(0.0000000000000001)
#5e-17
h(0.000000000000000000000000000000000000001)
#5e-40
これで計算できるようになりました。
さっき設定したgetcontext().precはなんでしょうか？
getcontext()
#Context(prec=100, rounding=ROUND_HALF_EVEN, Emin=-999999, Emax=999999, capitals=1, clamp=0, flags=[Inexact, FloatOperation, Rounded], traps=[InvalidOperation, DivisionByZero, Overflow])
#タイプを調べます。
type(getcontext())
#<class 'decimal.Context'>
prec以外の設定もここから変えられるようです。
意識していませんでしたが、浮動小数値を扱う時に、仮数部が長いと丸められてしまいます。
#ぎりぎりアウト
a = 0.99999999999999999
print(a)
#1.0
#ぎりぎりセーフ
b = 0.9999999999999999
print(b)
# 0.9999999999999999
Decimalに渡すときには注意が必要のようです。文字列で渡してあげないと、丸められたものが渡されてしまうようです。
Decimal(0.9999999999999999999999999)
#Decimal('1')
Decimal('0.9999999999999999999999999')
#Decimal('0.9999999999999999999999999')
終わりに
Decimalは奥が深そうですが、一応計算できるようになったと思われますのでこの記事はこれで終わりにしたいと思います。