　
三次方程式・四次方程式を直接解きたいような場面はしばしば生じうると思います。例えばシミュレーション中に何度も四次方程式を解く必要があるときなどには、ニュートン法などの反復による数値解法では時間がかかりすぎてしまいます。
そこでC言語でこれらの解を求める関数を作りました。数式の部分はほとんど
https://onihusube.hatenablog.com/entry/2018/10/08/140426
のページを参考にさせてもらっています。またC++のコードも非常に参考になります(私はC++はよくわからないのでCで実装することにしました)。

#0 サンプルコード
とりあえず使用例を示します。これだけ見ればだいたい何をしているのか、どのように使うのかは分かると思います。

equation.c

#include <stdio.h>
#include <math.h>
#include <complex.h>

void SolveCubicEquation(double complex x[3],double a,double b,double c,double d);
void SolveQuarticEquation(double complex x[4],double a,double b,double c,double d,double e);
double Cuberoot(double x);
double complex Squreroot(double complex x);

int main(void){
	double complex x[3],y[4];
	SolveCubicEquation(x,1.0,2.0,3.0,4.0);//x^3+2x^2+3x+4=0を解く,結果はxに格納
	printf("==Cubic Equation==\n");
	for(int i=0;i<3;i++){
		printf("%lf + %lf i\n",creal(x[i]),cimag(x[i]));//実部と虚部を出力
	}
	printf("\n");
	printf("==Quartic Equation==\n");
	SolveQuarticEquation(y,1.0,2.0,3.0,4.0,5.0);//x^4+2x^3+3x^2+4x+5=0を解く,結果はyに格納
	for(int i=0;i<4;i++){
		printf("%lf + %lf i\n",creal(y[i]),cimag(y[i]));
	}
	return 0;
}

void SolveCubicEquation(double complex x[3],double a,double b,double c,double d){
	if(a == 0.0){
		printf("Error:a = 0.0\n");
		printf("This equation is NOT Cubic.\n");
	}
	double A,B,C,p,q,D;
	A = b/a;
	B = c/a;
	C = d/a;
	p = B-A*A/3.0;
	q = 2.0*A*A*A/27.0-A*B/3.0+C;
	D = q*q/4.0+p*p*p/27.0;
	if(D < 0.0){//three real solutions
		double theta = atan2(sqrt(-D),-q*0.5);
		x[0] = 2.0*sqrt(-p/3.0)*cos(theta/3.0)-A/3.0;
		x[1] = 2.0*sqrt(-p/3.0)*cos((theta+2.0*M_PI)/3.0)-A/3.0;
		x[2] = 2.0*sqrt(-p/3.0)*cos((theta+4.0*M_PI)/3.0)-A/3.0;
	}else{//single real solution and two imaginary solutions(c.c)
		double u = Cuberoot(-q*0.5+sqrt(D)),v = Cuberoot(-q*0.5-sqrt(D));
		x[0] = u+v-A/3.0;
		x[1] = -0.5*(u+v)+sqrt(3.0)*0.5*1i*(u-v)-A/3.0;
		x[2] = -0.5*(u+v)-sqrt(3.0)*0.5*1i*(u-v)-A/3.0;
	}
}

void SolveQuarticEquation(double complex x[4],double a,double b,double c,double d,double e){
	if(a == 0.0){
		printf("Error:a = 0.0\n");
		printf("This equation is NOT Quartic.\n");
	}
	double A,B,C,D,p,q,r,epsilon_temp = 0.000001;
	A = b/a;
	B = c/a;
	C = d/a;
	D = e/a;
	p = -6.0*pow(A/4.0,2.0)+B;
	q = 8.0*pow(A/4.0,3.0)-2.0*B*A/4.0+C;
	r = -3.0*pow(A/4.0,4.0)+B*pow(A/4.0,2.0)-C*A/4.0+D;
	double complex t_temp[3];
	SolveCubicEquation(t_temp,1.0,-p,-4.0*r,4.0*p*r-q*q);
	double t = creal(t_temp[0]);
	double complex m = Squreroot(t-p);
	x[0] = (-m+Squreroot(-t-p+2.0*q/m))*0.5-A/4.0;
	x[1] = (-m-Squreroot(-t-p+2.0*q/m))*0.5-A/4.0;
	x[2] = (m+Squreroot(-t-p-2.0*q/m))*0.5-A/4.0;
	x[3] = (m-Squreroot(-t-p-2.0*q/m))*0.5-A/4.0;
}



