今日は数値計算法の授業があったのだが、そこで bairstow 法なる手法が紹介された。これは多項式の解を複素数解まで含めてすべて近似的に求めることができる手法である。授業では「多項式を二次式で因数分解していく」「ニュートン法でいい感じにすると二次式の係数が求まる」ということを述べたのみで、後は写経の時間になった。
なるほど写経は概念を習得していない場合に有効だと思うけど私は概念を習得して自分で考えてプログラムを書きたかったので頭をひねったのだけど理解できなかった。プログラム系の授業で敗北したのは初めてだと思う。ちっぽけなプライドに傷をつけられたまま引き下がることはできないのできちんと理解してやる。
まず、対象にするのは n次の多項式だ。これの解をすべて求めることを目標にする。
この式(1) を適当な係数を持つ二次式 
で割るとこうなる。
こうすると の解は解の公式から求められる。n-2次の多項式は同じように2次ずつ次元を減らしていけば解けることになる、というのが bairstow 法のコンセプトだ。
ここで問題にするのは2つだ。一つは、どのような p, q を選べば余り（
）がないように割れるかということ。もう一つは、商の項の係数 
の値は何になるのかということ。
2つめの問題はかんたんなので先に取り組むことにする。 式(1)、(2) で 各項の係数を比較すれば良い。例えば、(2) 式中で n-2 次の項は 
である。これが、(1)式中での n-2次の項 
と恒に等しいので、
を導ける。添字を修正して
である。p, q がわかれば係数列 b を次数の高い方から求めることができる。ただし、
としておく。
続いて、余りを0にするためには何がどうなればよいかを考える。先程のように (1)、 (2) 式からなる恒等式を見れば、 
であり、 
である。(3) を用いて
を得る。
を 係数列bの項とp, qを使って表すことができた。係数列bの値も p, qが決まれば決まるので、 は二変数 p, q で決まる値である。 をそれぞれ0にするような p, qの値を、多変数のニュートン法によって求める。
少々脱線するがニュートン法についても解説しておく。ニュートン法は多項式の解の近似値を一つ求める手法である。
原理は省略するが、手法としては
の解 x' を求めるときには、
とする。nが大きいほど解としての精度はあがる。ちなみに 
の値は適当に決めておく。この値によっては異なる解が得られることがある。
二変数のニュートン法では 
を満たすような x, y を探す。
これについても原理を省いて手法のみ説明する。偏導関数 
が求まっていて、適当な 
があるとき、
とする。*1
以上が bairstow 法の理論である。理論がわかったので実装を行う。実装の方針は次の通り
- 数値微分を実装する
- 二変数のニュートン法を実装する
- 二次方程式の解の公式を実装する
- bairstow 法を実装する（
を求める関数を実装してニュートン方にかける ）
数値微分の実装
数値微分は、二変数のニュートン法において、偏微分係数を求めるために必要となる。実装は微分の定義式に従えば良い。
ただし実装にはより精度の良い三点差分による近似を用いた。導出は面倒そうなので省略するが直感的な理解があるので使って良いことにした。
D言語による実装ではこうなる。
auto differential(F)(F f, double h = float.epsilon)
{
return (double x) {
return (f(x + h) - f(x - h)) / (h * 2);
};
}
ここで 変数 hの値を小さくしすぎる（ double.epsilon など）と失敗する。細かい理由は追いかけていないがそうなりそうということはわかっているので良いことにした。
二変数のニュートンはそのまま実装する。関数名が気に食わないのでよりよい名前をご存知という方はお知らせ下さい。
import std.typecons;
Tuple!(double, double) newton2(F)(F f, F g, double x, double y, double epsilon=float.epsilon)
{
import std.math;
while (true) {
auto f_x = differential((double x) => f(x, y));
auto f_y = differential((double y) => f(x, y));
auto g_x = differential((double x) => g(x, y));
auto g_y = differential((double y) => g(x, y));
auto d = f_x(x) * g_y(y) - f_y(y) * g_x(x);
auto dx = (f(x, y) * g_y(y) - f_y(y) * g(x, y)) / d;
auto dy = (f_x(x) * g(x, y) - f(x, y) * g_x(x)) / d;
x = x - dx;
y = y - dy;
if (abs(dx) <= epsilon && abs(dy) <= epsilon) {
break;
}
}
return tuple(x, y);
}
先程定義した微分関数 differential を用いて偏導関数を導いている。 二変数の偏導関数の定義を書いておくと
である*2
やるだけ。これについても良い名前をご存じの方はお知らせ下さい。
import std.complex;
Complex!double[] solve2d(double a, double b, double c)
{
import std.math;
auto d = b * b - 4 * a * c;
if (d > 0) {
auto ans1 = (-b + sqrt(d)) / (2*a);
auto ans2 = (-b - sqrt(d)) / (2*a);
return [complex(ans1), complex(ans2)];
}
else if (d == 0) {
auto ans = (-b + sqrt(d)) / (2 * a);
return [complex(ans)];
}
else {
auto ans1 = complex(-b, sqrt(-d)) / (2 * a);
auto ans2 = ans1;
ans2.im = -ans2.im;
return [ans1, ans2];
}
}
bairstow法を実装する
実装はこうなる。
import std.complex;
Complex!double[] bairstow(double[] as, double epsilon = float.epsilon)
in {
assert(as.length > 1);
}
do
{
import std.range;
import std.typecons : tuple, Tuple;
import std.meta : AliasSeq;
auto n = as.length;
if (n == 2) {
return [complex(-as[0]/ as[1])];
}
if (n == 3) {
return solve2d(as[2], as[1], as[0]);
}
alias memo_key = Tuple!(ulong, double, double);
double[memo_key] memo;
double b_i(ulong i, double p, double q) {
auto key = tuple(i, p, q);
if (key in memo) { return memo[key]; }
if (i >= n - 2) { return 0; }
auto r = as[i+2] - b_i(i+1, p, q) * p - b_i(i+2, p, q) * q;
memo[key] = r;
return r;
}
auto alpha = (double p, double q) {
return as[1] - b_i(0, p, q) * p - b_i(1, p, q) * q;
};
auto beta = (double p, double q) {
return as[0] - b_i(0, p, q) * q;
};
double p = 1.0, q = 1.0;
AliasSeq!(p, q) = newton2(alpha, beta, p, q);
double[] bs = new double[](n-2);
foreach (i; 0..n-2) {
bs[i] = b_i(i, p, q);
}
return solve2d(1, p, q) ~ bairstow(bs, epsilon);
}
以下、各部分に分解して解説をはさむ。
一次式、二次式が渡されたときはそれぞれそのまま解を返している。
auto n = as.length;
if (n == 2) {
return [complex(-as[0]/ as[1])];
}
if (n == 3) {
return solve2d(as[2], as[1], as[0]);
}
係数列 b の各係数を求める関数は以下のとおりである。
alias memo_key = Tuple!(ulong, double, double);
double[memo_key] memo;
double b_i(ulong i, double p, double q) {
auto key = tuple(i, p, q);
if (key in memo) { return memo[key]; }
if (i >= n - 2) { return 0; }
auto r = as[i+2] - b_i(i+1, p, q) * p - b_i(i+2, p, q) * q;
memo[key] = r;
return r;
}
ニュートン法の中で、何度も呼ばれることになるため、メモ化を行っている。値を求める方法は愚直に (3) 式に従っている。
alpha, beta は b_i を用いて愚直に計算すれば良い。
auto alpha = (double p, double q) {
return as[1] - b_i(0, p, q) * p - b_i(1, p, q) * q;
};
auto beta = (double p, double q) {
return as[0] - b_i(0, p, q) * q;
};
こうして得られた alpha, beta をそれぞれ 0にするような p, qをニュートン法を用いて探索する。その後、 
の解を求め、残りの式を再度 bairstow 法で解くことで解を得る。
double p = 1.0, q = 1.0;
AliasSeq!(p, q) = newton2(alpha, beta, p, q);
double[] bs = new double[](n-2);
foreach (i; 0..n-2) {
bs[i] = b_i(i, p, q);
}
return solve2d(1, p, q) ~ bairstow(bs, epsilon);
動作例
試しに適当な多項式を解いてみる。ここでは授業中に課題として出された 
を解く*3。
void main()
{
writeln(bairstow([1, -3, 0, 0, 1]));
}
実行結果は次の通り
[-0.822576+1.26032i, -0.822576-1.26032i, 1.30749+0i, 0.337667+0i]
Wolfram Alpha に解かせてみると Real Solutions として x = 0.33767, x = 1.3075 、 Complex Solutionsとして x = -0.8226 +- 1.2603i を出してきたので多分正しい。
どうか、人々が bairstow 法に負けないように。
間違いを見つけたら教えてください。
