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

「二次方程式の解を求めるプログラムのレビュー」で，実数解の場合に解の公式だけで解を求めてはならないと述べた。

前の記事で，Python の sympy がどうやっているのかレビューした結果を示した。

今回は，R において，sympy 相当の解を求める方法を示す。

sympy の Rational() に相当するのは，gmp ライブラリにおける，bigq クラスの小数である。しかし，許されている演算は，四則演算，比較演算，abs() 程度である。四則演算の中には ^ もあるが，平方根を取る '^0.5' はインプリメントされていない。そこで，以下のような Sqrt() 関数を自分で定義して使わなければならない。

> library(gmp)

> # ニュートン・ラフソン法で bigq クラスの小数の平方根を求める関数
> Sqrt = function(n) {
+   n = as.bigq(n)
+   x = as.bigq(1)
+   EPSILON = as.bigq(x, 1e20)
+   while (TRUE) {
+     x2 = (x + n / x) / 2
+     if (abs(x - x2) <= EPSILON) {
+       return(x2)
+     }
+     x = x2
+   }
+ }


ちゃんと動くか試しておこう
	
> a = Sqrt(1000)
> print(as.numeric(a*a))
[1] 1000

x^2 + 1.0000000000000001*x + 0.0000000000000001 = 0 を解く

> b = as.bigq("10000000000000001", "10000000000000000")
> c = as.bigq( "1", "10000000000000000")	
> x1 = (-b - Sqrt(b^2 - 4*c)) / 2
> print(x1)
Big Rational ('bigq') :
[1] -159999999999999984000000000000000800000000000000000000000000000001/159999999999999984000000000000000800000000000000000000000000000000
> print(as.numeric(x1))
[1] -1
	
> x2 = (-b + Sqrt(b^2 - 4*c)) / 2
> print(x2)
Big Rational ('bigq') :
[1] -15999999999999998400000000000000079999999999999999/159999999999999984000000000000000800000000000000000000000000000000
> print(as.numeric(x2))
[1] -0.0000000000000001
	
> x3 = c / x1
> print(x3)
Big Rational ('bigq') :
[1] -15999999999999998400000000000000080000000000000000/159999999999999984000000000000000800000000000000000000000000000001
> print(as.numeric(x3))
[1] -0.0000000000000001

解の公式だけで両方の実数解を正しく求めることができた。 当然ではあるが，大きい方の実数解を解の公式で求めて，小さい方の実数解は解と係数の関係から求めるというやり方でも正しく求める事が出来る。

二次方程式の解を求める - 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