はじめに
前回の記事からの続きである。前回の記事では実数係数の三次方程式の実数解を得るところまで辿りつけた。今回は目標の四次方程式の実数解を得るところまで到達したい。
四次方程式の解法
四次方程式の解法として最初に発見された手法は発見者にちなんでフェラーリの解法と言われる。名前自体がカッコイイので基本的にそのまま踏襲したい。
チルンハウス変換
$x$ に関する四次方程式の一般形式を以下に示す。係数 $A$ ～ $E$ はすべて実数とする。
\def\dfrac{\displaystyle\frac}
Ax^4 + Bx^3 + Cx^2 + Dx + E = 0
\tag{1}
最大次数の係数 $A \ne 0$ として，以下のように最大次数の係数を１に正規化したところからスタートしても構わないだろう。
x^4 + \underbrace{\frac{B}{A}}_ax^3 + \underbrace{\frac{C}{A}}_bx^2 + \underbrace{\frac{D}{A}}_cx + \underbrace{\frac{E}{A}}_d = 0
\tag{2}
さらに $y = x + a/4$ とおいて三次の係数をゼロにする。※チルンハウス変換と呼ぶらしい。
\left(y - \frac{a}{4}\right)^4
+ a\left(y - \frac{a}{4}\right)^3
+ b\left(y - \frac{a}{4}\right)^2
+ c\left(y - \frac{a}{4}\right) + e = 0
\begin{align}
y^4 -\cancel{ay^3} + \frac{6a^2}{16}y^2 - \frac{a^3}{16}y + \frac{a^4}{256} \\
+ a\left(\cancel{y^3} - \frac{3a}{4}y^2 + \frac{3a^2}{16}y - \frac{a^3}{64} \right) \\
+ b\left(y^2 - \frac{a}{2}y + \frac{a^2}{16}\right) \\
+ c\left(y - \frac{a}{4}\right) + e = 0
\end{align}
y^4 + \underbrace{\left(- \frac{3a^2}{8} + b\right)}_py^2 + \underbrace{\left(\frac{a^3}{8} - \frac{ab}{2} + c\right)}_q y + \underbrace{\left(- \frac{3a^4}{256} + \frac{a^2b}{16} - \frac{ac}{4} + d\right)}_r = 0
こうして以下のようにすっきりした状態からスタートすることにする。
y^4 + py^2 + q y + r = 0
\tag{3}
複二次式として解ける場合
式 $(3)$ において $q = 0$ のとき，複二次式として簡単に解ける。$t = y^2$ とおくと
t^2 + pt + r = 0
\tag{4}
と二次方程式になるからだ。判別式 $m = p^2 - 4r$ とおいて，$m < 0$ のときは $t$ の実数解が一つも存在しない。$t$ が実数でないときは $y$ も実数にならないので考えなくてよい。$m \ge 0$ のとき $t$ の二つの実数解をそれぞれ $u$，$v$ とおくと
\left\{\begin{align}
u &= \frac{-p + \sqrt{m}}{2} \\
v &= \frac{-p - \sqrt{m}}{2}
\end{align}\right.
\tag{5}
となる。さらに
y = \left\{\begin{align}
\pm\sqrt{u} & \quad (u \ge 0) \\
\pm\sqrt{v} & \quad (v \ge 0)
\end{align}\right.
\tag{6}
として最大４個の実数解を得る。
ここで $q = 0$ の条件を取り除いておかないと，次に紹介するフェラーリの解法では解けない場合が生じるので注意すること。
フェラーリの解法
式 $(3)$ において $q \ne 0$ のときである。式 $(3)$ の両辺に $ty^2 + t^2/4$ を加える。
y^4 + ty^2 + \frac{t^2}{4} = (t - p)y^2 - q y + \frac{t^2}{4} - r
\left(y^2 + \frac{t}{2}\right)^2 = (t - p)\left\{y^2 - \frac{q}{t - p} y + \frac{t^2 - 4r}{4(t - p)} \right\}
\tag{7}
ここで下記の関係が成立したとする。
\frac{t^2 - 4r}{4(t - p)} = \frac{q^2}{4(t - p)^2} = \left\{\frac{q}{2(t - p)}\right\}^2
\tag{8}
こうすると式 $(7)$ の右辺も平方完成する。
\left(y^2 + \frac{t}{2}\right)^2 = (t - p)\left(y - \frac{q}{2(t-p)} \right)^2
\tag{9}
これより，
\left(y^2 - y \sqrt{t - p} + \frac{t}{2} + \frac{q}{2\sqrt{t - p}}\right)\!\left(y^2 + y \sqrt{t - p} + \frac{t}{2} - \frac{q}{2\sqrt{t - p}}\right) = 0
\tag{10}
と二次式の積に因数分解できる。このような条件を満たす $t$ は式$(8)$を変形して得られる下記の三次方程式を解けば求められる。
(t - p)(t^2 - 4r) - q^2 = 0
\tag{11}
さて，本記事ではあくまで実数解を得たいので $t - p > 0$ になって貰わなくては困るが，式 $(11)$ の三次方程式において $q \ne 0$ の場合は $p$ より大きい実数解を持つことを簡単に示すことができる。式 $(11)$ より
f(t) = (t - p)(t^2 - 4r) - q^2
\tag{12}
とおくと，$f(p) = -q^2 < 0$ となり，必ず $t > p$ において $f(t) = 0$ となる $t$ が存在することが分かる。また逆に $q = 0$ のときは $t - p = 0$ となる可能性がある（または $t^2 - 4r = 0$ の可能性もある）ため，事前に排除しておかないと $t - p = 0$ となってゼロ除算エラーに陥る場合がある。
式 $(11)$ の三次方程式を解くと最大３つの実数解を得られるが，このうち最大値を用いるようにしたほうが良い。ちなみに筆者の前記事で作成した三次方程式の実数解ソルバーは得られた複数の実数解を配列に格納して返すが，配列の先頭が必ず最大値になっている。別にソートしている訳ではなく，偶々そのような特性になっている。
式 $(10)$ の二つの二次式の判別式をそれぞれ $u$，$v$ とおく。
\left\{\begin{align}
u &= t - p - 4\left(\frac{t}{2} + \frac{q}{2\sqrt{t - p}}\right)
= -t - p - \frac{2q}{\sqrt{t - p}}\\
v &= t - p - 4\left(\frac{t}{2} - \frac{q}{2\sqrt{t - p}}\right)
= -t - p + \frac{2q}{\sqrt{t - p}}
\end{align}\right.
\tag{13}
４個の実数解はぞれぞれ
y = \left\{\begin{align}
\dfrac{ \sqrt{t - p} \pm \sqrt{u}}{2} & \quad (u \ge 0) \\
\dfrac{-\sqrt{t - p} \pm \sqrt{v}}{2} & \quad (v \ge 0)
\end{align}\right.
\tag{14}
となる。さらに
\left\{\begin{align}
m &= \sqrt{t - p} \\
w &= \frac{q}{\sqrt{t - p}}
\end{align}\right.
\tag{15}
とおくと二つの判別式 $u$，$v$ はそれぞれ
\left\{\begin{align}
u &= -t - p - 2w\\
v &= -t - p + 2w
\end{align}\right.
\tag{16}
と書け，４個の実数解はぞれぞれ
y = \left\{\begin{align}
\dfrac{ m \pm \sqrt{u}}{2} & \quad (u \ge 0) \\
\dfrac{-m \pm \sqrt{v}}{2} & \quad (v \ge 0)
\end{align}
\right.
\tag{17}
となる。
ここで判別式 $u,\ v$ および実数解の算出式を中間変数 $m,\ w$ を用いて書き直したのは後述の形式と合わせるためである。
精度劣化対策
解析的に解くのであれば前項までの解法そのままでよい。ところが実際にプログラムに実装すると思わぬ精度劣化に悩む場合がある。それは浮動小数点同士の減算を行うとき，二つの値が極めて近接していると多くの有効数字を失ってしまうことだ。
具体的に言うと上記の $\sqrt{t - p}$ に当たる。$t$ と $p$ が極めて近接しているとき，$\sqrt{t - p}$ は
多くの有効数字を失ってしまう。幸いなことに式 $(11)$ より
\sqrt{t - p} = \frac{|q|}{\sqrt{t^2 - 4r}}
\tag{18}
という関係があるので置き換えることができる。これより
\left\{\begin{align}
m &= \sqrt{t - p} = \frac{|q|}{\sqrt{t^2 - 4r}}\\
w &= \frac{q}{m} = \frac{q}{|q|} \sqrt{t^2 - 4r}
\end{align}\right.
\tag{19}
となる。これをそのまま用いても良いが，もう少しだけ計算の最適化を図ろう。
まず $q > 0$ のとき
\left\{\begin{align}
m &= \frac{q}{\sqrt{t^2 - 4r}}\\
w &= \sqrt{t^2 - 4r}
\end{align}\right.
\tag{20}
となるから４つの解
y = \left\{\begin{align}
\dfrac{m \pm \sqrt{-t - p - 2w}}{2}
&= \dfrac{q}{2\sqrt{t^2 - 4r}} \pm \dfrac{\sqrt{-t - p - 2\sqrt{t^2 - 4r}}}{2}\\
\dfrac{-m \pm \sqrt{-t - p + 2w}}{2}
&= -\dfrac{q}{2\sqrt{t^2 - 4r}} \pm \dfrac{\sqrt{-t - p + 2\sqrt{t^2 - 4r}}}{2}
\end{align}\right.
\tag{21}
を得る。一方 $q < 0$ のとき
\left\{\begin{align}
m &= -\frac{q}{\sqrt{t^2 - 4r}}\\
w &= -\sqrt{t^2 - 4r}
\end{align}\right.
\tag{22}
となるから，４つの解
y = \left\{\begin{align}
\dfrac{m \pm \sqrt{-t - p - 2w}}{2}
&= -\dfrac{q}{2\sqrt{t^2 - 4r}} \pm \dfrac{\sqrt{-t - p + 2\sqrt{t^2 - 4r}}}{2} \\
\dfrac{- m \pm \sqrt{-t - p + 2w}}{2}
&= \dfrac{q}{2\sqrt{t^2 - 4r}} \pm \dfrac{\sqrt{-t - p - 2\sqrt{t^2 - 4r}}}{2}
\end{align}\right.
\tag{23}
を得る。ここで $q$ の正負にかかわらず，共通の解を得ていることに気づく。
y = \left\{\begin{align}
-\dfrac{q}{2\sqrt{t^2 - 4r}} \pm \dfrac{\sqrt{-t - p + 2\sqrt{t^2 - 4r}}}{2} \\
\dfrac{q}{2\sqrt{t^2 - 4r}} \pm \dfrac{\sqrt{-t - p - 2\sqrt{t^2 - 4r}}}{2}
\end{align}\right.
\tag{24}
ここで改めて
\left\{\begin{align}
m' &= \frac{q}{\sqrt{t^2 - 4r}}\\
w' &= \sqrt{t^2 - 4r}
\end{align}\right.
\tag{25}
とおくと $q$ の正負にかかわらず，４つの実数解を
y = \left\{\begin{align}
\dfrac{-m' \pm \sqrt{-t - p + 2w'}}{2} \\
\dfrac{ m' \pm \sqrt{-t - p - 2w'}}{2}
\end{align}\right.
\tag{26}
と書ける。
解法まとめ
以上をまとめると下記のようになる。チルンハウス変換以降のみを示す。なお，根号の中身が負になる場合は計算しないうえ，その結果を用いる計算も行わない。
| $q = 0$ のとき | $q \ne 0$ のとき | ||
|---|---|---|---|
| $t - p > t^2 - 4r$ のとき | $t - p \le t^2 - 4r$ のとき | ||
| 中間変数① | $m = p^2 - 4r$ | $m = \sqrt{t - p}$ | $m = \dfrac{q}{\sqrt{t^2 - 4r}}$ |
| 中間変数② | $w = \dfrac{q}{\sqrt{t - p}}$ | $w = \sqrt{t^2 - 4r}$ | |
| 判別式① | $u = \dfrac{-p + \sqrt{m}}{2}$ | $u = -t - p - 2w$ | $u = -t - p - 2w$ |
| 判別式② | $v = \dfrac{-p - \sqrt{m}}{2}$ | $v = -t - p + 2w$ | $v = -t - p + 2w$ |
| 判別式① による解 |
$y = \pm\sqrt{u}$ | $y = \dfrac{m \pm \sqrt{u}}{2}$ | $y = \dfrac{m \pm \sqrt{u}}{2}$ |
| 判別式② による解 |
$y = \pm\sqrt{v}$ | $y = \dfrac{-m \pm \sqrt{v}}{2}$ | $y = \dfrac{-m \pm \sqrt{v}}{2}$ |
$t - p$ と $t^2 - 4r$ の大小によって計算を分岐せざるを得ないが，中間変数 $m$，$w$ の採用により計算の共通化を図り，再度合流できるようにした。
実装コード
実装コードを以下に示す。四次方程式のソルバー内で三次方程式のソルバーを使うので，前回の記事の成果をまるまる盛り込んでいる。
//------------------------------------------------------------------------------
// 四次方程式ソルバーにて複二次式と判断するための閾値：|q|＜ε
//------------------------------------------------------------------------------
var EPS = 1e-12;
//------------------------------------------------------------------------------
// 実数の立方根を求める
//------------------------------------------------------------------------------
function cbrt(x) {
if(x >= 0) {
return Math.pow( x, 1 / 3);
} else {
return -Math.pow(-x, 1 / 3);
}
}
//------------------------------------------------------------------------------
// 実数係数の三次方程式の実数解を得る
//------------------------------------------------------------------------------
function cubic_equation(aa, bb, cc, dd) {
var a = bb / aa;
var b = cc / aa;
var c = dd / aa;
var p = b - a * a / 3;
var q = 2 * a * a * a / 27 - a * b / 3 + c;
var d = q * q / 4 + p * p * p / 27;
var ret = [];
if(d > 0) {
var u = cbrt(-q / 2 + Math.sqrt(d));
var v = cbrt(-q / 2 - Math.sqrt(d));
ret.push(u + v - a / 3);
} else {
var th = Math.atan2(Math.sqrt(-d), -q / 2);
ret.push(2 * Math.sqrt(-p / 3) * Math.cos(th / 3) - a / 3);
ret.push(2 * Math.sqrt(-p / 3) * Math.cos(th / 3 - 2 * Math.PI / 3) - a / 3);
ret.push(2 * Math.sqrt(-p / 3) * Math.cos(th / 3 + 2 * Math.PI / 3) - a / 3);
}
return ret;
}
//------------------------------------------------------------------------------
// 実数係数の四次方程式の実数解を得る
//------------------------------------------------------------------------------
function quartic_equation(aa, bb, cc, dd, ee) {
var a = bb / aa;
var b = cc / aa;
var c = dd / aa;
var d = ee / aa;
var p = -3 * a * a / 8 + b;
var q = a * a * a / 8 - b * a / 2 + c;
var r = -3 * a * a * a * a / 256 + b * a * a / 16 - c * a / 4 + d;
var ret = [];
if(Math.abs(q) < EPS) {
var m = p * p - 4 * r;
if(m >= 0) {
var u = (-p + Math.sqrt(m)) / 2;
var v = (-p - Math.sqrt(m)) / 2;
if(u >= 0) {
ret.push( Math.sqrt(u) - a / 4);
ret.push(-Math.sqrt(u) - a / 4);
}
if(v >= 0) {
ret.push( Math.sqrt(v) - a / 4);
ret.push(-Math.sqrt(v) - a / 4);
}
}
} else {
var e = cubic_equation(1, -p, -4 * r, 4 * p * r - q * q);
var t = e[0];
if(t - p > t * t - 4 * r) {
var m = Math.sqrt(t - p);
var w = q / m;
} else {
var w = Math.sqrt(t * t - 4 * r);
var m = q / w;
}
var u = -t - p + 2 * w;
var v = -t - p - 2 * w;
if(u >= 0) {
ret.push((-m + Math.sqrt(u)) / 2 - a / 4);
ret.push((-m - Math.sqrt(u)) / 2 - a / 4);
}
if(v >= 0) {
ret.push((m + Math.sqrt(v)) / 2 - a / 4);
ret.push((m - Math.sqrt(v)) / 2 - a / 4);
}
}
return ret;
}
テスト用コード
テスト用コードを示す。先ほどの real_solver.js を参照しているので同じフォルダに置くこと。
なお，得られた実数解に加えて，得られた実数解を元の方程式に代入した検証結果も出力するようにした。理想的にはもちろんゼロになることが望ましい。
<package>
<job>
<script src="real_solver.js"></script>
<script>
var args = WScript.Arguments.Unnamed;
var ret = main(args);
try {
WScript.Quit(ret);
} catch(e) {
/* 何もしない */
}
function main(args) {
if(args.Count < 5) {
WScript.StdErr.WriteLine("実数係数の四次方程式の実数解を求めます。");
WScript.StdErr.WriteLine("");
WScript.StdErr.WriteLine("TEST(.WSF) [A] [B] [C] [D] [E]");
WScript.StdErr.WriteLine("");
WScript.StdErr.WriteLine("四次方程式の係数を以下のように定義します。");
WScript.StdErr.WriteLine("A･x^4 + B･x^3 + C･x^2 + D･x + E = 0");
return -1;
}
var k = [];
for(var i = 0; i < 5; i++)
k.push(parseFloat(args(i)));
if(k[0] == 0) {
WScript.StdErr.WriteLine("四次の係数は非ゼロにして下さい!!");
return -1;
}
var ret = quartic_equation(k[0], k[1], k[2], k[3], k[4]);
for(var i = 0; i < ret.length; i++) {
var x = ret[i];
var y = k[0] * x * x * x * x
+ k[1] * x * x * x
+ k[2] * x * x
+ k[3] * x
+ k[4];
WScript.Echo([x, y].join("\t"));
}
return 0;
}
</script>
</job>
</package>
テスト結果
引数なし，または引数が４個に満たない場合はヘルプメッセージを表示する。
ｃ:\Qiita>test
実数係数の四次方程式の実数解を求めます。
TEST(.WSF) [A] [B] [C] [D] [E]
四次方程式の係数を以下のように定義します。
A･x^4 + B･x^3 + C･x^2 + D･x + E = 0
四次の係数にゼロを与えるとエラーになる。
ｃ:\Qiita>test 0 0 0 0 0
四次の係数は非ゼロにして下さい!!
複二次式で解けるルート：$q = 0$ のとき
$x^4 -5y^2 + 4 = (x^2 - 4)(x^2 - 1) = (x - 2)(x + 2)(x - 1)(x + 1)$
$x = \pm1,\ \pm2$
ｃ:\Qiita>test 1 0 -5 0 4
2 0
-2 0
1 0
-1 0
フェラーリの解法の通常ルート：$q \ne 0$，$t - p > t^2 - 4r$ のとき
$x^4 -6x^2 + 3x + 2 = \left(x^2 + \dfrac{3}{2}\right)^2 - 3^2\left(x - \dfrac{1}{6}\right)^2 = \left(x^2 - 3x + 2\right)\left(x^2 + 3x + 1\right)$
$x = \dfrac{-3 \pm \sqrt{5}}{2},\ 1,\ 2$
$(t,p,q,r) = (3, -6, 3, 2)$ より $t - p = 9$，$t^2 - 4r = 1$
ｃ:\Qiita>test 1 0 -6 3 2
-0.3819660112501051 4.440892098500626e-16
-2.618033988749895 -5.329070518200751e-15
2 0
1 0
精度劣化対策版のルートその１：$q > 0$，$t - p \le t^2 - 4r$ のとき
$x^4 + 3x - 2 = \left(x^2 + \dfrac{1}{2}\right)^2 - 1^2\left(x - \dfrac{3}{2}\right)^2 = (x^2 - x + 2)(x^2 + x - 1)$
$x = \dfrac{-1 \pm \sqrt{5}}{2}$
$(t,p,q,r) = (1, 0, 3, -2)$ より $t - p = 1$，$t^2 - 4r = 9$
ｃ:\Qiita>test 1 0 0 3 -2
0.6180339887498949 0
-1.618033988749895 0
精度劣化対策版のルートその２：$q < 0$，$t - p \le t^2 - 4r$ のとき
$x^4 - 3x - 2 = \left(x^2 + \dfrac{1}{2}\right)^2 - 1^2\left(x + \dfrac{3}{2}\right)^2 = (x^2 - x - 1)(x^2 + x + 2)$
$x = \dfrac{1 \pm \sqrt{5}}{2}$
$(t,p,q,r) = (1, 0, -3, -2)$ より $t - p = 1$，$t^2 - 4r = 9$
ｃ:\Qiita>test 1 0 0 -3 -2
1.618033988749895 0
-0.6180339887498949 0
失敗及び隘路事項
完成したと思ってテスト運用に入っていたところ，得られた四次方程式の解を元の方程式に代入してもゼロにならないという不具合事例を散見した。一例を示すと下記のような方程式だ。
$x^4 + 485x^3 + 59200x^2 + 84400x -123600000000 = 0$
開発当初のバージョンで実行すると下記の結果となった。得られた実数解は真の値にそれなりに近似しているが，何よりも得られた解を元の方程式に代入してもゼロにならない。
c:\Qiita>test 1 485 59200 84400 -12360000000
233.28071342832106 -118.06695938110351
-475.73086001324475 -118.06695747375488
これに対し，精度劣化対策版の結果を示す。得られた解を元の方程式に代入すると，ほぼゼロとみなせるようになった。
c:\Qiita>test 1 485 59200 84400 -12360000000
233.28071417716222 -0.000003814697265625
-475.73086076232397 -0.0000095367431640625
ちなみに高精度計算サイトの結果も併せて示す。不一致部分を朱記した。
\begin{array}{|c|c|c|}
\hline
\textsf{ソルバー} & \textsf{実数解①} & \textsf{実数解②} \\
\hline
\textsf{当初バージョン} & 233.28071\color{red}{342832106} & -475.730860\color{red}{01324475} \\
\hline
\textsf{精度劣化対策版} & 233.2807141771622\color{red}{2} & -475.73086076232\color{red}{397} \\
\hline
\textsf{高精度計算サイト} & 233.28071417716224 & -475.73086076232402 \\
\hline
\end{array}
精度劣化の要因は本文でも述べたように $t$ と $p$ が極めて近い値になり，$t - p$ の有効数字が大幅に失われたことにある。上記の例の場合を示すと下記のように10進数で上位７桁まで一致しており，$t - p$ はその上位７桁がゼロになって一気に有効数字を失ったのだ。
\left\{\begin{align}
t &= -29009.372514637005 \\
p &= -29009.375
\end{align}\right.
参考文献