三次方程式の解 - ゲームプログラミングというゲームを攻略するWIKI

三次方程式の解だ。実数解が１つのときはカルダノ式、３つのときはビエト式で計算する。

変数名はなるべく日本語版Wikipediaの「三次方程式」に合わせたので、なんでこれで計算できるかはそっちを読んでくれ。俺に聞くな。

なんかやたら面倒な公式だなと思ったけど、プログラム化したらそれほど複雑ではないな。

#include <math.h>
 
int cubeq(double a3, double a2, double a1, double a0, double * x)
{
	double A2 = a2 / a3;
	double A1 = a1 / a3;
	double A0 = a0 / a3;
	double p = A1 - A2 * A2 / 3;
	double q = A0 - A1 * A2 / 3 + 2 * A2 * A2 * A2 / 27;
	double d = q * q / 4 + p * p * p / 27;
	if (0 <= d)
	{
		double rt = sqrt(d);
		double y = cbrt(-q / 2 + rt) + cbrt(-q / 2 - rt);
		x[0] = y - A2 / 3;
		x[1] = NAN;
		x[2] = NAN;
		return 1;
	}
	else
	{
		double rt = 2 * sqrt(-p / 3);
		double ac = acos(3 * q * sqrt(-3 / p) / (2 * p));
		for (int k = 0; k < 3; k++)
		{
			double y = rt * cos((ac - 2 * M_PI * k) / 3);
			x[k] = y - A2 / 3;
		}
		return 3;
	}
}