3次方程式(\ref{eq.problem})の3つの解(複素数)を並べた配列を関数内で作成し、
その先頭アドレスを返す。
The first address of an array created in this function
which consists of the three solutions (complex numbers)
of eq. (\ref{eq.problem}).
◆使用例(Example)
struct im ∗solution=eqsol3(1.2,3.4,5.6,7.8);
◆計算式とアルゴリズム
(Formula and algorithm)
1. \(a=0\)の判定と処理
(Judgement and processings for \(a=0\))
\(a=0\)のときは2次方程式として解くことも可能であるが、
気づかず意図しない結果を生じるのを防ぐために
\(a=0\)の場合はエラー終了する仕様にしている。
Although it is possible to solve a quadratic equation
in case of \(a=0\),
this function treats \(a=0\) as an error
to avoid obtaining unexpected results without awaring it.
\(a=0\)であるか否かの判定には関数doublecmpを使用し、
\(a\)と\(\max\{|b|,|c|,|d|\}\)の比較によって判定している。
Function doublecmp is used for the judgement of \(a=0\),
where \(a\) is compared with \(\max\{|b|,|c|,|d|\}\).
2. \(a\neq 0\)のときの解
(Solution in case of \(a \neq 0\))
3次方程式(\ref{eq.problem})の解\(x^0\), \(x^1\), \(x^2\)は
以下のように与えられる。
\[\begin{equation}
x^m= \sqrt[3]{r_1}e^{i(2\theta_2-\theta_1+2m\pi)/3}
+\sqrt[3]{r_2}e^{-i(2\theta_1-\theta_2+2m\pi)/3}
-\frac{b}{3a} \hspace{1em} (m=0,1,2)
\label{eq.solution}
\end{equation}\]
ここで\(r_1\), \(\theta_1\), \(r_2\), \(\theta_2\)は
以下の式によって定義される。
\[\begin{equation}
A \equiv \frac{3ac-b^2}{3a^2}
\label{eq.A}
\end{equation}\]
\[\begin{equation}
B \equiv \frac{2b^3-9abc+27a^2d}{27a^3}
\label{eq.B}
\end{equation}\]
\[\begin{equation}
C \equiv \frac{B^2}{4}+\frac{A^3}{27}
\equiv R e^{i\phi}
\hspace{1em} (R \geq 0; -\pi < \phi \leq \pi)
\label{eq.C}
\end{equation}\]
\[\begin{equation}
D_1 \equiv -\frac{B}{2}+\sqrt{R}e^{i\phi/2}
\equiv r_1 e^{i\theta_1}
\hspace{1em} (r_1 \geq 0; -\pi < \theta_1 \leq \pi)
\label{eq.D1}
\end{equation}\]
\[\begin{equation}
D_2 \equiv -\frac{B}{2}-\sqrt{R}e^{i\phi/2}
\equiv r_2 e^{-i\theta_2}
\hspace{1em} (r_2 \geq 0; -\pi < \theta_2 \leq \pi)
\label{eq.D2}
\end{equation}\]
The three solutions \(x^0\), \(x^1\), and \(x^2\) for
the cubic equation of (\ref{eq.problem})
are given by eq. (\ref{eq.solution}),
where \(r_1\), \(\theta_1\), \(r_2\), and \(\theta_2\) are defined by
eqs (\ref{eq.A})-(\ref{eq.D2}).
解の導出
(Derivation of the solutions)
解(\ref{eq.solution})の導出過程は以下の通り。
The solution (\ref{eq.solution}) is obtained as follows.
方程式(\ref{eq.problem})の両辺を\(a\)で割って
\[\begin{equation}
x^3+\frac{b}{a}x^2+\frac{c}{a}x+\frac{d}{a}=0
\label{eq.problem.normalized}
\end{equation}\]
この式の左辺は
\[\begin{equation}
\left(x+\frac{b}{3a}\right)^3
+\frac{3ac-b^2}{3a^2}\left(x+\frac{b}{3a}\right)
+\frac{2b^3-9abc+27a^2d}{27a^3}
\label{eq.problem.left_side.arranged}
\end{equation}\]
と変形できるので、定数\(A\), \(B\)を(\ref{eq.A})(\ref{eq.B})のようにおき、
\[\begin{equation}
X \equiv x+\frac{b}{3a}
\label{eq.Xdef}
\end{equation}\]
と変数変換すれば
\[\begin{equation}
X^3+AX+B=0
\label{eq.X}
\end{equation}\]
を得る。これを解くために更に
\[\begin{equation}
X=Y-\frac{A}{3Y}
\label{eq.Ydef}
\end{equation}\]
と変数変換すると
\[\begin{equation}
X^3=Y^3-\frac{A^3}{27Y^3}-A\left(Y-\frac{A}{3Y}\right)
=Y^3-\frac{A^3}{27Y^3}-AX
\label{eq.Y.arranged}
\end{equation}\]
であるので、これらを(\ref{eq.X})式に代入して整理すれば
\[\begin{equation}
Y^6+BY^3-\frac{A^3}{27}=0
\label{eq.Y}
\end{equation}\]
が得られる。
Dividing the both sides of eq. (\ref{eq.problem}) by \(a\) results in
eq. (\ref{eq.problem.normalized}),
and the left hand side of this equation
is arranged as (\ref{eq.problem.left_side.arranged}).
Defining constants \(A\) and \(B\) as eqs (\ref{eq.A})(\ref{eq.B})
and converting the variable from \(x\) to \(X\) as (\ref{eq.Xdef})
yield (\ref{eq.X}).
A further conversion of the variable from \(X\) to \(Y\) as (\ref{eq.Ydef})
gives eq. (\ref{eq.Y.arranged}).
Inserting eq. (\ref{eq.X}) into (\ref{eq.Y.arranged}),
one finally obtains (\ref{eq.Y}).
(\ref{eq.Y})式は\(Y^3\)についての2次方程式になっており、その解は
\[\begin{equation}
Y^3
=\frac{-B\pm\sqrt{B^2+\frac{4A^3}{27}}}{2}
=-\frac{B}{2}\pm\sqrt{\frac{B^2}{4}+\frac{A^3}{27}}
\label{eq.Y3.original}
\end{equation}\]
と書ける。
\[\begin{equation}
C \equiv \frac{B^2}{4}+\frac{A^3}{27}
\label{eq.C.def}
\end{equation}\]
とおき、この\(C\)を
\[\begin{equation}
C=R e^{i\phi}
\hspace{1em} (R \geq 0; -\pi < \phi \leq \pi)
\label{eq.R.phi.def}
\end{equation}\]
と極座標表示すれば(\ref{eq.Y3.original})は
\[\begin{equation}
Y^3=-\frac{B}{2}\pm\sqrt{R}e^{i\phi/2}
\label{eq.Y3}
\end{equation}\]
と書き直せる。\(Y^3\)に関するこれら2つの解を
\[\begin{equation}
D_1 \equiv -\frac{B}{2}+\sqrt{R}e^{i\phi/2}
\label{eq.D1.def}
\end{equation}\]
\[\begin{equation}
D_2 \equiv -\frac{B}{2}-\sqrt{R}e^{i\phi/2}
\label{eq.D2.def}
\end{equation}\]
とおき、これらを
\[\begin{equation}
D_1=r_1 e^{i\theta_1}
\hspace{1em} (r_1 \geq 0; -\pi < \theta_1 \leq \pi)
\label{eq.r1.theta1.def}
\end{equation}\]
\[\begin{equation}
D_2=r_2 e^{-i\theta_2}
\hspace{1em} (r_2 \geq 0; -\pi < \theta_2 \leq \pi)
\label{eq.r2.theta2.def}
\end{equation}\]
と極座標表示すると、\(Y\)に関する解として
\[\begin{equation}
Y=Y_1^m, Y_2^m \hspace{1em} (m=0,1,2)
\label{eq.Y.solution}
\end{equation}\]
\[\begin{equation}
Y_1^m\equiv \sqrt[3]{r_1}e^{i(\theta_1+2m\pi)/3}
\label{eq.Y1}
\end{equation}\]
\[\begin{equation}
Y_2^m\equiv \sqrt[3]{r_2}e^{-i(\theta_2+2m\pi)/3}
\label{eq.Y2}
\end{equation}\]
が得られる。これら\(Y\)に関する解を(\ref{eq.Ydef})に代入すると
\(X\)についての解として
\[\begin{equation}
X=X_1^m, X_2^m
\label{eq.X.solution}
\end{equation}\]
\[\begin{equation}
X_1^m\equiv Y_1^m-\frac{A}{3Y_1^m}
=\sqrt[3]{r_1}e^{i(\theta_1+2m\pi)/3}
-\frac{A}{3\sqrt[3]{r_1}}e^{-i(\theta_1+2m\pi)/3}
\label{eq.X1.original}
\end{equation}\]
\[\begin{equation}
X_2^m\equiv Y_2^m-\frac{A}{3Y_2^m}
=\sqrt[3]{r_2}e^{-i(\theta_2+2m\pi)/3}
-\frac{A}{3\sqrt[3]{r_2}}e^{i(\theta_2+2m\pi)/3}
\label{eq.X2.original}
\end{equation}\]
が得られる。
Eq. (\ref{eq.Y}) is a quadratic equation for \(Y^3\),
the solution of which is given by eq. (\ref{eq.Y3.original}).
We now define \(C\) as (\ref{eq.C.def})
and represent it in a polar coordinate format as (\ref{eq.R.phi.def}).
Then eq. (\ref{eq.Y3.original}) reduces to (\ref{eq.Y3}).
Let \(D_1\) and \(D_2\) be these two solutions,
defined by eqs (\ref{eq.D1.def}) and (\ref{eq.D2.def}),
and represent them in polar coordinate formats as
eqs (\ref{eq.r1.theta1.def}) and (\ref{eq.r2.theta2.def}).
Then the solution for \(Y\) is given by eq. (\ref{eq.Y.solution}),
where \(Y_1^m\) and \(Y_2^m\) are given by
eqs (\ref{eq.Y1}) and (\ref{eq.Y2}), respectively.
Inserting these solutions for \(Y\) into (\ref{eq.Ydef}) gives
the solutions for \(X\) represented by
eqs (\ref{eq.X.solution})-(\ref{eq.X2.original}).
ところで(\ref{eq.D1.def})-(\ref{eq.r2.theta2.def})より
\[\begin{eqnarray}
r_1r_2 e^{i(\theta_1-\theta_2)}
&=& D_1D_2
\nonumber \\
&=& \left(-\frac{B}{2}+\sqrt{R}e^{i\phi/2}\right)
\left(-\frac{B}{2}-\sqrt{R}e^{i\phi/2}\right)
\nonumber \\
&=& \frac{B^2}{4}-Re^{i\phi}
\nonumber \\
&=& \frac{B^2}{4}-C
\nonumber \\
&=& -\frac{A^3}{27}
\label{eq.D1D2}
\end{eqnarray}\]
が成り立つ。\(\theta_1\), \(\theta_2\)の定義範囲より
\[\begin{equation}
-2\pi < \theta_1-\theta_2 < 2\pi
\label{eq.theta_diff}
\end{equation}\]
である。これらの性質を用いると(\ref{eq.X1.original})(\ref{eq.X2.original})を
より見通しの良い形に変形でき、両者を比較することができる。
以下では\(A\)の符号に応じて場合分けしてこの比較を行う。
From (\ref{eq.D1.def})-(\ref{eq.r2.theta2.def}),
one obtains eq. (\ref{eq.D1D2}),
and the definition ranges of \(\theta_1\) and \(\theta_2\) mean that
the range of \(\theta_1-\theta_2\) is limited as (\ref{eq.theta_diff}).
Using these relations, eqs (\ref{eq.X1.original}) and (\ref{eq.X2.original})
can be arranged to more comprehensive forms and can be easily compared.
These forms depend on the sign of \(A\).
まず\(A\leq 0\)のとき、
(\ref{eq.D1D2})より\(e^{i(\theta_1-\theta_2)}=1\)となり、
(\ref{eq.theta_diff})の条件と合わせて\(\theta_1-\theta_2=0\)を得る。
またこのとき(\ref{eq.D1D2})より
\(\sqrt[3]{r_1}\sqrt[3]{r_2}=-\frac{A}{3}\)であるので、
これらの関係を用いて(\ref{eq.X1.original})(\ref{eq.X2.original})式は
\[\begin{equation}
X_1^m
=\sqrt[3]{r_1}e^{i(\theta_1+2m\pi)/3}
+\sqrt[3]{r_2}e^{-i(\theta_2+2m\pi)/3}
\label{eq.X1.negativeA}
\end{equation}\]
\[\begin{equation}
X_2^m
=\sqrt[3]{r_2}e^{-i(\theta_2+2m\pi)/3}
+\sqrt[3]{r_1}e^{i(\theta_2+2m\pi)/3}
\label{eq.X2.negativeA}
\end{equation}\]
と書き直せる。
これより\(X_1^m=X_2^m\)であり、
独立な解としては\(X_1^m\)のみを考えれば良いことが分かる。
In case of \(A\leq 0\),
(\ref{eq.D1D2}) gives \(e^{i(\theta_1-\theta_2)}=1\),
from which one obtains \(\theta_1-\theta_2=0\)
taking into account eq. (\ref{eq.theta_diff}).
In this case eq. (\ref{eq.D1D2}) reduces to
\(\sqrt[3]{r_1}\sqrt[3]{r_2}=-\frac{A}{3}\).
Using these relations, eqs (\ref{eq.X1.original}) and (\ref{eq.X2.original})
can be rewritten as (\ref{eq.X1.negativeA}) and (\ref{eq.X2.negativeA}),
respectively.
Therefore \(X_1^m=X_2^m\) holds, so that
only \(X_1^m\) needs to be considered as the independent solutions.
次に\(A>0\)のとき、(\ref{eq.D1D2})より
\(e^{i(\theta_1-\theta_2)}=-1\)となり、
(\ref{eq.theta_diff})の条件と合わせて\(\theta_1-\theta_2=\pm\pi\)を得る。
ところで\(A>0\)のときは
(\ref{eq.C.def})より\(C>\frac{B^2}{4}>0\)であるので
(\ref{eq.R.phi.def})において\(R>\frac{B^2}{4}\), \(\phi=0\)となる。
したがって(\ref{eq.D1.def})(\ref{eq.D2.def})式から
\(D_1=-\frac{B}{2}+\sqrt{R}>0\), \(D_2=-\frac{B}{2}-\sqrt{R}<0\)となり、
\(\theta_1=0\), \(\theta_2=\pi\), \(\theta_1-\theta_2=-\pi\)であることが分かる。
またこのとき(\ref{eq.D1D2})より
\(\sqrt[3]{r_1}\sqrt[3]{r_2}=\frac{A}{3}\)となるので
(\ref{eq.X1.original})(\ref{eq.X2.original})式は
\[\begin{eqnarray}
X_1^m
&=& \sqrt[3]{r_1}e^{i(\theta_1+2m\pi)/3}
-\sqrt[3]{r_2}e^{-i(\theta_2-\pi+2m\pi)/3}
\nonumber \\
&=& \sqrt[3]{r_1}e^{i(\theta_1+2m\pi)/3}
+\sqrt[3]{r_2}e^{i\pi}e^{-i(\theta_2-\pi+2m\pi)/3}
\nonumber \\
&=& \sqrt[3]{r_1}e^{i(\theta_1+2m\pi)/3}
+\sqrt[3]{r_2}e^{-i[\theta_2+2(m-2)\pi]/3}
\label{eq.X1.positiveA}
\end{eqnarray}\]
\[\begin{eqnarray}
X_2^m
&=& \sqrt[3]{r_2}e^{-i(\theta_2+2m\pi)/3}
-\sqrt[3]{r_1}e^{i(\theta_1+\pi+2m\pi)/3}
\nonumber \\
&=& \sqrt[3]{r_2}e^{-i(\theta_2+2m\pi)/3}
+\sqrt[3]{r_1}e^{i\pi}e^{i(\theta_1+\pi+2m\pi)/3}
\nonumber \\
&=& \sqrt[3]{r_2}e^{-i(\theta_2+2m\pi)/3}
+\sqrt[3]{r_1}e^{i[\theta_1+2(m+2)\pi]/3}
\label{eq.X2.positiveA}
\end{eqnarray}\]
と書き直せる。
これより\(X_1^0=X_2^1\), \(X_1^1=X_2^2\), \(X_1^2=X_2^0\)となり、
独立な解としては\(X_1^m\)のみを考えれば良いことが分かる。
また、その解については(\ref{eq.X1.positiveA})を対称性の良い形に変形すると
\[\begin{equation}
X_1^m
=\sqrt[3]{r_1}e^{2\pi i/3+i[\theta_1+2(m-1)\pi]/3}
+\sqrt[3]{r_2}e^{2\pi i/3-i[\theta_2+2(m-1)\pi]/3}
\label{eq.X1.positiveA.arrange1}
\end{equation}\]
となるが、\(m\)の定義を変えて
\[\begin{equation}
X_1^m
=\sqrt[3]{r_1}e^{2\pi i/3+i(\theta_1+2m\pi)/3}
+\sqrt[3]{r_2}e^{2\pi i/3-i(\theta_2+2m\pi)/3}
\label{eq.X1.positiveA.arrange2}
\end{equation}\]
としても3つの解の順番が変わるだけで解自体は変わらない。
We next consider a case where \(A>0\).
In this case, eq. (\ref{eq.D1D2}) gives \(e^{i(\theta_1-\theta_2)}=-1\),
from which \(\theta_1-\theta_2=\pm\pi\) is obtained based on
eq. (\ref{eq.theta_diff}).
Indeed, when \(A>0\),
eq. (\ref{eq.C.def}) indicates \(C>\frac{B^2}{4}>0\),
so that \(R>\frac{B^2}{4}\)
and \(\phi=0\) hold in eq. (\ref{eq.R.phi.def}).
Then eqs (\ref{eq.D1.def}) and (\ref{eq.D2.def}) indicate
\(D_1=-\frac{B}{2}+\sqrt{R}>0\)
and \(D_2=-\frac{B}{2}-\sqrt{R}<0\), respectively,
from which \(\theta_1=0\), \(\theta_2=\pi\),
and thus \(\theta_1-\theta_2=-\pi\) are obtained.
Eq. (\ref{eq.D1D2}) also indicates
\(\sqrt[3]{r_1}\sqrt[3]{r_2}=\frac{A}{3}\),
so that eqs (\ref{eq.X1.original}) and (\ref{eq.X2.original}) are rewritten as
(\ref{eq.X1.positiveA}) and (\ref{eq.X2.positiveA}), respectively.
These equations indicate
\(X_1^0=X_2^1\), \(X_1^1=X_2^2\), and \(X_1^2=X_2^0\),
so that only \(X_1^m\) are needed as the independent solutions.
Eq. (\ref{eq.X1.positiveA}) can further be arranged to more symmetric form
of (\ref{eq.X1.positiveA.arrange1}),
and shifting the definition of \(m\) by 1,
which does not change the solutions but changes only the order of them,
gives a simpler form of eq. (\ref{eq.X1.positiveA.arrange2}).
\(A\leq 0\)の場合に\(\theta_1-\theta_2=0\)であったこと、
\(A>0\)の場合には\(\theta_1-\theta_2=-\pi\)であったことを踏まえると、
(\ref{eq.X1.negativeA})(\ref{eq.X1.positiveA.arrange2})はまとめて
\[\begin{eqnarray}
X_1^m
&=& \sqrt[3]{r_1}e^{2i(\theta_2-\theta_1)/3+i(\theta_1+2m\pi)/3}
+\sqrt[3]{r_2}e^{2i(\theta_2-\theta_1)/3-i(\theta_2+2m\pi)/3}
\nonumber \\
&=& \sqrt[3]{r_1}e^{i(2\theta_2-\theta_1+2m\pi)/3}
+\sqrt[3]{r_2}e^{-i(2\theta_1-\theta_2+2m\pi)/3}
\label{eq.X1}
\end{eqnarray}\]
と書ける。これが\(X\)に関する解であり、(\ref{eq.Xdef})に代入すれば
\(x\)に関する解として(\ref{eq.solution})が得られる。
Since \(\theta_1-\theta_2=0\) and \(\theta_1-\theta_2=-\pi\) hold
in cases of \(A\leq 0\) and \(A>0\), respectively,
eqs (\ref{eq.X1.negativeA}) and (\ref{eq.X1.positiveA.arrange2})
can be unified as (\ref{eq.X1}).
This is the solution for \(X\),
and inserting it into (\ref{eq.Xdef}) results in
the solution for \(x\), given by eq. (\ref{eq.solution}).