用盛金公式求解一元三次方程

 

    解一元三次方程一般用盛金公式求解，算法高效且求出来的解精确。
    百度百科关于盛金公式有如下解释：

盛金公式
　　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;
}