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

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