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

wangxuwei

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实现一元三次方程求解

### wangxuwei

10071 - Back to High School Physics

sailtseng
2013/06/04
216
0
Eran/BezierMathUtils

#BezierMathUtils 介绍 BezierMathUtils是用来提供计算贝塞尔曲线时候需要用到的数学公式。 比如存在如下需求： 有两点坐标A和B，希望某一个物品从A点移动到B点. 移动方式希望采用贝塞尔曲线...

Eran
2015/02/18
0
0

davelet
2014/05/11
0
0

fcsong000833
2013/01/24
375
1
error C2143: syntax error : missing ';' before ...

dalu~
2012/12/06
297
0

iOS Xcode升级包地址（感谢大神）

_____1____
10分钟前
4
0
Qt编写自定义控件71-圆弧进度条

11分钟前
3
0

22分钟前
4
0
【1015】LNMP架构二

【1015】LNMP架构二 三、PHP安装 PHP安装和LAMP安装PHP方法有差别，需要开启php-fpm服务 1、下载PHP7至/usr/local/src/ 切换目录：cd /usr/local/src 2、解压缩 tar -jxvf php-7.3.0.tar.gz...

56分钟前
5
0

6
0