用盛金公式解三次方程（ansi c版）

ansi風格： 

/*
   cc cubic.c -lm
   gcc cubic.c -lm

Shengjin's Formulas Univariate cubic equation aX ^ 3 + bX ^ 2 + cX + d = 0, (a, b, c, d < R, and a!= 0).
       Multiple root discriminant: delta1 = b^2-3*a*c; delta2 = b*c-9*a*d; delta3 = c^2-3*b*d,
           The total discriminant is delta=delta2^2-4*delta1*delta3.
           When delta1 = delta2 = 0, Shengjin Formula (1): X1=X2=X3=-b/(3*a)=-c/b=-3d/c.
           When delta=B^2-4*A*C>0, Shengjin Formula II:
                 Y1= delta1*b + 3*a *((-B + (delta)^1/2))/ 2.
                 Y2= delta1*b + 3*a *((-B - (delta)^1/2))/ 2.
                 x1=(-b-Y1^(1/3) - Y1^(1/3)/(3*a);
                 X2=(-2*b+Y1^(1/3)+Y2^(1/3)/(6*a)+3^(1/2)* (Y1^(1/3)-Y2^(1/3)/(6a)i,
                 X3=(-2*b+Y1^(1/3)+Y2^(1/3)/(6*a)-3^(1/2)* (Y1^(1/3)-Y2^(1/3)/(6a)i,

           When delta=B^2-4AC=0, Shengjin Formula 3:
                 X1=-b/a+K; X2=X3=-K/2,            K = delta2/delta1, (A<>0).
           When delta=B^2-4AC<0, Shengjin Formula 4:
                 X1= (-b-2*sqrt(delta1)*cos(theta/3))/(3*a);
                 X2= (-b+sqrt(delta1)*(cos(theta/3)+sqrt(3)*sin(theta/3)))/(3*a);
                 X3= (-b+sqrt(delta1)*(cos(theta/3)-sqrt(3)*sin(theta/3)))/(3*a)
                 theta=arccosT,T=(2Ab-3aB)/(2A^(3/2))
    Shengjin's Distinguishing Means
        (1)A = B = 0, the equation has a triple real root.
        (2) When delta =B^2-4AC>0, the equation has a real root and a pair of conjugate complex roots.
        (3) When delta=B^2-4AC=0, the equation has three real roots, one of which has two double roots.
        (4) When delta=B^2-4AC<0, the equation has three unequal real roots.
*/

#include <stdio.h>
#include <math.h>
enum roottype{UniReal,OneRPairComplex,TwoReal,UnequalReal,NotCubic};

int main(){
    const double PI=3.14159265359;
    double a,b,c,d;    /*Coefficient of cubic Equation*/
    double delta1,delta2,delta3,delta;
    double Y1,Y2,expY1,expY2;
    double K,theta,T;
    enum roottype rt;

    printf("Input Coefficient of cubic Equation:a b c d.\n");
    scanf("%lf%lf%lf%lf",&a,&b,&c,&d);

    delta1 = b*b-3*a*c;
    delta2 = b*c-9*a*d;
    delta3 = c*c-3*b*d;
    delta  = delta2*delta2-4*delta1*delta3;

    if(delta1==0&&delta2==0) rt=UniReal;
    else if(delta >0 )  rt=OneRPairComplex;
    else if(delta ==0)  rt=TwoReal;
    else if(delta <0 )  rt=UnequalReal;
    if (a==0) rt=NotCubic;
    switch(rt){
      case UniReal:
         printf("the equation has a triple real root.\n");
         printf("x1=x2=x3=%f\n",-b/3/a);
         break;
      case OneRPairComplex:
         printf("the equation has a real root and a pair of conjugate complex roots.\n");
         Y1 = delta1*b + 3*a *((-delta2 + sqrt(delta)))/2;
         Y2 = delta1*b + 3*a *((-delta2 - sqrt(delta)))/2;

         if (Y1>0)  expY1 = pow(Y1,1.0/3.0); else expY1 = (-1)*pow(fabs(Y1),1.0/3.0);
         if (Y2>0)  expY2 = pow(Y2,1.0/3.0); else expY2 = (-1)*pow(fabs(Y2),1.0/3.0);
         printf("x1=%f\n",(-b-expY1-expY2)/(3*a));
         printf("x2=%f+%fi\n",(-2*b+expY1+expY2)/(6*a),sqrt(3.0)* (expY1-expY2)/(6*a));
         printf("x3=%f-%fi\n",(-2*b+expY1+expY2)/(6*a),sqrt(3.0)* (expY1-expY2)/(6*a));
         break;
      case TwoReal:
         printf("the equation has three real roots, one of which has two double roots.\n");
         K =delta2/delta1;
         printf("X1=X2=%f\n",-K/2);
         printf("X3=%f\n",-b/a+K);
         break;
      case UnequalReal:
         printf("the equation has three unequal real roots.\n");
         T=(2.0*delta1*b-3.0*a*delta2)/(2.0*sqrt(delta1*delta1*delta1));
         theta=acos(T);
         printf("X1=%f\n",(-b-2.0*sqrt(delta1)*cos(theta/3.0))/(3.0*a));
         printf("X2=%f\n",(-b+sqrt(delta1)*(cos(theta/3.0)+sqrt(3.0)*sin(theta/3.0)))/(3.0*a));
         printf("X3=%f\n",(-b+sqrt(delta1)*(cos(theta/3.0)-sqrt(3.0)*sin(theta/3.0)))/(3.0*a));
         break;
      case NotCubic:
         printf("the equation is  not a cubic.\n");
         break;
    } /* end of switch*/

   return(0);
}

