java 一元三次方程_用盛金公式解三次方程（ansi c版）

/*

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

#include

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

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;

switch(rt){

case UniReal:

printf("the equation has a triple real root.\n");

printf("x1=x2=x3=%f\n",-c/b);

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)* (expY1-expY2)/(6*a));

printf("x3=%f-%fi\n",(-2*b+expY1+expY2)/(6*a),sqrt(3)* (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)*sin(theta/3.0)))/(3.0*a));

printf("X3=%f\n",(-b+sqrt(delta1)*(cos(theta/3.0)-sqrt(3)*sin(theta/3.0)))/(3.0*a));

break;

} /* end of switch*/

return(0);

}