カルダノの公式
カルダノさんによる三次方程式の解の公式のことをカルダノの公式と呼びます。
まずはこれを導出してみます。以下、の形の実数係数の三次方程式の解を考えていきます。
1. 3次の係数を1にする
まずは3次の項の係数を1にするように式を変形します。この操作によって解は変化しません。
3次の係数aで全体を割ってやります。残った2次以下の係数をそれぞれA,B,Cと置き、次に進みます。
2. 2次の項を削除する
次は、2次の項を削除します。式から消してしまうのです・・・。
そのために、と置換してやります。
代入し、展開し、項を集めていきます、めんどい。しかし見事に消えてしまいました。残ったyの係数と切片をそれぞれp,qと置き、次に行きましょう。
ちなみに、この操作をn次方程式に一般化したものをチルンハウス変換と呼びます。2次なら平方完成、この場合は立法完成でしょうか。n次のチルンハウス変換はn-1次の項の係数を0にします。どうやら方程式の変曲点(微分係数が減少から増大に転じる点、もしくはその逆)をy=0の点に平行移動させているようです。
3. 式の次数を落とす(2次式にして考える)
3次の係数を1にし、2次の項を消したとしても、まだ式は3次式のままです。ですので、次はこの式の次数を落としてやります。
と置いて式を変形していきます。
ここで、が成り立つ時を考えると
という連立方程式が得られるので、この式を解けばu(v)が求められそうです。
まず二つ目の式からvを求め
これを一つ目の式に代入してやります
この式はの2次式になっている事が分かるので、2次方程式の解の公式を使えばを求められます。
また、全く同じようにを求めることもでき、ルートの前の符号をそれぞれに割り当ててを求めることができました。
4. 得られた各値より解を求める
先ほど得られたの立方根を取れば、が求まります。が得られれば、よりが得られ、よりがめでたく求まる筈です。
しかし、実数/複素数の立方根は3つあるのでそれを考慮しなければなりません。とはいえ、1の3つの立方根()を使えば解決できます。一般にあるの立方根はの3つある事から、yは以下のようになります。
ここで、です。
とyが3つ求まりました。これを使えばxも求まります。
こうしてようやくxにたどり着くことが出来ました。
しかし、解の公式というほどピシッと一つの式で出るわけではなく、何とか代数的に求められるよ!という事くらいしか分からない式ですね・・・
3次方程式の判別式
得られたカルダノの公式をC++コードにコピペする前に、解のパターンを考えてみましょう。2次方程式の解が重解だったり虚数解だったりしたように、3次方程式でも同じようになりそうです。そして、式中にが出て着ていたり、ルートがあることからもそれは間違いなさそうです。
2次方程式ではそれを判別式によって場合分けしていましたので、3次方程式でも判別式を使って解を見てみましょう。
2次方程式の判別式は解の公式中にあるルートの中身です。なぜなら、ルートの中身の符号によって虚数が出てくるかが決まるうえ、ルートの中が0ならば±で変化する項が消えることから重解となる事が分かります。
先ほど求めた3次方程式の解の公式中のの式を見てみると、ここにいかにも分岐しそうな±付きのルートがあり、この中身を判別式として利用できそうです。
p,qは元の方程式の各係数から計算される値なので、任意の値が入るものとしましょう。すると、やはりこの式には3つのパターンがある事が分かります。つまり、、、のいずれかです。この時のu,vの様子を考えてみましょう。
D = 0 となる場合
判別式D=0となる場合のの式はかなり簡単になります。
とuを一つ決めます、qが実数であることからこの立方根も実数なものを選択できます。カルダノの公式導出中(3step目)に出てきた連立式の一部()から対応するvを求めます。
ここで、判別式の中()に目を向けます。
判別式が0になるという事は、これらの値の絶対値は等しいはずです。そして、qが実数であることからの符号は必ず正になるので、が負になっており、pの符号が負である事が分かります。
それらのことから
となるので、先ほどの式に入れてやれば
と、u=vとなる事が分かりました。ここからyを求めるのですが、の関わってくる部分を先に見ておきましょう。
となり、y3も同様になる事も分かるでしょう。結果各yの値は
となります。
qが実数であることからこれらはすべて実数であり、y2とy3は明らかに重解です。また、q=0のときすべての値が等しくなり、3重解となる事が分かります。
また、このyからxを求めてもこれらの関係は変化しないことは明らかでしょう(より)。
これらのことから、判別式の時の三次方程式は重解となり、解は全て実数となる事が分かりました。
D > 0 となる場合
判別式が正となる場合、ルートを開いた結果は実数となり、ともに実数となります。
先ほどと同様、uを一つ決め、対応するvを求めてやります。
のように、がともに求まりました。先ほどと同じようにが出てくるところを先に見ておきます。
y2の導出ではy1の結果を使って式を省略しています(uとvの位置が入れ替わるだけ)。結果、各yの値は
となります。、であり、はともに実数なのでu,vも実数となることからy1は実数で、y2とy3は複素数になる事が分かります。そして、y2とy3は明らかに共役複素数の関係にあります。また、xを求めるにあたってもこの関係性は変化しません。
これらの事から、判別式の時の3次方程式の解は実数解一つと、互いに共役な複素数解2つとなる事が分かりました。
D < 0 となる場合
判別式が負となる場合、ルートを開くと複素数が出てくるため、も複素数となります。
そのまま解の公式に当てはめて計算していく事も出来ますが、負の数の平方根と複素数の立方根を考えなければならず、複雑になってしまいます。そこで、複素数の極形式に変形して考えてみます。
まずは、ノルムである絶対値。これは実部と虚部の値の二乗和平方根(つまりはユークリッド距離)となり、共役な複素数との積のルートでもあります。
となります。この時、Dの値は負であり、qは実数である事からの項の符号は必ず正です。なので、の項が負(pが負)になる事が分かります。
と、絶対値がきれいな形に求まりました。この絶対値と偏角(arctanで求まりますが、ここでは単にとしておきます)を使えば、複素数の立方根の一つは以下のように表すことが出来ます。
まだこれは複素数となります。なぜ偏角がとなるのかというと、極形式で複素数の積を考えると、偏角は和を取る形になるため、偏角の3乗根のうち一つはになります(しかし、残り2つの立方根はそこから±120°の所にあります)。
そして例のごとく、このuに対応するvを求めます。
と得られ、u,vの関係は共役となる事が分かりました。そして、例のごとくが出てくるyを求めます。
となるので、を求めますと
となるので、を求めると
これらからy2を求めます。なぜ約分したり進めないのかはすぐにわかります。
加法定理と三角関数の周期性を利用して式変形を行い、比較的すっきりとした形になりました。y3も同じように求めましょう。
同じような形で得られ、結果各yの値は
となります。途中で複素数が出て着ましたが、pが実数であるので最終的には全て実数となる事が分かるでしょう。そして、明らかに全て異なる値となります。また、当然ですがxを求めるにあたってもこれらの関係は変化しません。
これらの事から、判別式の時の3次方程式の解は異なる3つの実数解となる事が分かりました。
ちなみに、この場合の三次方程式の事を還元不能と言います。複素数を考えないように解こうとすると別の三次式が出現してしまうためです。還元不能な三次式を解くためにビエタの方法というcosの3倍角の公式を利用した別の方法もあります。
まとめ
-
- 重解(すべて実数解)
-
- 実数解一つと、2つの共役な複素数解
-
- 異なる3つの実数解
C++実装
さて、無事に各式が得られたのでC++コードにコピペしていきます。公式通り順番に進めていきましょう。
なお、floatでもdoubleでも同じコードになるのでテンプレートでジェネリックに書いていきます。そして、構造化束縛というC++17の機能を用いているのでC++17未満のコンパイラでは書き換えが必要です、VS2015とか。その場合はstd::tieを使うといいでしょう。
template<typename T> auto SolveCubicEquation(const T a, const T b, const T c, const T d){ return SolveCubicEquation(b/a, c/a, d/a); }
まずは第一段階、カルダノの公式step1。3次の係数を1にし、次に行きます。何のことはなくただの割り算です。
template<typename T> auto SolveCubicEquation(const T A, const T B, const T C){ //A/3 auto A_div3 = A / T(3.0); const T p = B - A_div3*A; //A^3/27 auto q_t = A_div3 * A_div3 * A_div3; //2A^3/27 q_t += q_t; const T q = q_t - A_div3*B + C; auto [x1, x2, x3] = SolveCubicEquation(p, q); //xn = yn - A/3 return std::make_tuple(x1 - A_div3, x2 - A_div3, x3 - A_div3); }
第2step、チルンハウス変換により2次の項を削除します。そして次のステップに渡し、帰ってきたからを引いて最終的な解とします。
template<typename T> auto SolveCubicEquation(const T p, const T q){ //p/3 const auto p_d3 = p / T(3.0); //q/2 const auto q_half = q / T(2.0); //D = (q/2)^2 + (p/3)^3 const auto D = q_half*q_half + p_d3*p_d3*p_d3; T x1{}; std::complex<T> x2{}, x3{}; if(std::fabs(D) < 1.0E-8) { // D == 0 重解 auto tmp = std::cbrt(q_half); x1 = T(-2.0) * tmp; x2 = x3 = tmp; } else if(D < 0.0) { // D < 0 異なる3つの実数解 //2√(-p/3) const auto sqrt_p_d3 = T(2.0) * std::sqrt(-p_d3); //θ/3 const auto arg = std::arg(std::complex<T>{ -q_half, std::sqrt(-D) }) / T(3.0); //2π/3 constexpr T pi2d3 = T(2.0 * 3.1415926535897932384626433832795 / 3.0); x1 = sqrt_p_d3*std::cos(arg); x2 = sqrt_p_d3*std::cos(arg + pi2d3); x3 = sqrt_p_d3*std::cos(arg + pi2d3 + pi2d3); } else { // 0 < D 1つの実数解と二つの共役な複素数解 //√D const auto D_root = std::sqrt(D); //u = v = ∛(-q/2 ± √D) const auto u = std::cbrt(-q_half + D_root); const auto v = std::cbrt(-q_half - D_root); x1 = u + v; //-(u + v)/2 const auto real = T(-0.5)*x1; //√3 constexpr T root3 = T(1.7320508075688772935274463415059); //√3(u - v)/2 const auto imag = root3 * (u - v) * T(0.5); x2 = {real, imag}; x3 = {real, -imag}; } return std::make_tuple(x1, x2, x3); }
キモとなる判別式による分岐とそれに応じた求解部分。基本的には公式をコピペしただけです。
異なる3つの実数解を求めるところでは、あえてとを分けて計算することで割り算の回数を1回に減らしています(constexprがコンパイル時評価される場合)。
template<typename T> auto SolveQuadraticEquation(const T a, const T b, const T c){ const auto D = b*b - T(4.0)*a*c; std::complex<T> x1{}, x2{}; if(std::fabs(D) < 1.0E-6) { x1 = x2 = -b / (a + a); } else if(D < 0.0) { const auto a2 = a + a; const auto real = -b / a2; const auto imag = std::sqrt(D) / a2; x1 = {real, imag}; x2 = {real, -imag}; } else { //-b - √D const auto tmp = -b - std::sqrt(D); //2c/(-b - √D) x1 = (c + c) / tmp; //(-b - √D)/2a x2 = tmp / (a + a); } return std::make_tuple(T(0.0), x1, x2); } template<typename T> auto SolveCubicEquation(const T a, const T b, const T c, const T d){ if(a != 0.0) { return SolveCubicEquation(b/a, c/a, d/a); } else if(b != 0.0) { //bx^2 + cx + d = 0 //x = (-b ± √(b^2 - 4ac))/2a return SolveQuadraticEquation(b, c, d); } else { constexpr std::complex<T> z = T(0.0); //cx + d = 0 //x = -d / c return std::make_tuple(-d/c, z, z); } }
おまけで、最初の3次式の係数が0である場合に次数を落とした式として解くように修正してやります。本当はc=0も書くべきですので、必要ならば追加してください。
あとは数値計算用に桁落ちとかの配慮が必要な気がしますがそれはまたの機会に・・・
コード全文とテスト実行
[Wandbox]三へ( へ՞ਊ ՞)へ ハッハッ
解の正しさはこちらのサイトで
三次方程式の解 - 高精度計算サイト
なんだかめっちゃ長くなりました、何か役に立っていれば幸いです・・・。