一、问题描述
已知一元三次方程 a x 3 + b x 2 + c x + d = 0 ( a ≠ 0 ) ax^3 + bx^2 + cx+d = 0(a\neq0) ax3+bx2+cx+d=0(a=0),求方程的实数根。笔者用C语言实现了著名的盛金公式求解一元三次方程高精度的实数根。一元三次方程实数根求解过程会调用一元二次方程高精度实数根的求解函数,参见另外一篇博文:一元二次方程高精度实数根(C语言)。
二、 C C C代码
#include <math.h>
#include <float.h>
#include <stdio.h>
#define UINT32 unsigned int
#define ERR_NO_ERROR 0x00000000
#define ERR_NAN 0x00000001
#define ERR_INF 0x00000002
#define MAX(a, b) ((a) > (b)) ? (a) : (b)
#define MIN(a, b) ((a) < (b)) ? (a) : (b)
/*************************************************
Function: is_number
Description: 判断浮点数是否为nan
Input: 浮点数x
Output: 无
Return: 若浮点数x为nan返回0,否则返回1
Author: Marc Pony(marc_pony@163.com)
*************************************************/
int is_number(float x)
{
return (x == x);
}
/*************************************************
Function: is_finite_number
Description: 判断浮点数是否为inf
Input: 浮点数x
Output: 无
Return: 若浮点数x为inf返回0,否则返回1
Author: Marc Pony(marc_pony@163.com)
*************************************************/
int is_finite_number(float x)
{
return (x >= -FLT_MAX && x <= FLT_MAX);
}
/*************************************************
Function: solve_quadratic_equation
Description: 求一元二次方程(a*x^2 + b*x + c = 0)的所有实数根
Input: 方程的系数 p = {c, b, a}
Output: 方程的所有实数根x, 实数根的个数rootCount
Return: 错误号
Author: Marc Pony(marc_pony@163.com)
*************************************************/
UINT32 solve_quadratic_equation(float p[], float x[], int *rootCount)
{
int i;
float a, b, c, delta, sqrtDelta;
const float ZERO = FLT_MIN; // min normalized positive value(1.175494351e-38F)
const float EPS = FLT_MIN;
UINT32 errNo = ERR_NO_ERROR;
*rootCount = 0;
for (i = 0; i < 3; i++)
{
if (!is_number(p[i]))
{
errNo = ERR_NAN;
return errNo;
}
if (!is_finite_number(p[i]))
{
errNo = ERR_INF;
return errNo;
}
}
a = p[2];
b = p[1];
c = p[0];
if (fabs(a - 0.0) < EPS)
{
if (fabs(b - 0.0) > EPS)
{
x[0] = -c / b;
*rootCount = 1;
}
}
else
{
b /= a;
c /= a;
a = 1.0;
delta = b * b - 4.0 * a * c;
if (delta > ZERO)
{
if (fabs(c - 0.0) < EPS) //若c = 0,由于计算误差,sqrt(b*b - 4*a*c)不等于|b|
{
x[0] = 0.0;
x[1] = -b / a;
}
else
{
sqrtDelta = sqrt(delta);
if (b > 0.0)
{
x[0] = (-2.0 * c) / (b + sqrtDelta); //避免两个很接近的数相减,导致精度丢失
x[1] = (-b - sqrtDelta) / (2.0 * a);
}
else
{
x[0] = (-b + sqrtDelta) / (2.0 * a);
x[1] = (-2.0 * c) / (b - sqrtDelta); //避免两个很接近的数相减,导致精度丢失
}
}
*rootCount = 2;
}
else if (fabs(delta - 0.0) < EPS)
{
x[0] = x[1] = -b / (2.0 * a);
*rootCount = 2;
}
else
{
*rootCount = 0;
}
}
return errNo;
}
/*************************************************
Function: solve_cubic_equation
Description: 盛金公式求一元三次方程(a*x^3 + b*x^2 + c*x + d = 0)的所有实数根
A = b * b - 3.0 * a * c;
B = b * c - 9.0 * a * d;
C = c * c - 3.0 * b * d;
(1)当A = B = 0时,方程有一个三重实根
(2)当Δ = B^2-4 * A * C > 0时,方程有一个实根和一对共轭虚根
(3)当Δ = B^2-4 * A * C = 0时,方程有三个实根,其中有一个两重根
(4)当Δ = B^2-4 * A * C < 0时,方程有三个不相等的实根
Input: 方程的系数 p = {d, c, b, a}
Output: 方程的所有实数根x,实数根的个数rootCount
Return: 错误号
Author: Marc Pony(marc_pony@163.com)
*************************************************/
UINT32 solve_cubic_equation(float p[], float x[], int *rootCount)
{
int i;
float a, b, c, d, A, B, C, delta;
float Y1, Y2, Z1, Z2, K, parm[3], roots[2], theta, T;
const float ZERO = FLT_MIN; // min normalized positive value(1.175494351e-38F)
const float EPS = FLT_MIN;
const float CALCULATE_ERROR = 1.0e-7;
UINT32 errNo = ERR_NO_ERROR;
*rootCount = 0;
for (i = 0; i < 4; i++)
{
if (!is_number(p[i]))
{
errNo = ERR_NAN;
return errNo;
}
if (!is_finite_number(p[i]))
{
errNo = ERR_INF;
return errNo;
}
}
a = p[3];
b = p[2];
c = p[1];
d = p[0];
if (fabs(a - 0.0) < EPS)
{
parm[2] = b;
parm[1] = c;
parm[0] = d;
errNo = solve_quadratic_equation(parm, x, rootCount);
}
else
{
b /= a;
c /= a;
d /= a;
a = 1.0;
A = b * b - 3.0 * a * c;
B = b * c - 9.0 * a * d;
C = c * c - 3.0 * b * d;
delta = B * B - 4.0 * A * C;
if (fabs(A - 0.0) < EPS && fabs(B - 0.0) < EPS)
{
x[0] = x[1] = x[2] = -b / (3.0 * a);
*rootCount = 3;
return errNo;
}
if (delta > ZERO)
{
parm[2] = 1.0;
parm[1] = B;
parm[0] = A * C;
errNo = solve_quadratic_equation(parm, roots, rootCount);
if (errNo != ERR_NO_ERROR)
{
return errNo;
}
Z1 = roots[0];
Z2 = roots[1];
Y1 = A * b + 3.0 * a * Z1;
Y2 = A * b + 3.0 * a * Z2;
if (Y1 < 0.0 && Y2 < 0.0) //pow函数的底数必须为非负数,必须分类讨论
{
x[0] = (-b + pow(-Y1, 1.0 / 3.0) + pow(-Y2, 1.0 / 3.0)) / (3.0*a);
}
else if (Y1 < 0.0 && Y2 > 0.0)
{
x[0] = (-b + pow(-Y1, 1.0 / 3.0) - pow(Y2, 1.0 / 3.0)) / (3.0*a);
}
else if (Y1 > 0.0 && Y2 < 0.0)
{
x[0] = (-b - pow(Y1, 1.0 / 3.0) + pow(-Y2, 1.0 / 3.0)) / (3.0*a);
}
else
{
x[0] = (-b - pow(Y1, 1.0 / 3.0) - pow(Y2, 1.0 / 3.0)) / (3.0*a);
}
*rootCount = 1;
}
else if (fabs(delta - 0.0) < EPS)
{
if (fabs(A - 0.0) > EPS)
{
K = B / A;
x[0] = -b / a + K;
x[1] = x[2] = -0.5 * K;
*rootCount = 3;
}
}
else
{
if (A > 0.0)
{
T = (2.0 * A * b - 3.0 * a * B) / (2.0 * pow(A, 3.0 / 2.0));
if (T > 1.0) //由于计算误差,T的值可能略大于1(如1.0000001)
{
if (T < 1.0 + CALCULATE_ERROR)
{
T = 1.0;
}
else
{
return errNo;
}
}
if (T < -1.0)
{
if (T > -1.0 - CALCULATE_ERROR)
{
T = -1.0;
}
else
{
return errNo;
}
}
theta = acos(T);
x[0] = (-b - 2.0 * sqrt(A) * cos(theta / 3.0)) / (3.0 * a);
x[1] = (-b + sqrt(A) * (cos(theta / 3.0) + sqrt(3.0) * sin(theta / 3.0))) / (3.0 * a);
x[2] = (-b + sqrt(A) * (cos(theta / 3.0) - sqrt(3.0) * sin(theta / 3.0))) / (3.0 * a);
*rootCount = 3;
}
}
}
return errNo;
}
void main(void)
{
float x[3], p[4];
int rootCount;
float a, b, c, d;
UINT32 errNo = ERR_NO_ERROR;
//一元三次方程测试
//(1)(x - 1) * (x^2 + 1) = 0 (x^3 - x^2 + x - 1 = 0)
a = 1;
b = -1;
c = 1;
d = -1;
//(2) (x - 1)^3 = 0 (x^3 - 3*x^2 + 3*x - 1 = 0)
//a = 1;
//b = -3;
//c = 3;
//d = -1;
//(3) (x - 1)^2 * (x - 2) = 0 (x^3 - 4*x^2 + 5*x - 2 = 0)
//a = 1;
//b = -4;
//c = 5;
//d = -2;
//(4) (x - 1) * (x - 2) * (x - 3) = 0 (x^3 - 6*x^2 + 11*x - 6 = 0)
//a = 1;
//b = -6;
//c = 11;
//d = -6;
//(5) 0*x^3 + x^2 - 2*x + 1 = 0
//a = 0;
//b = 1;
//c = -2;
//d = 1;
//(6) 0*x^3 + 0*x^2 - 2*x + 1 = 0
//a = 0;
//b = 0;
//c = -2;
//d = 1;
//(7) 0*x^3 + 0*x^2 + 0*x + 1 = 0
//a = 0;
//b = 0;
//c = 0;
//d = 1;
p[0] = d;
p[1] = c;
p[2] = b;
p[3] = a;
errNo = solve_cubic_equation(p, x, &rootCount);
}