double Cuberoot(double x){
	if(x > 0.0){
		return pow(x,1.0/3.0);
	}else{
		return -pow(-x,1.0/3.0);
	}
}

double complex Squreroot(double complex x){
	double complex y;
	double r = sqrt(creal(x)*creal(x)+cimag(x)*cimag(x)),theta = atan2(cimag(x),creal(x));
	if(cimag(x) == 0.0){
		if(creal(x) > 0.0){
			y = sqrt(r);
		}else{
			y = sqrt(r)*1i;
		}
	}else{
		if(theta < 0.0){
			theta += 2.0*M_PI;
		}
		double complex y = sqrt(r)*(cos(theta*0.5)+1i*sin(theta*0.5));
	}
	return y;
}

以上を実行すると以下のような出力が得られます。

==Cubic Equation==
-1.650629 + 0.000000 i
-0.174685 + 1.546869 i
-0.174685 + -1.546869 i

==Quartic Equation==
-1.287815 + 0.857897 i
-1.287815 + -0.857897 i
0.287815 + 1.416093 i
0.287815 + -1.416093 i

確かに正しい解が出力されています。いろいろな値を入れて確かめると面白いと思います。もしうまくいかないような例があれば報告していただけるとありがたいです。

#1 解の公式
まずは三次方程式・四次方程式の解の公式をそれぞれ示します。
##1.1 三次方程式
実数係数の三次方程式$ax^3+bx^2+cx+d = 0$を考えることにします。
このとき以下のようにパラメータを導入します：
$A = b/a,~ B = c/a,~ C = d/a,~ p = B-\dfrac{A^2}{3},~ q = \dfrac{2A^3}{27}-\dfrac{AB}{3}+C$
また判別式$D = \dfrac{q^2}{4}+\dfrac{p^3}{27}$を導入します。判別式の値に応じて以下のようになります：

(1) $D>0$のとき
このとき一つの実数解と二つの複素共役な複素数解を持ちます。二次方程式の判別式とは異なり、値が正の時に複素数解を持つことに注意してください。
実数$u = \sqrt[3]{-\dfrac{q}{2}+\sqrt{D} },~ v = \sqrt[3]{-\dfrac{q}{2}-\sqrt{D}}$を利用すると、三次方程式の解$x_1,x_2,x_3$は以下で与えられます：

\begin{eqnarray*}
x_1 &=& u+v-\dfrac{A}{3} \\
x_2 &=& -\dfrac{1}{2}(u+v)+\dfrac{\sqrt3 i}{2}(u-v)-\dfrac{A}{3} \\
x_3 &=& -\dfrac{1}{2}(u+v)-\dfrac{\sqrt3 i}{2}(u-v)-\dfrac{A}{3}
\end{eqnarray*}

(2) $D < 0$のとき
このとき三つの実数解を持ちます。ここで複素数$-\dfrac{q}{2}+i\sqrt{-D}$の偏角$\theta$を導入すれば、先ほど定義したパラメータ$p$を用いることで三次方程式の解が求まります：

\begin{eqnarray*}
x_1 &=& 2 \sqrt{ -\dfrac{p}{3} } \cos \dfrac{\theta}{3} \\
x_2 &=& 2 \sqrt{ -\dfrac{p}{3} } \cos \dfrac{\theta+2\pi}{3} \\
x_3 &=& 2 \sqrt{ -\dfrac{p}{3} } \cos \dfrac{\theta+4\pi}{3} 
\end{eqnarray*}

ただし判別式$D$の定義と$D<0$から$p<0$となるので、確かに$x_1,x_2,x_3$は実数となることに注意してください。実数解であるにも関わらず複素数を経由しなければならないということは興味深いです。

