解一元三次方程一般用盛金公式求解，算法高效且求出来的解精确。
百度百科关于盛金公式有如下解释：
盛金公式
Shengjin's Formulas
一元三次方程aX^3+bX^2+cX+d=0，（a，b，c，d∈R，且a≠0）。
重根判别式：A=b^2－3ac；B=bc－9ad；C=c^2－3bd，
总判别式：Δ=B^2－4AC。
当A=B=0时，盛金公式①：
X⑴=X⑵=X⑶=－b/(3a)=－c/b=－3d/c。
当Δ=B^2－4AC>0时，盛金公式②：
X⑴=(－b－Y⑴^(1/3)－Y⑵^(1/3))/(3a)；
X(2,3)=(－2b+Y⑴^(1/3)+Y⑵^(1/3))/(6a)±i3^(1/2)(Y⑴^(1/3)－Y⑵^(1/3))/(6a)，
其中Y(1, 2)=Ab+3a(－B±(B^2－4AC)^(1/2))/2，i^2=－1。
当Δ=B^2－4AC=0时，盛金公式③：
X⑴=－b/a+K；X⑵=X⑶=－K/2，
其中K=B/A，(A≠0)。
当Δ=B^2－4AC<0时，盛金公式④：
X⑴=(－b－2A^(1/2)cos(θ/3))/(3a)；
X(2,3)=(－b+A^(1/2)(cos(θ/3)±3^(1/2)sin(θ/3)))/(3a)，
其中θ=arccosT，T=(2Ab－3aB)/(2A^(3/2))，(A>0，－1<T<1)。
盛金判别法
Shengjin's Distinguishing Means
①：当A=B=0时，方程有一个三重实根；
②：当Δ=B^2－4AC>0时，方程有一个实根和一对共轭复根；
③：当Δ=B^2－4AC=0时，方程有三个实根，其中有一个两重根；
④：当Δ=B^2－4AC<0时，方程有三个不相等的实根。
个人实现的代码实现如下：
#include <cmath>
using namespace std;
#define ZERO 1e-6
#define ISZERO(Value) ((-1.0f * ZERO < (Value)) && ((Value) < 1.0f * ZERO))
#define LESSTHANZERO(Value) ((Value) < ZERO)
#define GREATERTHANZERO(Value) (ZERO < (Value))
#define ISBETWEEN(Value, Value1, Value2)
(GREATERTHAN(Value, MIN(Value1, Value2)) && LESSTHAN(Value,
MAX(Value1, Value2)))
typedef double RealNum;
typedef struct tagComplex
{
RealNum real;
RealNum image;
}Complex;
//rA * (X ^3) + rB * (X ^ 2) + rC * (X) + rD = 0
int Shengjin(RealNum rA, RealNum rB, RealNum rC, RealNum rD, RealNum* rX1, RealNum* rX2, RealNum* rX3)
{
RealNum fDA = rB * rB - 3.0f * rA * rC;
RealNum fDB = rB * rC - 9.0f * rA * rD;
RealNum fDC = rC * rC - 3.0f * rB * rD;
RealNum rTempA = ABSZERO;
RealNum rTempB = ABSZERO;
RealNum fDelta = fDB * fDB - 4.0f * fDA * fDC;
if (ISZERO(fDA) && ISZERO(fDB))
{
if (!ISZERO(rA))
{
*rX1 = -1.0f * rB / (3.0f * rA);
}
else if (!ISZERO(rB))
{
*rX1 = -1.0f * rC / rB;
}
else if (!ISZERO(rC))
{
*rX1 = -1.0f * rD / rC;
}
else// if (!ISZERO(rD))
{
return 0;
}
*rX2 = *rX1;
*rX3 = *rX1;
return 1;
}
if (GREATERTHANZERO(fDelta))
{
RealNum rY1 = fDA * rB + 3.0f * rA * (0.5f * (-1.0f * fDB + sqrtf(fDelta) ));
RealNum rY2 = fDA * rB + 3.0f * rA * (0.5f * (-1.0f * fDB - sqrtf(fDelta) ));
if (GREATERTHANZERO(rY1))
{
rTempA = powf(rY1, 1.0f / 3.0f);
}
else
{
rTempA = -1.0f * powf(-1.0f * rY1, 1.0f / 3.0f);
}
if (GREATERTHANZERO(rY2))
{
rTempB = powf(rY2, 1.0f / 3.0f);
}
else
{
rTempB = -1.0f * powf(-1.0f * rY2, 1.0f / 3.0f);
}
*rX1 = -1.0f * (rB + rTempA + rTempB) * (1.0f / (3.0f * rA));
Complex cX2;
cX2.real = (1.0f / (6.0f * rA)) * (-2.0f * rB + rTempA + rTempB);
cX2.image = (1.0f / (6.0f * rA)) * sqrtf(3.0f) * (rTempA -rTempB);
Complex cX3;
cX3.real = cX2.real;
cX3.image = -1.0f * cX2.image;
return 1;
}
else if (ISZERO(fDelta))
{
if (ISZERO(fDA))
{
return 0;
}
RealNum rK = fDB / fDA;
*rX1 = -1.0f * rB / rA + rK;
*rX2 = -0.5f * rK;
*rX3 = *rX2;
return 2;
}
else
{
if (LESSTHANZERO(fDA))
{
return 0;
}
if (GREATERTHANZERO(fDA))
{
rTempA = powf(fDA, 3.0f);
}
else
{
rTempA = -1.0f * powf(-1.0f * fDA, 3.0f);
}
RealNum rT = (2.0f * fDA * rB - 3.0f * rA * fDB) / (2.0f * sqrtf(rTempA));
if (!ISBETWEEN(rT, -1.0f, 1.0f))
{
return 0;
}
RealNum rTheta = acosf(rT);
*rX1 = -1.0f * (rB + 2.0f * sqrtf(fDA) * cos(rTheta / 3.0f)) * (1.0f / (3.0f * rA));
*rX2 = (-1.0f * rB + sqrtf(fDA) * (cos(rTheta / 3.0f) + sqrtf(3.0f) * sin(rTheta / 3.0f))) * (1.0f / (3.0f * rA));
*rX3 = (-1.0f * rB + sqrtf(fDA) * (cos(rTheta / 3.0f) - sqrtf(3.0f) * sin(rTheta / 3.0f))) * (1.0f / (3.0f * rA));
return 3;
}
return 0;
}