二次方程式を習うのって中学三年らしい・・・
はじめに
中学で二次方程式を習ったとき，繰り返しの計算練習に飽き飽きした記憶はないだろうか？二次方程式があるのなら三次方程式，四次方程式もあるはずで，これらの解法（解の公式）がどんなものなのか興味本位で調べた人も多いだろう。
続けて五次方程式には解の公式が存在しないこと，さらに欲張ってガロア理論まで学ぼうとして挫折した人もいるだろう・・・自分も同類なので安心されたい。
課題の再定義
仕事で三次方程式あるいは四次方程式を解く必要が生じた人っているだろうか？ちなみに筆者はある・・・というか今だ。
係数が具体的に与えられた方程式を解くだけなら，そのようなサービスを提供している WEB サイトがあると思うので，それを利用すればいい。Python などにもそのようなパッケージがあると思う。
問題は自社ツール内に三次方程式や四次方程式のソルバーを組み込みたいとき，仕様や精度・権利関係のよく分からない外部ライブラリに依存したくない，ささっと書ける分量なら自前で実装したいという人も少なくはないだろう。
ということで，今回のお題は「実数係数の三次方程式の実数解を得るソルバー関数を自作する」ということである。プログラミング言語は Windows を買うとタダで付いてくる WSH / JavaScript とする。
言語ネイティブに複素数型をサポートしていないのがポイントである。
三次方程式の定義
$x$ に関する三次方程式の一般形式を以下に示す。係数 $A$ ～ $D$ はすべて実数とする。
Ax^3 + Bx^2 + Cx + D = 0 \qquad(A \ne 0)
\tag{1}
とはいえ最大次数の係数 $A \ne 0$ として，以下の形式からスタートしても問題ないだろう。
x^3 + \underbrace{\frac{B}{A}}_ax^2 + \underbrace{\frac{C}{A}}_bx + \underbrace{\frac{D}{A}}_c = 0
x^3 + ax^2 + bx + c = 0
\tag{2}
チルンハウス変換
三次方程式の解法の多くが二次の係数をゼロにするように変数変換を行う。一般的にはチルンハウス変換と呼ぶらしい。具体的には式$(2)$に $x = y - a / 3$ と代入する。実際に確かめてみよう。
\left(y - \frac{a}{3}\right)^3 + a\left(y - \frac{a}{3}\right)^2 + b\left(y - \frac{a}{3}\right) + c = 0
y^3 - \cancel{ay^2} + \frac{a^2}{3} y - \frac{a^3}{27} + \cancel{ay^2} - \frac{2a^2}{3} y + \frac{a^2}{9} + by - \frac{ab}{3} + c = 0
y^3 + \underbrace{\left(b - \frac{a^2}{3}\right)}_py + \underbrace{\frac{2a^2}{27} - \frac{ab}{3} + c}_q = 0
y^3 + py + q = 0 \tag{3}
こうして変換された式$(3)$を改めてスタート位置として再定義する。
実数解の個数と判別式
実数係数の三次方程式は少なくとも一つ以上の実数解を持つ。重解の場合を除き，実数解の個数は１個，もしくは３個のいずれかになり，実数解２個と複素数解１個のような組み合わせは起こらない。
後で詳しく述べるが，実数解の個数に応じて解法を切り替えているので，できるだけ早く，できるだけ簡便な判別式で実数解の個数を得たい。
まず三次関数 $f(y)$ を以下のように定義する。
f(y) = y^3 + py + q
\tag{4}
導関数 $f'(y)$ は二次式となる。
f'(y) = 3y^2 + p
\tag{5}
導関数 $f'(y) = 0$ となる $y$ を $\alpha$，$\beta$ とおく。なお $\alpha < \beta$ とする。
\begin{array}{c|c|c|c|c|c}
y & \cdots & \alpha & \cdots & \beta & \cdots \\
\hline
f'(y) & + & 0 & - & 0 & + \\
\hline
f(y) & \nearrow & \textsf{極大} & \searrow & \textsf{極小} & \nearrow
\end{array}
三次方程式 $f(y) = 0$ が相異なる３つの実数解を持つ条件は $f(\alpha) > 0$ かつ $f(\beta) < 0$ である。これを一つの条件式にまとめると $f(\alpha) f(\beta) < 0$ となる。
\begin{align}
f(\alpha) f(\beta) &= \left(\alpha^3 + p\alpha +q\right) \cdot \left(\beta^3 + p\beta + q\right) \\
&= \alpha^3 \beta^3 + p\alpha^3\beta + q\alpha^3 + p\alpha\beta^3 + p^2\alpha\beta + pq\alpha + q\beta^3 + pq\beta + q^2 \\
&= \alpha^3\beta^3 + q(\alpha^3 + \beta^3) + p\alpha\beta(\alpha^2 + \beta^2) + p^2\alpha\beta + pq(\alpha + \beta) + q^2
\end{align}
\tag{6}
解と係数の関係より，
\left\{\begin{align}
\alpha + \beta &= 0 \\
\alpha \beta &= \frac{p}{3}
\end{align}\right.
\tag{7}
であるから
\left\{\begin{align}
\alpha^2 + \beta^2 &= (\underbrace{\alpha + \beta}_0)^2 - 2\alpha\beta = -2\alpha\beta = -\frac{2p}{3} \\
\alpha^3 + \beta^3 &= (\underbrace{\alpha + \beta}_0)^3 - 3\alpha\beta(\underbrace{\alpha + \beta}_0) = 0
\end{align}\right.
\tag{8}
となるので，
\begin{align}
f(\alpha) f(\beta) &= \left(\frac{p}{3}\right)^3 + p\left(\frac{p}{3}\right)\left(-\frac{2p}{3}\right) + p^2\left(\frac{p}{3}\right) + q^2 \\
&= \frac{1 - 6 + 9}{27}\cdot p^3 + q^2 \\
&= q^2 + \frac{4p^3}{27} < 0
\end{align}
\tag{9}
となる。判別式 $D = q^2 + 4p^3/27$ とおくと下記のようになる。
\begin{array}{c|c}
\textsf{判別式} & \textsf{実数解の個数} \\
\hline
D > 0 & 1 \\
\hline
D = 0 & 1 \sim 3 \\
\hline
D < 0 & 3 \\
\end{array}
なお，$D < 0$ のとき，$(p / 3)^3 < -(q / 2)^2 \le 0$ となるから必ず $p < 0$ であるとも言える。
カルダノの解法
まずはカルダノの解法を紹介する。$y = u + v$ とおいて変数を二つに増やす。
(u + v)^3 + p(u + v) + q = 0
\tag{10}
u^3 + v^3 + q + (3uv + p)(u + v) = 0
\tag{11}
ここで二つに増やした変数 $u$，$v$ が下記の連立方程式を満たせば，三次方程式の解となっている。
\left\{ \begin{align}
u^3 + v^3 + q = 0 \\
3uv + p = 0
\end{align} \right.
\tag{12}
これを $u^3$ と $v^3$ の連立方程式に変形する。
\left\{ \begin{align}
u^3 + v^3 &= -q \\
u^3 \cdot v^3 &= -\frac{p^3}{27}
\end{align} \right.
\tag{13}
解と係数の関係より，$u^3$ と $v^3$ は下記の二次方程式の解となっている。
t^2 + qt - \frac{p^3}{27} = 0
\tag{14}
これを解けば，$u^3$ と $v^3$ の値が得られる。これらの立方根より $u$ と $v$ が得られ，元の三次方程式の解 $y$ は $u$ と $v$ の和になるという流れである。
なお，解の個数に注意する必要があることを予め断っておく。式$(14)$を解くと（複素数解も許容すると）二つの解を得るのでそれぞれ $u^3$ と $v^3$ に割り当てれば良い。問題は $u$ と $v$ である。 三乗して $u^3$ と $v^3$ になる値はそれぞれ３つずつあるのだ。元の三次方程式の解 $y$ は $u$ と $v$ の和であるから合計９個の組み合わせが存在し得るが，もちろん正しい解はこのうち３個のみである。
この選別が非常に面倒である。式$(14)$の判別式 $D = q^2 + 4p^3 / 27$ を考えると，$D > 0$ のときは $u^3$ ，$v^3$ ともに実数解が得られ，その立方根の $u$ と $v$ にも実数が必ず含まれるので，実数同士の和が元の三次方程式の唯一の実数解 $y$ となる。
ところが $D < 0$ のとき，$u^3$ と $v^3$ はともに複素数になるにもかかわらず，それらの立方根 $u$，$v$ の和の組み合わせ９通りの中に実数となるものが３組生じるのだ。なぜなら $D < 0$ とは元の三次方程式が相異なる三つの実数解を持つ条件だからである。
すなわち実数解を得ることができると分かっている場合でも，複素数を経由して計算を進めなくてはならないのだ。複素数型をネイティブにサポートしている言語でなければ面倒である。というわけでカルダノの解法をそのまま実装するのは止めた。
非カルダノ解法に挑戦
敢えてカルダノの解法を用いず，ちょっとオリジナルっぽい解法に挑戦してみよう。まずは作為的だが，変数 $y$ を次のように置き換える。
y = z - \frac{p}{3z}
\tag{15}
\left(z - \frac{p}{3z}\right)^3 + p\left(z - \frac{p}{3z}\right) + q = 0
\tag{16}
z^3 - \cancel{pz} + \cancel{\frac{p^2}{3z}} - \frac{p^3}{27z^3} + \cancel{pz} - \cancel{\frac{p^2}{3z}} + q = 0
\tag{17}
こうすると次数が増えて六次方程式になってしまったが，もう解ける形だ
z^6 + qz^3 - \frac{p^3}{27} = 0
\tag{18}
$t = z^3$ とおくと $t$ に関する二次方程式に落とし込むことができる。よく見るとカルダノ解法における式 $(14)$ と全く同じ形だ。
t^2 + qt - \frac{p^3}{27} = 0
\tag{19}
ここで $t$ が相異なる二つの実数解を持つという条件として判別式 $D$ を考える。
D = q^2 + \frac{4p^3}{27} > 0
\tag{20}
これは元の三次方程式において実数解が一つしかない条件と一致することに気づく。そう，これは元の三次方程式の唯一の実数解を得られる条件なのだ。
あ れ れ ？ お か し い ぞ ？
ここで英邁な皆さんは謎に気づいただろうか？$t$ に関する二次方程式$(19)$を解けば，相異なる実数解二つを得られる。その立方根を求めると，$z$ についても相異なる実数値二つを得る。これらを式$(15)$に代入すると $y$ についても相異なる実数値二つを得てしまうのではないだろうか？
だがしかし，そんなことは起こらない。確認してみよう。$t$ に関する二次方程式$(19)$を解いて得られた相異なる実数解二つを $t_1$，$t_2$ とおくと解と係数の関係から
t_1 \cdot t_2 = -\frac{p^3}{27}
\tag{21}
という関係になる。$t$ の立方根である $z$ についてもそれぞれ $z_1$，$z_2$ とおくと
z_1 \cdot z_2 = -\frac{p}{3}
\tag{22}
となる。これらを式$(15)$に代入して得られる $y$ についても $y_1$，$y_2$ とおくと
y_1 = z_1 - \frac{p}{3z_1} = -\frac{p}{3z_2} + z_2 = y_2
\tag{23}
となるので $y$ の値まで変換するとどちらも同じ値になり，別に実数解の数が増える訳ではないのだ。
なお，立方根の有理化を図ると
\begin{align}
y_1 &= y_2 = z_1 + z_2 = \sqrt[3]{t_1} + \sqrt[3]{t_2} \\
&= \sqrt[3]{-\frac{q}{2} + \sqrt{\left(\frac{q}{2}\right)^2 + \left(\frac{p}{3}\right)^3}}
+ \sqrt[3]{-\frac{q}{2} - \sqrt{\left(\frac{q}{2}\right)^2 + \left(\frac{p}{3}\right)^3}}
\end{align}
\tag{24}
となる。
三倍角公式による解法
さっき非カルダノ解法に挑戦するといったが，賢明な読者はカルダノの解法と実質変わらないことに気づいたかもしれない。結局のところ，複素数が登場しないよう元の三次方程式が相異なる三つの実数解を持つ条件（判別式 $D < 0$ のとき）を避けている。
三つの実数解を持つ条件でも複素数の登場を避けるための手段として，三角関数の三倍角公式を用いる方法が知られている。三倍角公式とは
\begin{align}
\cos{3\theta} &= \cos{(2\theta + \theta)} \\
&= \cos{2\theta}\cos\theta - \sin{2\theta}\sin\theta \\
&= (2\cos^2\theta - 1)\cos\theta - 2\sin^2\theta\cos\theta \\
&= 2\cos^3\theta - \cos\theta - 2(1 - \cos^2\theta)\cos\theta \\
&= 4\cos^3\theta - 3\cos\theta \tag{25}
\end{align}
という関係を用いる解法だ。本解法を用いる大前提として元の三次方程式が相異なる三つの実数解を持つ（判別式 $D < 0$）とする。
これも作為的だが
y = 2\sqrt{-\frac{p}{3}}\cos\theta
\tag{26}
とする。なお判別式 $D < 0$ のときは必ず $p < 0$ なので，根号の中身は必ず正の値になっているので安心されたい。で，これを式$(3)$ に代入する。
8\left(-\frac{p}{3}\right)\sqrt{-\frac{p}{3}}\cos^3\theta + 2p\sqrt{-\frac{p}{3}}\cos\theta + q = 0
4\left(-\frac{p}{3}\right)\sqrt{-\frac{p}{3}}\cos^3\theta + p\sqrt{-\frac{p}{3}}\cos\theta + \frac{q}{2} = 0
4\left(-\frac{p}{3}\right)\sqrt{-\frac{p}{3}}\cos^3\theta - 3\left(-\frac{p}{3}\right)\sqrt{-\frac{p}{3}}\cos\theta + \frac{q}{2} = 0
\sqrt{\left(-\frac{p}{3}\right)^3}\underbrace{\left(4\cos^3\theta - 3\cos\theta\right)}_{\cos{3\theta}} + \frac{q}{2} = 0
\def\dfrac{\displaystyle\frac}
\cos{3\theta}= \frac{-\dfrac{q}{2}}{\sqrt{\left(-\dfrac{p}{3}\right)^3}}
以上より，以下の逆三角関数を解けば $\theta$ を求めることができる。
\theta = \frac{1}{3} \arccos \frac{-\dfrac{q}{2}}{\sqrt{\left(-\dfrac{p}{3}\right)^3}}
\tag{27}
逆正接関数を用いて書き直すと下記のようになる。なお，引数の順番は①Ｙ軸座標，②Ｘ軸座標とする。
\theta = \frac{1}{3} \arctan2\left(\sqrt{-\left\{\left(\frac{q}{2}\right)^2 + \left(\frac{p}{3}\right)^3\right\}}, -\frac{q}{2}\right)
\tag{28}
直角三角形を描いてみると分かり易い。
ここでＹ軸座標は常にゼロ以上であることから $0 \le 3\theta \le \pi$，すなわち $0 \le \theta \le \pi / 3$ であるから $\theta$ に加えて $\theta \pm 2\pi /3$ も式$(27)$を満たすので，これらを式$(26)$に代入すれば三つの実数解を得られる。
実装コード
JavaScript コードを以下に示す。なお，WSH の JavaScript エンジンは立方根を求めるメソッド Math.cbrt() をサポートしていないので代替品を用意した。べき乗を求めるメソッド Math.pow() を使えば代替できるが，引数に負の値を取れないことに注意する。
//------------------------------------------------------------------------------
// 実数の立方根を求める
//------------------------------------------------------------------------------
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;
}
テスト用コード
テスト用コードを示す。先ほどの 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 < 4) {
WScript.StdErr.WriteLine("実数係数の三次方程式の実数解を求めます。");
WScript.StdErr.WriteLine("");
WScript.StdErr.WriteLine("TEST(.WSF) [A] [B] [C] [D]");
WScript.StdErr.WriteLine("");
WScript.StdErr.WriteLine("三次方程式の係数を以下のように定義します。");
WScript.StdErr.WriteLine("A･x^3 + B･x^2 + C･x + D = 0");
return -1;
}
var k = [];
for(var i = 0; i < 4; i++)
k.push(parseFloat(args(i)));
if(k[0] == 0) {
WScript.StdErr.WriteLine("三次の係数は非ゼロにして下さい!!");
return -1;
}
var ret = cubic_equation(k[0], k[1], k[2], k[3]);
for(var i = 0; i < ret.length; i++)
WScript.Echo(ret[i]);
return 0;
}
</script>
</job>
</package>
テスト結果
引数なし，または引数が４個に満たない場合はヘルプメッセージを表示する。
C:\Qiita>test
実数係数の三次方程式の実数解を求めます。
TEST(.WSF) [A] [B] [C] [D]
三次方程式の係数を以下のように定義します。
A･x^3 + B･x^2 + C･x + D = 0
三次の係数にゼロを与えるとエラーになる。
C:\Qiita>test 0 0 0 0
三次の係数は非ゼロにして下さい!!
実数解が１個の場合を試す。$x^3 + 1 = (x + 1)(x^2 - x + 1) = 0$，すなわち $x = -1$ となる。
C:\Qiita>test 1 0 0 1
-1
実数解が３個の場合を試す。$x^3 - x = x(x - 1)(x + 1) = 0$，すなわち $x = 0, \pm 1$ となる。
C:\Qiita>test 1 0 -1 0
1
7.07050159149938E-17
-1
三重解の場合を試す。$x^3 = 0$，すなわち $x = 0$ となる。
C:\Qiita>test 1 0 0 0
-0
0
-0
次回予告
残念ながら四次方程式の解法（フェラーリの解法）まで到達できなかった。あと，Excel VBA 版も作りたい。
参考文献