(3) $D = 0$のとき
このとき1つの実数解と、実数の重解を持ちます。三次関数のグラフで考えると、極値が$x$軸に接しているような状況になります。三次方程式の解は以下になります：

\begin{eqnarray*}
x_1 &=& -2 \sqrt[3]{\dfrac{p}{2} }  \\
x_2 &=& x_3 =   \sqrt[3]{\dfrac{p}{2} }   
\end{eqnarray*}

ただし実装にあたっては数値誤差の影響により$D=0$を正確に判定できないので特別な処理をしていません。それでも(1)や(2)の特殊な場合に含まれるので、解の値としては(数値誤差を除き)正しいものが得られます。

#1.2 四次方程式
実数係数の三次方程式$ax^4+bx^3+cx^2+dx + e = 0$を考えることにします。
このとき以下のようにパラメータを導入します：
$A = b/a,~ B = c/a,~ C = d/a,~ D = e/a$
$p = -6(A/4)^2+B,~q = 8(A/4)^3-2B(A/4)+C,~r = -3(A/4)^3+B(A/4)^2-C(A/4)+D$
また$t^3-pt^2-4rt+(4pr-q^2) = 0$という三次方程式を考えて、実数解$t$を利用することにします。なお、前述のように三次方程式は必ず一つ以上の実数解を持ちます。このとき$t-p$も実数になるので、$m = \sqrt{t-p}$というパラメータを導入すればこれは実数または純虚数になっています。
以上のパラメータを利用すれば、四次方程式の解は以下のように記述されます：

\begin{eqnarray*}
x_1 &=& \dfrac{1}{2} \left( -m+ \sqrt{ -t-p+\dfrac{2q}{m} } \right)-\dfrac{A}{4}  \\
x_2 &=& \dfrac{1}{2} \left( -m- \sqrt{ -t-p+\dfrac{2q}{m} } \right)-\dfrac{A}{4}  \\
x_3 &=& \dfrac{1}{2} \left( m+ \sqrt{ -t-p-\dfrac{2q}{m} } \right)-\dfrac{A}{4}  \\
x_4 &=& \dfrac{1}{2} \left( m-\sqrt{ -t-p-\dfrac{2q}{m} } \right)-\dfrac{A}{4}  \\ 
\end{eqnarray*}

#2 Cによる実装
前述の解の公式を利用して、それぞれの場合に解を計算するプログラムを示します。
#2.1 複素数の記述
実はCにおいてはcomplex.hをincludeすることにより複素数を扱うことができるようになります。

complex_sample.c


#include <stdio.h>
#include <complex.h> //これにより複素数が利用できるようになる

int main(void){
	double a = 1.0,b = 2.0,c = 1.0,d = -1.0;
	double complex z;//複素数の宣言
	double complex x[2];//複素数の配列の宣言はこのようにする
	x[0] = a+b*1i;//1iが虚数単位iを表す
	x[1] = c+d*1i;
	z = x[0]*x[1];//通常の積をするだけで複素数としての積になる
	printf("z = %lf + %lf i\n",creal(z),cimag(z));//crealとcimagで実部と虚部を取り出せる
	return 0;
}

また、sqrt()を複素数に拡張した関数Squrerootを以下に示します。

double complex Squreroot(double complex x){
	double complex y;
	double r = sqrt(creal(x)*creal(x)+cimag(x)*cimag(x)),theta = atan2(cimag(x),creal(x));
	if(cimag(x) == 0.0){//xが実数のとき
		if(creal(x) > 0.0){//xが正なら普通のsqrt()を返す
			y = sqrt(r);
		}else{//xが負なら純虚数を返す
			y = sqrt(r)*1i;
		}
	}else{
		if(theta < 0.0){//この処理がないとtheta*0.5がうまくいかないはず
			theta += 2.0*M_PI;
		}
		double complex y = sqrt(r)*(cos(theta*0.5)+1i*sin(theta*0.5));
	}
	return y;
}

