pythonで2次方程式の解の公式を解く

Question

前提・実現したいこと

pythonで二次方程式の解の公式を解いた時に出る誤差をなくしたいです。

発生している問題・エラーメッセージ

a=1.0,b=1.0000000000000001,c=0.0000000000000001としたときに
ax^2+bx+c=0
を解くと、x=-1,-10^(-16)となるはずですが
x1 = -0.9999999999999999
x2 = -1.0000000000000001e-16
となってしまいます。
これは何による誤差なのでしょうか。
修正する部分をを教えてほしいです。

補足情報（FW/ツールのバージョンなど）

該当のコードです

import math

a = float(input('aの値:  '))
b = float(input('bの値:  '))
c = float(input('cの値:  '))

D = b ** 2 - 4.0 * a * c

if D >= 0.0:
    x1 = (abs(b) + math.sqrt(D)) / (2.0 * a)
    if b > 0:
        x1 = -x1
    x2 = c / (a * x1)
    print('x1 = ', x1)
    print('x2 = ', x2)
else:
    print('実数解はありません')

Answer 1

問題があるのはbという変数です。以下のソースコードを試して下さい。

b = 1.0000000000000001
print("b={:.16f}".format(b))
# 結果
# b=1.00000000000000000

bの真の値は1+1E-16ですが1E-16の部分がなくなっています。これはコンピュータの小数の表し方が浮動小数点という方式を採用しているためです(現在ほとんどのコンピュータやプログラミング言語で採用されている方式はIEEE-754という浮動小数点方式になります)。

pythonでfloatは64bitsの倍精度になります(倍というのは単精度が32bitsですのでその倍という意味です)。簡単にいうと最上位桁の(bの例だと1)と最下位桁(bの例だと1E-16)は2進数で53桁までしか表現できない方式です。log10(2**53)はおよそ16弱になりますので10進数だと最上位の桁と最下位の桁は15桁までということになります。bは1.0と1E-16が16桁離れているので表現できずに誤差になるという訳です。

改善提案ですが64bitsの倍精度浮動小数点数を単純に128bitsの4倍精度浮動小数点を用いれば良いです。Pythonではnumpyという数値計算に必須のパッケージがありますのでこちらを使います。

import numpy as np
a = np.float128(input('aの値:  '))
b = np.float128(input('bの値:  '))
c = np.float128(input('cの値:  '))
# 以下同じ

他にも任意精度の浮動小数点ライブラリMPFRを用いたgmpy2パッケージや、なるべく二次元方程式に値を代入せずに変数のまま方程式として解くsimpyがあります。これらは用途に応じて使い分けることになります。

参考URL

• ウィキペディア IEEE 754
• docs.python.org 15. 浮動小数点演算、その問題と制限

Answer 2

浮動小数点数の有効数値は約15桁で、それ以上の桁は意味がありません。