一、问题描述
已知一元四次方程 a x 4 + b x 3 + c x 2 + d x + e = 0 ( a ≠ 0 ) ax^4 + bx^3 + cx^2+dx+e = 0(a\neq0) ax4+bx3+cx2+dx+e=0(a=0),求方程的实数根。笔者用C语言实现了费拉里法求解一元四次方程高精度的实数根。一元四次方程实数根求解过程会调用一元二次方程和一元三次方程高精度实数根的求解函数,参见另外一篇博文:一元二次方程高精度实数根(C语言),一元三次方程高精度实数根(C语言)
二、算法步骤
将一元四次方程最高次项系数化为1,则:
x
4
+
b
x
3
+
c
x
2
+
d
x
+
e
=
0
(1)
x^4+bx^3+cx^2+dx+e=0 \tag 1
x4+bx3+cx2+dx+e=0(1)
移项得到:
x
4
+
b
x
3
=
−
c
x
2
−
d
x
−
e
(2)
x^4+bx^3=-cx^2-dx-e \tag 2
x4+bx3=−cx2−dx−e(2)
式(2)两边同时加上
(
1
2
b
x
)
2
(\dfrac{1}{2}bx)^2
(21bx)2,可将式(2)左边配成完全平方式:
(
x
2
+
1
2
b
x
)
2
=
(
1
4
b
2
−
c
)
x
2
−
d
x
−
e
(3)
(x^2+\dfrac{1}{2}bx)^2=(\dfrac{1}{4}b^2-c)x^2-dx-e \tag 3
(x2+21bx)2=(41b2−c)x2−dx−e(3)
式(3)两边同时加上
(
x
2
+
1
2
b
x
)
y
+
1
4
y
2
(x^2+\dfrac{1}{2}bx)y+\dfrac{1}{4}y^2
(x2+21bx)y+41y2,可得:
[
(
x
2
+
1
2
b
x
)
+
1
2
y
]
2
=
(
1
4
b
2
−
c
+
y
)
x
2
+
(
1
2
b
y
−
d
)
x
+
1
4
y
2
−
e
(4)
[(x^2+\dfrac{1}{2}bx) + \dfrac{1}{2}y]^2=(\dfrac{1}{4}b^2-c+y)x^2+(\dfrac{1}{2}by-d)x+ \dfrac{1}{4}y^2-e\tag 4
[(x2+21bx)+21y]2=(41b2−c+y)x2+(21by−d)x+41y2−e(4)
式(4)中的y是一个参数。当式(4)中的x为原方程的根时,不论y取什么值,式(4)都应该成立。特别地,如果所取得y值使得式(4)右边关于x的二次多项式也能变成一个完全平方式,则对式(4)两边同时开方可以得到次数较低的方程。
为了使式(4)右边关于x的二次多项式也能变成一个完全平方式,只需使它的判别式变成0,即:
(
1
2
b
y
−
d
)
2
−
4
(
1
4
b
2
−
c
+
y
)
(
1
4
y
2
−
e
)
=
0
(5)
(\dfrac{1}{2}by-d)^2-4(\dfrac{1}{4}b^2-c+y)(\dfrac{1}{4}y^2-e)=0 \tag 5
(21by−d)2−4(41b2−c+y)(41y2−e)=0(5)
化简得到:
y
3
−
c
y
2
+
(
b
d
−
4
e
)
y
+
(
4
c
−
b
2
)
e
−
d
2
=
0
(6)
y^3 -cy^2 + (bd - 4e)y + (4c - b^2)e - d^2=0 \tag 6
y3−cy2+(bd−4e)y+(4c−b2)e−d2=0(6)
式(6)为关于y的一元三次方程,求出其任意一实数根后代入式(4),式(4)两边都成为完全平方式,两边开方,可以得到两个关于x的一元二次方程,解这两个一元二次方程,就得到原方程的实数根。
求解式(6)的一元三次方程得到的实数根后,需要分情况讨论:
a、若式(4)右边关于x的二次项系数
1
4
b
2
−
c
+
y
=
0
\dfrac{1}{4}b^2-c+y=0
41b2−c+y=0,根据式(5),一次项系数
1
2
b
y
−
d
\dfrac{1}{2}by-d
21by−d也必然为0。当
1
4
y
2
−
e
<
0
\dfrac{1}{4}y^2-e <0
41y2−e<0时,原方程无解;当
1
4
y
2
−
e
≥
0
\dfrac{1}{4}y^2-e \ge0
41y2−e≥0时,原方程的解为这两个一元二次方程的解:
x
2
+
1
2
b
x
+
1
2
y
=
±
1
4
y
2
−
e
(7)
x^2+\dfrac{1}{2}bx + \dfrac{1}{2}y = \pm\sqrt{\dfrac{1}{4}y^2-e} \tag 7
x2+21bx+21y=±41y2−e(7)
化简得到:
2
x
2
+
b
x
+
y
±
y
2
−
4
e
=
0
(8)
2x^2+bx + y \pm\sqrt{y^2-4e}= 0 \tag 8
2x2+bx+y±y2−4e=0(8)
b、若式(4)右边关于x的二次项系数
1
4
b
2
−
c
+
y
≠
0
\dfrac{1}{4}b^2-c+y\ne0
41b2−c+y=0,式(4)右边写成:
(
1
4
b
2
−
c
+
y
)
x
2
+
(
1
2
b
y
−
d
)
x
+
1
4
y
2
−
e
=
(
1
4
b
2
−
c
+
y
)
[
x
2
+
1
2
b
y
−
d
1
4
b
2
−
c
+
y
x
+
1
4
y
2
−
e
1
4
b
2
−
c
+
y
]
=
(
1
4
b
2
−
c
+
y
)
[
x
+
1
2
b
y
−
d
2
(
1
4
b
2
−
c
+
y
)
]
2
(9)
(\dfrac{1}{4}b^2-c+y)x^2+(\dfrac{1}{2}by-d)x+ \dfrac{1}{4}y^2-e\\ = (\dfrac{1}{4}b^2-c+y)[x^2+ \dfrac{\dfrac{1}{2}by-d}{\dfrac{1}{4}b^2-c+y}x+\dfrac{\dfrac{1}{4}y^2-e}{\dfrac{1}{4}b^2-c+y} ]\\=(\dfrac{1}{4}b^2-c+y)[x+\dfrac{\dfrac{1}{2}by-d}{2(\dfrac{1}{4}b^2-c+y)}]^2\tag 9
(41b2−c+y)x2+(21by−d)x+41y2−e=(41b2−c+y)[x2+41b2−c+y21by−dx+41b2−c+y41y2−e]=(41b2−c+y)[x+2(41b2−c+y)21by−d]2(9)
注:因为上式中
x
2
+
1
2
b
y
−
d
1
4
b
2
−
c
+
y
x
+
1
4
y
2
−
e
1
4
b
2
−
c
+
y
x^2+ \dfrac{\dfrac{1}{2}by-d}{\dfrac{1}{4}b^2-c+y}x+\dfrac{\dfrac{1}{4}y^2-e}{\dfrac{1}{4}b^2-c+y}
x2+41b2−c+y21by−dx+41b2−c+y41y2−e可以写成完全平方式。
故式(4)可以写成:
[
(
x
2
+
1
2
b
x
)
+
1
2
y
]
2
=
(
1
4
b
2
−
c
+
y
)
[
x
+
1
2
b
y
−
d
2
(
1
4
b
2
−
c
+
y
)
]
2
(10)
[(x^2+\dfrac{1}{2}bx) + \dfrac{1}{2}y]^2=(\dfrac{1}{4}b^2-c+y)[x+\dfrac{\dfrac{1}{2}by-d}{2(\dfrac{1}{4}b^2-c+y)}]^2\tag {10}
[(x2+21bx)+21y]2=(41b2−c+y)[x+2(41b2−c+y)21by−d]2(10)
式(10)两边开方并化简得到两个关于x的一元二次方程:
2
x
2
+
(
b
±
M
)
x
+
y
±
N
M
=
0
(11)
2x^2+(b\pm M)x + y \pm \dfrac{N}{M}= 0 \tag {11}
2x2+(b±M)x+y±MN=0(11)
其中,
M
=
b
2
+
4
(
y
−
c
)
,
N
=
b
y
−
2
d
M=\sqrt{b^2+4(y-c)},N=by-2d
M=b2+4(y−c),N=by−2d。
在计算M的值时,通常取式(6)中一元三次方程的根y,使得
b
2
+
4
(
y
−
c
)
b^2+4(y-c)
b2+4(y−c)最大,以保证根式里值为较大正数(由于计算误差,可能导致
b
2
+
4
(
y
−
c
)
b^2+4(y-c)
b2+4(y−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(double 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(double 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(double p[], double x[], int *rootCount)
{
int i;
double a, b, c, delta, sqrtDelta;
const double ZERO = FLT_MIN; // min normalized positive value(1.175494351e-38F)
const double 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(double p[], double x[], int *rootCount)
{
int i;
double a, b, c, d, A, B, C, delta;
double Y1, Y2, Z1, Z2, K, parm[3], roots[2], theta, T;
const double ZERO = FLT_MIN; // min normalized positive value(1.175494351e-38F)
const double EPS = FLT_MIN;
const double 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;
}
/*************************************************
Function: solve_quartic_equation
Description: 费拉里法求一元四次方程(a*x^4 + b*x^3 + c*x^2 + d*x + e = 0)的所有实数根
Input: 方程的系数 p = {e, d, c, b, a}
Output: 方程的所有实数根x,实数根的个数rootCount
Return: 错误号
Author: Marc Pony(marc_pony@163.com)
*************************************************/
UINT32 solve_quartic_equation(double p[], double x[], int *rootCount)
{
double a, b, c, d, e;
double parm[4], roots[3];
double y, M, N;
double x1[2], x2[2];
int rootCount1, rootCount2, i;
double MSquare, MSquareTemp, temp, yTemp;
const double EPS = FLT_MIN; //min normalized positive value(1.175494351e-38F)
UINT32 errNo = ERR_NO_ERROR;
*rootCount = 0;
for (i = 0; i < 5; i++)
{
if (!is_number(p[i]))
{
errNo = ERR_NAN;
return errNo;
}
if (!is_finite_number(p[i]))
{
errNo = ERR_INF;
return errNo;
}
}
a = p[4];
b = p[3];
c = p[2];
d = p[1];
e = p[0];
if (fabs(a - 0.0) < EPS)
{
if (fabs(b - 0.0) < EPS)
{
parm[2] = c;
parm[1] = d;
parm[0] = e;
errNo = solve_quadratic_equation(parm, x, rootCount);
}
else
{
parm[3] = b;
parm[2] = c;
parm[1] = d;
parm[0] = e;
errNo = solve_cubic_equation(parm, x, rootCount);
}
}
else
{
b /= a;
c /= a;
d /= a;
e /= a;
parm[3] = 1.0;
parm[2] = -c;
parm[1] = b * d - 4.0 * e;
parm[0] = (4 * c - b * b) * e - d * d;
errNo = solve_cubic_equation(parm, roots, rootCount);
if (*rootCount != 0)
{
y = roots[0];
MSquare = b * b + 4.0 * (y - c);
for (i = 1; i < *rootCount; i++)
{
yTemp = roots[i];
MSquareTemp = b * b + 4.0 * (yTemp - c);
if (MSquareTemp > MSquare)
{
MSquare = MSquareTemp;
y = yTemp;
}
}
if (MSquare > 0.0)
{
if (MSquare > 1.0e-8)
{
M = sqrt(MSquare);
N = b * y - 2.0 * d;
parm[2] = 2.0;
parm[1] = b + M;
parm[0] = y + N / M;
errNo = solve_quadratic_equation(parm, x1, &rootCount1);
parm[2] = 2.0;
parm[1] = b - M;
parm[0] = y - N / M;
errNo = solve_quadratic_equation(parm, x2, &rootCount2);
}
else
{
temp = y * y - 4.0 * e;
if (temp >= 0.0)
{
parm[2] = 2.0;
parm[1] = b;
parm[0] = y + sqrt(temp);
errNo = solve_quadratic_equation(parm, x1, &rootCount1);
parm[2] = 2.0;
parm[1] = b;
parm[0] = y - sqrt(temp);
errNo = solve_quadratic_equation(parm, x2, &rootCount2);
}
else
{
*rootCount = 0;
return errNo;
}
}
if (rootCount1 == 2)
{
x[0] = x1[0];
x[1] = x1[1];
x[2] = x2[0];
x[3] = x2[1];
}
else
{
x[0] = x2[0];
x[1] = x2[1];
x[2] = x1[0];
x[3] = x1[1];
}
*rootCount = rootCount1 + rootCount2;
}
else
{
*rootCount = 0;
return errNo;
}
}
}
return errNo;
}
//void main(void)
//{
// double x[2], p[3];
// int rootCount;
// double x1, x2, a, b, c;
// UINT32 errNo = ERR_NO_ERROR;
//
// //一元二次方程测试
// //(1)(x - 1000)*(x - 0.001) = x^2 - 1000.001*x + 1 = 0
// //a = 1;
// //b = -1000.001;
// //c = 1;
//
// //(2) 3*x ^ 2 - 1000000*x + 0 = 0
// //a = 3;
// //b = -1000000;
// //c = 0;
//
// //(3) 1.0e-10*x^2 - 2.0e-10*x + 1.0e-10 = 0
// a = 1.0e-10;
// b = -2.0e-10;
// c = 1.0e-10;
//
// p[0] = c;
// p[1] = b;
// p[2] = a;
// errNo = solve_quadratic_equation(p, x, &rootCount);
// x1 = (-b - sqrt(b * b - 4.0 * a * c)) / (2.0 * a);
// x2 = (-b + sqrt(b * b - 4.0 * a * c)) / (2.0 * a);
//
// printf("原始方法求得方程的根:x1=%f, x2=%f\n本文方法求得方程的根:x1=%f, x2=%f\n",
// MIN(x1, x2), MAX(x1, x2), MIN(x[0], x[1]), MAX(x[0], x[1]));
//}
//void main(void)
//{
// double x[3], p[4];
// int rootCount;
// double 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);
//}
void main(void)
{
double x[4], p[5];
int rootCount;
double a, b, c, d, e;
UINT32 errNo = ERR_NO_ERROR;
//一元四次方程测试
//(1)(x - 1) * (x - 2) * (x^2 + 1) = 0 (x^4 - 3*x^3 + 3*x^2 - 3*x + 2 = 0)
a = 1;
b = -3;
c = 3;
d = -3;
e = 2;
//(2) (x - 1)^2 * (x^2 + 1) = 0 (x^4 - 2*x^3 + 2*x^2 - 2*x + 1 = 0)
//a = 1;
//b = -2;
//c = 2;
//d = -2;
//e = 1;
//(3) (x - 1) * (x - 2) * (x - 3) * (x - 4) = 0 (x^4 - 10*x^3 + 35*x^2 - 50*x + 24 = 0)
//a = 1;
//b = -10;
//c = 35;
//d = -50;
//e = 24;
//(4) (x - 1)^2 * (x - 2)^2 = 0 (x^4 - 6*x^3 + 13*x^2 - 12*x + 4 = 0)
//a = 1;
//b = -6;
//c = 13;
//d = -12;
//e = 4;
//(5) 0*x^4 + x^3 - 3*x^2 + 3*x - 1 = 0
//a = 0;
//b = 1;
//c = -3;
//d = 3;
//e = -1;
//(6) 0*x^4 + 0*x^3 + x^2 - 2*x + 1 = 0
//a = 0;
//b = 0;
//c = 1;
//d = -2;
//e = 1;
//(7) 0*x^4 + 0*x^3 + 0*x^2 - 2*x + 1 = 0
//a = 0;
//b = 0;
//c = 0;
//d = -2;
//e = 1;
p[0] = e;
p[1] = d;
p[2] = c;
p[3] = b;
p[4] = a;
errNo = solve_quartic_equation(p, x, &rootCount);
p[0] = 0.14423869029206515;
p[1] = -0.8769437053260228;
p[2] = 2.0565485222122555;
p[3] = -2.2394836729320553;
p[4] = 0.9623610000000001;
errNo = solve_quartic_equation(p, x, &rootCount);
}