この関数は四次方程式を解く関数の実装の際に必要になります。

#2.2 三次方程式を解く関数の実装
以下の関数SolveCubicEquationが三次方程式を解く関数です。仕様としては、三次方程式$ax^3+bx^2+cx+d = 0$の係数$a,b,c,d$を入力とし、複素数の配列x[3]にこの三次方程式の解を格納します。なお、x[0]は必ず実数になるようになっていることがわかると思います。また実数に対し三乗根を計算するCuberootという関数を定義しています。これは負の値xに対してpow(x,1.0/3.0)がうまく動かないためです。

void SolveCubicEquation(double complex x[3],double a,double b,double c,double d){
	if(a == 0.0){
		printf("Error:a = 0.0\n");
		printf("This equation is NOT Cubic.\n");
	}
	double A,B,C,p,q,D;
	A = b/a;
	B = c/a;
	C = d/a;
	p = B-A*A/3.0;
	q = 2.0*A*A*A/27.0-A*B/3.0+C;
	D = q*q/4.0+p*p*p/27.0;
	if(D < 0.0){//three real solutions
		double theta = atan2(sqrt(-D),-q*0.5);
		x[0] = 2.0*sqrt(-p/3.0)*cos(theta/3.0)-A/3.0;
		x[1] = 2.0*sqrt(-p/3.0)*cos((theta+2.0*M_PI)/3.0)-A/3.0;
		x[2] = 2.0*sqrt(-p/3.0)*cos((theta+4.0*M_PI)/3.0)-A/3.0;
	}else{//single real solution and two imaginary solutions(c.c)
		double u = Cuberoot(-q*0.5+sqrt(D)),v = Cuberoot(-q*0.5-sqrt(D));
		x[0] = u+v-A/3.0;
		x[1] = -0.5*(u+v)+sqrt(3.0)*0.5*1i*(u-v)-A/3.0;
		x[2] = -0.5*(u+v)-sqrt(3.0)*0.5*1i*(u-v)-A/3.0;
	}
}

double Cuberoot(double x){//負の数に対しても三乗根を計算することができる
	if(x > 0.0){
		return pow(x,1.0/3.0);
	}else{
		return -pow(-x,1.0/3.0);
	}
}

#2.3 四次方程式を解く関数の実装
以下の関数SolveQuarticEquationが三次方程式を解く関数です。仕様としては、四次方程式$ax^4+bx^3+cx^2+dx+e = 0$の係数$a,b,c,d,e$を入力とし、複素数の配列x[4]にこの四次方程式の解を格納します。

void SolveQuarticEquation(double complex x[4],double a,double b,double c,double d,double e){
	if(a == 0.0){
		printf("Error:a = 0.0\n");
		printf("This equation is NOT Quartic.\n");
	}
	double A,B,C,D,p,q,r,epsilon_temp = 0.000001;
	A = b/a;
	B = c/a;
	C = d/a;
	D = e/a;
	p = -6.0*pow(A/4.0,2.0)+B;
	q = 8.0*pow(A/4.0,3.0)-2.0*B*A/4.0+C;
	r = -3.0*pow(A/4.0,4.0)+B*pow(A/4.0,2.0)-C*A/4.0+D;
	double complex t_temp[3];
	SolveCubicEquation(t_temp,1.0,-p,-4.0*r,4.0*p*r-q*q);
	double t = creal(t_temp[0]);
	double complex m = Squreroot(t-p);
	x[0] = (-m+Squreroot(-t-p+2.0*q/m))*0.5-A/4.0;
	x[1] = (-m-Squreroot(-t-p+2.0*q/m))*0.5-A/4.0;
	x[2] = (m+Squreroot(-t-p-2.0*q/m))*0.5-A/4.0;
	x[3] = (m-Squreroot(-t-p-2.0*q/m))*0.5-A/4.0;
}

#3 終わりに
いくつかの入力に対しては正しい解が得られることを確かめましたが、コーナーケースに対しての処理は甘い可能性はあるので、もし問題があればご連絡ください。
　　