二次方程式の解を求める - Python の sympy の場合

「二次方程式の解を求めるプログラムのレビュー」で，実数解の場合に解の公式だけで解を求めてはならないと述べたが，Python の sympy がどうやっているのかレビューした結果を示す。

以下の二次方程式を解け。
x**2 +1.0000000000000001*x + 0.0000000000000001 = 0

結論から述べると，正しいやり方は「二次方程式の係数を Rational() で与え，solve() で解を求める」ということである。

>>> from sympy import *
>>> var('a:z')
(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z)

係数を Rational() で定義する。

>>> b = Rational(10000000000000001, 10000000000000000)
>>> c = Rational( 1, 10000000000000000)
>>> print(b, c)
10000000000000001/10000000000000000 1/10000000000000000

Rational() ならば精度は完全に保証される。

>>> print(b.evalf(30))
1.00000000000000010000000000000
>>> print(c.evalf(30))
1.00000000000000000000000000000e-16

solve() で解を求める

>>> ans = solve(x**2 + b * x + c)

>>> print(ans)
[-1, -1/10000000000000000]

つまり，

>>> print(ans[0].evalf(30))
-1.00000000000000000000000000000
>>> print(ans[1].evalf(30))
-1.00000000000000000000000000000e-16

以下は読まなくて良い。

ダメなやり方

python の float を Float() で変換しても，evalf() で有効数字を 30 桁表示させれば精度が足りないのは明らかである。

>>> b = Float(1.0000000000000001)
>>> c = Float(0.0000000000000001)
>>> print(b.evalf(30), c.evalf(30))
1.00000000000000000000000000000 9.99999999999999979097786724035e-17

文字型で指定してやれば，精度は増す。

>>> b = Float('1.00000000000000010000000000000000000000000000')
>>> c = Float('0.00000000000000010000000000000000000000000000')
>>> print(b.evalf(30), c.evalf(30))
1.00000000000000010000000000000 1.00000000000000000000000000000e-16

十分（？）な精度の b, c で二次方程式を解くと以下のようになる。

>>> ans = solve(x**2 + b * x + c, x)

デフォルトでの表示は正確な解が得られたと見まがう。

>>> print(ans)
[-1.00000000000000, -1.00000000000000e-16]

しかし，個別の解を .evalf(30) で表示させると，精度が足りないことが明白になる。

>>> print(ans[0].evalf(30))
-0.999999999999999888977697537484
>>> print(ans[1].evalf(30))
-1.00000000000000010235730316482e-16

解の公式だけで実数解を求める場合には，大小どちらの解も精度は十分であることがわかる。（何故 solve() の場合と違うの か？）

>>> x1 = (-b - sqrt(b**2 - 4 * c)) / 2
>>> print(x1.evalf(30))
-1.00000000000000000000000000000
>>> x2 = (-b + sqrt(b**2 - 4 * c)) / 2
>>> print(x2.evalf(30))
-1.00000000000000000000000000000e-16