用法拉利方法找出四次方程的真根
我目前正试图用维基百科的法拉利方法解决一个四次方程式。 我想只检索真正的根，丢弃想象的根。 我的实现并没有为真正的根源带来好处。 我找不到公式中的错误。
我的CubicEquation按预期工作，我的双二次求解也是如此。 现在，我只是错过了法拉利的方法，但是我无法让它发挥作用！
这是我的class级：
public class QuarticFunction { private final double a; private final double b; private final double c; private final double d; private final double e; public QuarticFunction(double a, double b, double c, double d, double e) { this.a = a; this.b = b; this.c = c; this.d = d; this.e = e; } public final double[] findRealRoots() { if (a == 0) { return new CubicEquation(b, c, d, e).findRealRoots(); } if (isBiquadratic()) { //This part works as expected return solveUsingBiquadraticMethod(); } return solveUsingFerrariMethodWikipedia(); } private double[] solveUsingFerrariMethodWikipedia() { // http://en.wikipedia.org/wiki/Quartic_function#Ferrari.27s_solution // ERROR: Wrong numbers when Two Real Roots + Two Complex Roots // ERROR: Wrong numbers when Four Real Roots QuarticFunction depressedQuartic = toDepressed(); if (depressedQuartic.isBiquadratic()) { return depressedQuartic.solveUsingBiquadraticMethod(); } double y = findFerraryY(depressedQuartic); double originalRootConversionPart = -b / (4 * a); double firstPart = Math.sqrt(depressedQuartic.c + 2 * y); double positiveSecondPart = Math.sqrt(-(3 * depressedQuartic.c + 2 * y + 2 * depressedQuartic.d / Math.sqrt(depressedQuartic.c + 2 * y))); double negativeSecondPart = Math.sqrt(-(3 * depressedQuartic.c + 2 * y - 2 * depressedQuartic.d / Math.sqrt(depressedQuartic.c + 2 * y))); double x1 = originalRootConversionPart + (firstPart + positiveSecondPart) / 2; double x2 = originalRootConversionPart + (-firstPart + negativeSecondPart) / 2; double x3 = originalRootConversionPart + (firstPart - positiveSecondPart) / 2; double x4 = originalRootConversionPart + (-firstPart - negativeSecondPart) / 2; Set realRoots = findAllRealRoots(x1, x2, x3, x4); return toDoubleArray(realRoots); } private double findFerraryY(QuarticFunction depressedQuartic) { double a3 = 1; double a2 = 5 / 2 * depressedQuartic.c; double a1 = 2 * Math.pow(depressedQuartic.c, 2) - depressedQuartic.e; double a0 = Math.pow(depressedQuartic.c, 3) / 2 - depressedQuartic.c * depressedQuartic.e / 2 - Math.pow(depressedQuartic.d, 2) / 8; //CubicEquation works as expected! No need to worry! It returns either 1 or 3 roots. CubicEquation cubicEquation = new CubicEquation(a3, a2, a1, a0); double[] roots = cubicEquation.findRealRoots(); for (double y : roots) { if (depressedQuartic.c + 2 * y != 0) { return y; } } throw new IllegalStateException("Ferrari method should have at least one y"); } public final boolean isBiquadratic() { return Double.compare(b, 0) == 0 && Double.compare(d, 0) == 0; } private double[] solveUsingBiquadraticMethod() { //It works as expected! QuadraticLine quadraticEquation = new QuadraticLine(a, c, e); if (!quadraticEquation.hasRoots()) { return new double[] {}; } double[] quadraticRoots = quadraticEquation.findRoots(); Set roots = new HashSet(); for (double quadraticRoot : quadraticRoots) { if (quadraticRoot > 0) { roots.add(Math.sqrt(quadraticRoot)); roots.add(-Math.sqrt(quadraticRoot)); } else if (quadraticRoot == 0.00) { roots.add(0.00); } } return toDoubleArray(roots); } public QuarticFunction toDepressed() { // http://en.wikipedia.org/wiki/Quartic_function#Converting_to_a_depressed_quartic double p = (8 * a * c - 3 * Math.pow(b, 2)) / (8 * Math.pow(a, 2)); double q = (Math.pow(b, 3) - 4 * a * b * c + 8 * d * Math.pow(a, 2)) / (8 * Math.pow(a, 3)); double r = (-3 * Math.pow(b, 4) + 256 * e * Math.pow(a, 3) - 64 * d * b * Math.pow(a, 2) + 16 * c * a * Math.pow(b, 2)) / (256 * Math.pow(a, 4)); return new QuarticFunction(1, 0, p, q, r); } private Set findAllRealRoots(double... roots) { Set realRoots = new HashSet(); for (double root : roots) { if (!Double.isNaN(root)) { realRoots.add(root); } } return realRoots; } private double[] toDoubleArray(Collection values) { double[] doubleArray = new double[values.size()]; int i = 0; for (double value : values) { doubleArray[i] = value; ++i; } return doubleArray; } }
我在这里尝试了一个更简单的实现，但是我遇到了一个新问题。 现在，一半的根是好的，但另一半是错的。 这是我的试探性的：
private double[] solveUsingFerrariMethodTheorem() { // https://proofwiki.org/wiki/Ferrari's_Method CubicEquation cubicEquation = getCubicEquationForFerrariMethodTheorem(); double[] cubicRoots = cubicEquation.findRealRoots(); double y1 = findFirstNonZero(cubicRoots); double inRootOfP = Math.pow(b, 2) / Math.pow(a, 2) - 4 * c / a + 4 * y1; double p1 = b / a + Math.sqrt(inRootOfP); double p2 = b / a - Math.sqrt(inRootOfP); double inRootOfQ = Math.pow(y1, 2) - 4 * e / a; double q1 = y1 - Math.sqrt(inRootOfQ); double q2 = y1 + Math.sqrt(inRootOfQ); double x1 = (-p1 + Math.sqrt(Math.pow(p1, 2) - 8 * q1)) / 4; double x2 = (-p2 - Math.sqrt(Math.pow(p2, 2) - 8 * q1)) / 4; double x3 = (-p1 + Math.sqrt(Math.pow(p1, 2) - 8 * q2)) / 4; double x4 = (-p2 - Math.sqrt(Math.pow(p2, 2) - 8 * q2)) / 4; Set realRoots = findAllRealRoots(x1, x2, x3, x4); return toDoubleArray(realRoots); } private CubicEquation getCubicEquationForFerrariMethodTheorem() { double cubicA = 1; double cubicB = -c / a; double cubicC = b * d / Math.pow(a, 2) - 4 * e / a; double cubicD = 4 * c * e / Math.pow(a, 2) - Math.pow(b, 2) * e / Math.pow(a, 3) - Math.pow(d, 2) / Math.pow(a, 2); return new CubicEquation(cubicA, cubicB, cubicC, cubicD); } private double findFirstNonZero(double[] values) { for (double value : values) { if (Double.compare(value, 0) != 0) { return value; } } throw new IllegalArgumentException(values + " does not contain any non-zero value"); }
我不知道我错过了什么。 我花了几个小时调试试图查看错误，但现在我完全迷失了（加上一些令人头疼的事）。 我究竟做错了什么？ 我不关心使用哪些公式，因为我只想让其中一个公式起作用。
我有两个错误。
- 我所有的“整数常量”都应该是双倍的。
- 在维基百科方法中，当凹陷的四次方是双二次时，我们必须将根再转换为原始的四次方。 我正在回归沮丧的根源。
这是我最终的实现
public class QuarticFunction { private static final double NEAR_ZERO = 0.0000001; private final double a; private final double b; private final double c; private final double d; private final double e; public QuarticFunction(double a, double b, double c, double d, double e) { this.a = a; this.b = b; this.c = c; this.d = d; this.e = e; } public final double[] findRealRoots() { if (Math.abs(a) < NEAR_ZERO) { return new CubicFunction(b, c, d, e).findRealRoots(); } if (isBiquadratic()) { return solveUsingBiquadraticMethod(); } return solveUsingFerrariMethodWikipedia(); } private double[] solveUsingFerrariMethodWikipedia() { // http://en.wikipedia.org/wiki/Quartic_function#Ferrari.27s_solution QuarticFunction depressedQuartic = toDepressed(); if (depressedQuartic.isBiquadratic()) { double[] depressedRoots = depressedQuartic.solveUsingBiquadraticMethod(); return reconvertToOriginalRoots(depressedRoots); } double y = findFerraryY(depressedQuartic); double originalRootConversionPart = -b / (4.0 * a); double firstPart = Math.sqrt(depressedQuartic.c + 2.0 * y); double positiveSecondPart = Math.sqrt(-(3.0 * depressedQuartic.c + 2.0 * y + 2.0 * depressedQuartic.d / Math.sqrt(depressedQuartic.c + 2.0 * y))); double negativeSecondPart = Math.sqrt(-(3.0 * depressedQuartic.c + 2.0 * y - 2.0 * depressedQuartic.d / Math.sqrt(depressedQuartic.c + 2.0 * y))); double x1 = originalRootConversionPart + (firstPart + positiveSecondPart) / 2.0; double x2 = originalRootConversionPart + (-firstPart + negativeSecondPart) / 2.0; double x3 = originalRootConversionPart + (firstPart - positiveSecondPart) / 2.0; double x4 = originalRootConversionPart + (-firstPart - negativeSecondPart) / 2.0; Set realRoots = findOnlyRealRoots(x1, x2, x3, x4); return toDoubleArray(realRoots); } private double[] reconvertToOriginalRoots(double[] depressedRoots) { double[] originalRoots = new double[depressedRoots.length]; for (int i = 0; i < depressedRoots.length; ++i) { originalRoots[i] = depressedRoots[i] - b / (4.0 * a); } return originalRoots; } private double findFerraryY(QuarticFunction depressedQuartic) { double a3 = 1.0; double a2 = 5.0 / 2.0 * depressedQuartic.c; double a1 = 2.0 * Math.pow(depressedQuartic.c, 2.0) - depressedQuartic.e; double a0 = Math.pow(depressedQuartic.c, 3.0) / 2.0 - depressedQuartic.c * depressedQuartic.e / 2.0 - Math.pow(depressedQuartic.d, 2.0) / 8.0; CubicFunction cubicFunction = new CubicFunction(a3, a2, a1, a0); double[] roots = cubicFunction.findRealRoots(); for (double y : roots) { if (depressedQuartic.c + 2.0 * y != 0.0) { return y; } } throw new IllegalStateException("Ferrari method should have at least one y"); } private double[] solveUsingBiquadraticMethod() { QuadraticFunction quadraticFunction = new QuadraticFunction(a, c, e); if (!quadraticFunction.hasRoots()) { return new double[] {}; } double[] quadraticRoots = quadraticFunction.findRoots(); Set roots = new HashSet<>(); for (double quadraticRoot : quadraticRoots) { if (quadraticRoot > 0.0) { roots.add(Math.sqrt(quadraticRoot)); roots.add(-Math.sqrt(quadraticRoot)); } else if (quadraticRoot == 0.00) { roots.add(0.00); } } return toDoubleArray(roots); } public QuarticFunction toDepressed() { // http://en.wikipedia.org/wiki/Quartic_function#Converting_to_a_depressed_quartic double p = (8.0 * a * c - 3.0 * Math.pow(b, 2.0)) / (8.0 * Math.pow(a, 2.0)); double q = (Math.pow(b, 3.0) - 4.0 * a * b * c + 8.0 * d * Math.pow(a, 2.0)) / (8.0 * Math.pow(a, 3.0)); double r = (-3.0 * Math.pow(b, 4.0) + 256.0 * e * Math.pow(a, 3.0) - 64.0 * d * b * Math.pow(a, 2.0) + 16.0 * c * a * Math.pow(b, 2.0)) / (256.0 * Math.pow(a, 4.0)); return new QuarticFunction(1.0, 0.0, p, q, r); } private Set findOnlyRealRoots(double... roots) { Set realRoots = new HashSet<>(); for (double root : roots) { if (Double.isFinite(root)) { realRoots.add(root); } } return realRoots; } private double[] toDoubleArray(Collection values) { double[] doubleArray = new double[values.size()]; int i = 0; for (double value : values) { doubleArray[i] = value; ++i; } return doubleArray; } }
你提到的维基百科文章是正确的。 但是，由于有限的浮点算术精度，仅应用公式并不是那么简单。 如果将决定因素的值与零进行比较尤其如此，因为在某些情况下，精度不足会改变符号。
这是我天真的C ++实现（正在进行中），正如您所看到的，它充满了阈值。 我打算使用自适应多精度算术来摆脱未来的门槛。 这在某种程度上适用于[-100, 100]间隔中1-4个实根的方程：
template // lower than above class CQuarticEq { protected: T a; /**< @brief 4th order coefficient */ T b; /**< @brief 3rd order coefficient */ T c; /**< @brief the 2nd order coefficient */ T d; /**< @brief the 1st order coefficient */ T e; /**< @brief 0th order coefficient */ T p_real_root[4]; /**< @brief list of the roots (real parts) */ //T p_im_root[4]; // imaginary part of the roots size_t n_real_root_num; /**< @brief number of real roots */ public: /** * @brief default constructor; solves for roots of \f$ax^4 + bx^3 + cx^2 + dx + e = 0\f$ * * This finds roots of the given equation. It tends to find two identical roots instead of one, rather * than missing one of two different roots - the number of roots found is therefore orientational, * as the roots might have the same value. * * @param[in] _a is 4th order coefficient * @param[in] _b is 3rd order coefficient * @param[in] _c is the 2nd order coefficient * @param[in] _d is the 1st order coefficient * @param[in] _e is constant coefficient */ CQuarticEq(T _a, T _b, T _c, T _d, T _e) :a(_a), b(_b), c(_c), d(_d), e(_e), n_real_root_num(0) { if(fabs(_a) < f_Power_Static(10, n_4th_order_coeff_log10_thresh)) { // otherwise division by a yields large numbers, this is then more precise CCubicEq eq3(_b, _c, _d, _e); n_real_root_num = eq3.n_RealRoot_Num(); for(unsigned int i = 0; i < n_real_root_num; ++ i) { p_real_root[i] = eq3.f_RealRoot(i); //p_im_root[i] = eq3.f_ImagRoot(i); } return; } // the highest power is multiplied by 0, it is a cubic equation if(fabs(_e) == 0) { CCubicEq eq3(_a, _b, _c, _d); n_real_root_num = eq3.n_RealRoot_Num(); for(unsigned int i = 0; i < n_real_root_num; ++ i) { p_real_root[i] = eq3.f_RealRoot(i); //p_im_root[i] = eq3.f_ImagRoot(i); } return; } // the constant is 0, divide the whole equation by x to get a cubic equation if(fabs(b) == 0 && fabs(d) == 0) { // this possibly needs to be done with a threshold CQuadraticEq eq2(a, c, e); for(int i = 0, n = eq2.n_RealRoot_Num(); i < n; ++ i) { T y = eq2.f_RealRoot(i); if(y < 0) { if(y > -f_Power_Static(10, n_zero_discriminant_log10_thresh)) // only slightly negative p_real_root[n_real_root_num ++] = 0; continue; // this root would be imaginary } y = sqrt(y); if(y > 0) p_real_root[n_real_root_num ++] = -y; // generate the roots in approximately sorted order, avoid negative zeros p_real_root[n_real_root_num ++] = y; } // solve a biquadratic equation a*x^4 + c*x^2 + e = 0 } else { #if 0 // this first branch seems to be slightly less precise in the worst case _b /= _a; _c /= _a; _d /= _a; _e /= _a; // normalize T f_roots_offset = _b / 4; T f_alpha = _c - 3.0 / 8 * _b * _b; // p // ok T f_beta = _b * _b * _b / 8 - _b * _c / 2 + _d; // q // ok T f_gamma = -3.0 / 256 * _b * _b * _b * _b + _e - 64.0 / 256 * _d * _b + 16.0 / 256 * _b * _b * _c; // r #else // 0 const T a4 = _a, a3 = _b, a2 = _c, a1 = _d, a0 = _e; // just rename, no normalization T f_roots_offset = a3 / (4 * a4); T f_alpha = (8 * a2 * a4 - 3 * a3 * a3) / (8 * a4 * a4); T f_beta = (a3 * a3 * a3 - 4 * a2 * a3 * a4 + 8 * a1 * a4 * a4) / (8 * a4 * a4 * a4); T f_gamma = (-3 * a3 * a3 * a3 * a3 + 256 * a0 * a4 * a4 * a4 - 64 * a1 * a3 * a4 * a4 + 16 * a2 * a3 * a3 * a4) / (256 * a4 * a4 * a4 * a4); #endif // 0 // convert to a depressed quartic u^4 + f_alpha*u^2 + f_beta*u + f_gamma = 0 (where x = u - f_roots_offset) if(fabs(f_beta) < f_Power_Static(10, n_depressed_1st_order_coeff_log10_thresh)) { CQuadraticEq eq2(1, f_alpha, f_gamma); for(int i = 0, n = eq2.n_RealRoot_Num(); i < n; ++ i) { T y = eq2.f_RealRoot(i); if(y < 0) { if(y > -f_Power_Static(10, n_zero_discriminant_log10_thresh)) // only slightly negative p_real_root[n_real_root_num ++] = 0 - f_roots_offset; continue; // this root would be imaginary } y = sqrt(y); p_real_root[n_real_root_num ++] = -y - f_roots_offset; if(y > 0) p_real_root[n_real_root_num ++] = y - f_roots_offset; // generate the roots in approximately sorted order } // solve a depressed quartic u^4 + f_alpha*u^2 + f_gamma = 0 } else { T f_cub_a = 1, f_cub_b = 5.0 / 2 * f_alpha, f_cub_c = 2 * f_alpha * f_alpha - f_gamma, f_cub_d = (f_alpha * f_alpha * f_alpha - f_alpha * f_gamma) / 2 - f_beta * f_beta / 8; // form a cubic, which gets the parameter y for Ferrari's method CCubicEq eq3(f_cub_a, f_cub_b, f_cub_c, f_cub_d); for(int i = 0, n = eq3.n_RealRoot_Num(); i < n; ++ i) { const T y = eq3.f_RealRoot(i); _ASSERTE(fabs(f_alpha + 2 * y) > 0); // if this triggers, we have apparently missed to detect a biquadratic equation above (would cause a division by zero later on) // get a root of the cubic equation T f_a2y = f_alpha + 2 * y; if(f_a2y < 0) continue; // the root would be an imaginary one f_a2y = (T)sqrt(f_a2y); // get square root of alpha + 2y { T f_det0 = -(3 * f_alpha + 2 * y + 2 * f_beta / f_a2y); if(f_det0 >= 0) { double f_det_sqrt = sqrt(f_det0); p_real_root[n_real_root_num ++] = (f_a2y + f_det_sqrt) / 2 - f_roots_offset; if(f_det_sqrt > 0) // otherwise they are the same p_real_root[n_real_root_num ++] = (f_a2y - f_det_sqrt) / 2 - f_roots_offset; } else if(f_det0 > -f_Power_Static(10, n_zero_discriminant_log10_thresh)) { // only slightly negative double f_det_sqrt = 0; p_real_root[n_real_root_num ++] = (f_a2y + f_det_sqrt) / 2 - f_roots_offset; } } { T f_det1 = -(3 * f_alpha + 2 * y - 2 * f_beta / f_a2y); if(f_det1 >= 0) { double f_det_sqrt = sqrt(f_det1); p_real_root[n_real_root_num ++] = (-f_a2y + f_det_sqrt) / 2 - f_roots_offset; if(f_det_sqrt > 0) // otherwise they are the same p_real_root[n_real_root_num ++] = (-f_a2y - f_det_sqrt) / 2 - f_roots_offset; } else if(f_det1 > -f_Power_Static(10, n_zero_discriminant_log10_thresh)) { // only slightly negative double f_det_sqrt = 0; p_real_root[n_real_root_num ++] = (-f_a2y + f_det_sqrt) / 2 - f_roots_offset; } } break; // the same results are obtained for any value of y } // find roots using Ferrari's method } } if(b_sort_roots) std::sort(p_real_root, p_real_root + n_real_root_num); // sort the roots explicitly (we could use an ordered merge above, // but that would convolute the code somewhat) } /** * @brief gets number of real roots * @return Returns number of real roots (0 to 3). */ size_t n_RealRoot_Num() const { _ASSERTE(n_real_root_num >= 0); return n_real_root_num; } /** * @brief gets value of a real root * @param[in] n_index is zero-based index of the root * @return Returns value of the specified root. */ T f_RealRoot(size_t n_index) const { _ASSERTE(n_index < 4 && n_index < n_real_root_num); return p_real_root[n_index]; } /** * @brief evaluates the equation for a given argument * @param[in] f_x is value of the argument \f$x\f$ * @return Returns value of \f$ax^4 + bx^3 + cx^2 + dx + e\f$. */ T operator ()(T f_x) const { T f_x2 = f_x * f_x; T f_x3 = f_x2 * f_x; T f_x4 = f_x2 * f_x2; return f_x4 * a + f_x3 * b + f_x2 * c + f_x * d + e; } };
您可以将其用作调试实现的起点。 这给出了以下结果：
root response precision 1e-2147483648: 1774587 cases root response precision 1e-24: 5 cases root response precision 1e-23: 1 cases root response precision 1e-22: 10 cases root response precision 1e-21: 9 cases root response precision 1e-20: 16 cases root response precision 1e-19: 43 cases root response precision 1e-18: 84 cases root response precision 1e-17: 136 cases root response precision 1e-16: 447 cases root response precision 1e-15: 97499 cases root response precision 1e-14: 394130 cases root response precision 1e-13: 483321 cases root response precision 1e-12: 435378 cases root response precision 1e-11: 203487 cases root response precision 1e-10: 69011 cases root response precision 1e-9: 22429 cases root response precision 1e-8: 6961 cases root response precision 1e-7: 2368 cases root response precision 1e-6: 744 cases root response precision 1e-5: 208 cases root response precision 1e-4: 24 cases 4-root solution precision 1e-2147483648: 73928 cases 4-root solution precision 1e-23: 1876 cases 4-root solution precision 1e-22: 81123 cases 4-root solution precision 1e-21: 293684 cases 4-root solution precision 1e-20: 745332 cases 4-root solution precision 1e-19: 919578 cases 4-root solution precision 1e-18: 597248 cases 4-root solution precision 1e-17: 268523 cases 4-root solution precision 1e-16: 99841 cases 4-root solution precision 1e-15: 33995 cases 4-root solution precision 1e-14: 11962 cases 4-root solution precision 1e-13: 4366 cases 4-root solution precision 1e-12: 1609 cases 4-root solution precision 1e-11: 547 cases 4-root solution precision 1e-10: 194 cases 4-root solution precision 1e-9: 64 cases 4-root solution precision 1e-8: 4 cases 4-root solution precision 1e-7: 2 cases 3-root solution precision 1e-2147483648: 75013 cases 3-root solution precision 1e-23: 321 cases 3-root solution precision 1e-22: 3281 cases 3-root solution precision 1e-21: 4961 cases 3-root solution precision 1e-20: 5543 cases 3-root solution precision 1e-19: 3248 cases 3-root solution precision 1e-18: 1424 cases 3-root solution precision 1e-17: 654 cases 3-root solution precision 1e-16: 177 cases 3-root solution precision 1e-15: 96 cases 3-root solution precision 1e-14: 372 cases 3-root solution precision 1e-13: 874 cases 3-root solution precision 1e-12: 7575 cases 3-root solution precision 1e-11: 14236 cases 3-root solution precision 1e-10: 11076 cases 3-root solution precision 1e-9: 5020 cases 3-root solution precision 1e-8: 2202 cases 3-root solution precision 1e-7: 1044 cases 3-root solution precision 1e-6: 326 cases 3-root solution precision 1e-5: 39 cases 3-root solution precision 1e-4: 4 cases 3-root solution precision 1e-3: 1 cases 2-root solution precision 1e-2147483648: 35077 cases 2-root solution precision 1e-23: 399 cases 2-root solution precision 1e-22: 2054 cases 2-root solution precision 1e-21: 3016 cases 2-root solution precision 1e-20: 3236 cases 2-root solution precision 1e-19: 2340 cases 2-root solution precision 1e-18: 1337 cases 2-root solution precision 1e-17: 687 cases 2-root solution precision 1e-16: 356 cases 2-root solution precision 1e-15: 3596 cases 2-root solution precision 1e-14: 12063 cases 2-root solution precision 1e-13: 8287 cases 2-root solution precision 1e-12: 7902 cases 2-root solution precision 1e-11: 11479 cases 2-root solution precision 1e-10: 4458 cases 2-root solution precision 1e-9: 1212 cases 2-root solution precision 1e-8: 142 cases 2-root solution precision 1e-7: 21 cases 2-root solution precision 1e-6: 2 cases 2-root solution precision 1e-4: 1 cases 2-root solution precision 1e-3: 1 cases 1-root solution precision 1e-2147483648: 111525 cases 1-root solution precision 1e-23: 867 cases 1-root solution precision 1e-22: 3505 cases 1-root solution precision 1e-21: 2735 cases 1-root solution precision 1e-20: 1888 cases 1-root solution precision 1e-19: 686 cases 1-root solution precision 1e-18: 497 cases 1-root solution precision 1e-17: 138 cases 1-root solution precision 1e-16: 26 cases 1-root solution precision 1e-14: 2 cases had 0 quartic equations with 0 roots had 121869 quartic equations with 1 roots had 48833 quartic equations with 2 roots had 45829 quartic equations with 3 roots had 783469 quartic equations with 4 roots
其中“根响应精度”表示找到的根中的函数的绝对值（ y误差），“根解精度”是找到的根到地面真实根的距离（ x误差）。
请注意，我没有给出CCubicEq和CQuadraticEq ，因为你已经有了这些。