代码改变世界

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

2011-08-20 19:49  menggucaoyuan  阅读(929)  评论(0编辑  收藏

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

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