K&R風格

#include <stdio.h>
#include <math.h>
#define PI 3.14159265359

enum roottype{UniReal,OneRPairComplex,TwoReal,UnequalReal,NotCubic};

main(){
    double a,b,c,d;    /*Coefficient of cubic Equation*/
    double delta1,delta2,delta3,delta;
    double Y1,Y2,expY1,expY2;
    double K,theta,T;
    enum roottype rt;

    printf("Input Coefficient of cubic Equation:a b c d.\n");
    scanf("%lf%lf%lf%lf",&a,&b,&c,&d);

    delta1 = b*b-3*a*c;
    delta2 = b*c-9*a*d;
    delta3 = c*c-3*b*d;
    delta  = delta2*delta2-4*delta1*delta3;

    if(delta1==0&&delta2==0) rt=UniReal;
    else if(delta >0 )  rt=OneRPairComplex;
    else if(delta ==0)  rt=TwoReal;
    else if(delta <0 )  rt=UnequalReal;
    if (a==0) rt=NotCubic;

    switch(rt){
      case UniReal:
         printf("the equation has a triple real root.\n");
         printf("x1=x2=x3=%f\n",-b/3/a);
         break;
      case OneRPairComplex:
         printf("the equation has a real root and a pair of conjugate complex roots.\n");
         Y1 = delta1*b + 3*a *((-delta2 + sqrt(delta)))/2;
         Y2 = delta1*b + 3*a *((-delta2 - sqrt(delta)))/2;

         if (Y1>0)  expY1 = pow(Y1,1.0/3.0); else expY1 = (-1)*pow(fabs(Y1),1.0/3.0);
         if (Y2>0)  expY2 = pow(Y2,1.0/3.0); else expY2 = (-1)*pow(fabs(Y2),1.0/3.0);
         /*printf("expY1=%f,expY2=%f,x3=%f\n",expY1,expY2,a);*/		 
         printf("x1=%f\n",(-b-expY1-expY2)/(3*a));
         printf("x2=%f+%fi\n",(-2*b+expY1+expY2)/(6.0*a),sqrt(3.0)* (expY1-expY2)/(6.0*a));
         printf("x3=%f-%fi\n",(-2*b+expY1+expY2)/(6.0*a),sqrt(3.0)* (expY1-expY2)/(6.0*a));
         break;
      case TwoReal:
         printf("the equation has three real roots, one of which has two double roots.\n");
         K =delta2/delta1;
         printf("X1=X2=%f\n",-K/2);
         printf("X3=%f\n",-b/a+K);
         break;
      case UnequalReal:
         printf("the equation has three unequal real roots.\n");
         T=(2.0*delta1*b-3.0*a*delta2)/(2.0*sqrt(delta1*delta1*delta1));
         theta=acos(T);
         printf("X1=%f\n",(-b-2.0*sqrt(delta1)*cos(theta/3.0))/(3.0*a));
         printf("X2=%f\n",(-b+sqrt(delta1)*(cos(theta/3.0)+sqrt(3.0)*sin(theta/3.0)))/(3.0*a));
         printf("X3=%f\n",(-b+sqrt(delta1)*(cos(theta/3.0)-sqrt(3.0)*sin(theta/3.0)))/(3.0*a));
         break;
      case NotCubic:
         printf("the equation is  not a cubic.\n");
         break;
    } /* end of switch*/

}

參考

matlab實現一元三次方程求解