已知一元四次方程 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 (21​bx)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+21​bx)2=(41​b2−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+21​bx)y+41​y2,可得： [ ( 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+21​bx)+21​y]2=(41​b2−c+y)x2+(21​by−d)x+41​y2−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 (21​by−d)2−4(41​b2−c+y)(41​y2−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 41​b2−c+y=0，根据式(5)，一次项系数 1 2 b y − d \dfrac{1}{2}by-d 21​by−d也必然为0。当 1 4 y 2 − e < 0 \dfrac{1}{4}y^2-e return (x >= -FLT_MAX && x 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) 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) 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) *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([email protected]) *************************************************/ 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 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) 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) 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 x[0] = (-b + pow(-Y1, 1.0 / 3.0) - pow(Y2, 1.0 / 3.0)) / (3.0*a); } else if (Y1 > 0.0 && Y2 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) 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 return errNo; } } if (T 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([email protected]) *************************************************/ 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 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) 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 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); } 四、参考资料

费拉里与一元四次方程的解法