ceph/src/boost/libs/math/test/cubic_roots_test.cpp

   1 /*
   2  * Copyright Nick Thompson, 2021
   3  * Use, modification and distribution are subject to the
   4  * Boost Software License, Version 1.0. (See accompanying file
   5  * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   6  */
   7
   8 #include "math_unit_test.hpp"
   9 #include <boost/math/tools/cubic_roots.hpp>
  10 #include <random>
  11 #ifdef BOOST_HAS_FLOAT128
  12 #include <boost/multiprecision/float128.hpp>
  13 using boost::multiprecision::float128;
  14 #endif
  15
  16 using boost::math::tools::cubic_root_condition_number;
  17 using boost::math::tools::cubic_root_residual;
  18 using boost::math::tools::cubic_roots;
  19 using std::cbrt;
  20
  21 template <class Real> void test_zero_coefficients() {
  22     Real a = 0;
  23     Real b = 0;
  24     Real c = 0;
  25     Real d = 0;
  26     auto roots = cubic_roots(a, b, c, d);
  27     CHECK_EQUAL(roots[0], Real(0));
  28     CHECK_EQUAL(roots[1], Real(0));
  29     CHECK_EQUAL(roots[2], Real(0));
  30
  31     a = 1;
  32     roots = cubic_roots(a, b, c, d);
  33     CHECK_EQUAL(roots[0], Real(0));
  34     CHECK_EQUAL(roots[1], Real(0));
  35     CHECK_EQUAL(roots[2], Real(0));
  36
  37     a = 1;
  38     d = 1;
  39     // x^3 + 1 = 0:
  40     roots = cubic_roots(a, b, c, d);
  41     CHECK_EQUAL(roots[0], Real(-1));
  42     CHECK_NAN(roots[1]);
  43     CHECK_NAN(roots[2]);
  44     d = -1;
  45     // x^3 - 1 = 0:
  46     roots = cubic_roots(a, b, c, d);
  47     CHECK_EQUAL(roots[0], Real(1));
  48     CHECK_NAN(roots[1]);
  49     CHECK_NAN(roots[2]);
  50
  51     d = -2;
  52     // x^3 - 2 = 0
  53     roots = cubic_roots(a, b, c, d);
  54     CHECK_ULP_CLOSE(roots[0], cbrt(Real(2)), 2);
  55     CHECK_NAN(roots[1]);
  56     CHECK_NAN(roots[2]);
  57
  58     d = -8;
  59     roots = cubic_roots(a, b, c, d);
  60     CHECK_EQUAL(roots[0], Real(2));
  61     CHECK_NAN(roots[1]);
  62     CHECK_NAN(roots[2]);
  63
  64     // (x-1)(x-2)(x-3) = x^3 - 6x^2 + 11x - 6
  65     roots = cubic_roots(Real(1), Real(-6), Real(11), Real(-6));
  66     CHECK_ULP_CLOSE(roots[0], Real(1), 2);
  67     CHECK_ULP_CLOSE(roots[1], Real(2), 2);
  68     CHECK_ULP_CLOSE(roots[2], Real(3), 2);
  69
  70     // Double root:
  71     // (x+1)^2(x-2) = x^3 - 3x - 2:
  72     // Note: This test is unstable wrt to perturbations!
  73     roots = cubic_roots(Real(1), Real(0), Real(-3), Real(-2));
  74     CHECK_ULP_CLOSE(Real(-1), roots[0], 2);
  75     CHECK_ULP_CLOSE(Real(-1), roots[1], 2);
  76     CHECK_ULP_CLOSE(Real(2), roots[2], 2);
  77
  78     std::uniform_real_distribution<Real> dis(-2, 2);
  79     std::mt19937 gen(12345);
  80     // Expected roots
  81     std::array<Real, 3> r;
  82     int trials = 10;
  83     for (int i = 0; i < trials; ++i) {
  84         // Mathematica:
  85         // Expand[(x - r0)*(x - r1)*(x - r2)]
  86         // - r0 r1 r2 + (r0 r1 + r0 r2 + r1 r2) x
  87         // - (r0 + r1 + r2) x^2 + x^3
  88         for (auto &root : r) {
  89             root = static_cast<Real>(dis(gen));
  90         }
  91         std::sort(r.begin(), r.end());
  92         Real a = 1;
  93         Real b = -(r[0] + r[1] + r[2]);
  94         Real c = r[0] * r[1] + r[0] * r[2] + r[1] * r[2];
  95         Real d = -r[0] * r[1] * r[2];
  96
  97         auto roots = cubic_roots(a, b, c, d);
  98         // I could check the condition number here, but this is fine right?
  99         if (!CHECK_ULP_CLOSE(r[0], roots[0], 25)) {
 100             std::cerr << "  Polynomial x^3 + " << b << "x^2 + " << c << "x + "
 101                       << d << " has roots {";
 102             std::cerr << r[0] << ", " << r[1] << ", " << r[2]
 103                       << "}, but the computed roots are {";
 104             std::cerr << roots[0] << ", " << roots[1] << ", " << roots[2]
 105                       << "}\n";
 106         }
 107         CHECK_ULP_CLOSE(r[1], roots[1], 25);
 108         CHECK_ULP_CLOSE(r[2], roots[2], 25);
 109         for (auto root : roots) {
 110             auto res = cubic_root_residual(a, b, c, d, root);
 111             CHECK_LE(abs(res[0]), res[1]);
 112         }
 113     }
 114 }
 115
 116 void test_ill_conditioned() {
 117     // An ill-conditioned root reported by SATovstun:
 118     // "Exact" roots produced with a high-precision calcuation on Wolfram Alpha:
 119     // NSolve[x^3 + 10000*x^2 + 200*x +1==0,x]
 120     std::array<double, 3> expected_roots{-9999.97999997,
 121                                          -0.010010015026300100757327057,
 122                                          -0.009990014973799899662674923};
 123     auto roots = cubic_roots<double>(1, 10000, 200, 1);
 124     CHECK_ABSOLUTE_ERROR(expected_roots[0], roots[0],
 125                          std::numeric_limits<double>::epsilon());
 126     CHECK_ABSOLUTE_ERROR(expected_roots[1], roots[1], 1.01e-5);
 127     CHECK_ABSOLUTE_ERROR(expected_roots[2], roots[2], 1.01e-5);
 128     double cond =
 129         cubic_root_condition_number<double>(1, 10000, 200, 1, roots[1]);
 130     double r1 = expected_roots[1];
 131     // The factor of 10 is a fudge factor to make the test pass.
 132     // Nonetheless, it does show this is basically correct:
 133     CHECK_LE(abs(r1 - roots[1]) / abs(r1),
 134              10 * std::numeric_limits<double>::epsilon() * cond);
 135
 136     cond = cubic_root_condition_number<double>(1, 10000, 200, 1, roots[2]);
 137     double r2 = expected_roots[2];
 138     // The factor of 10 is a fudge factor to make the test pass.
 139     // Nonetheless, it does show this is basically correct:
 140     CHECK_LE(abs(r2 - roots[2]) / abs(r2),
 141              10 * std::numeric_limits<double>::epsilon() * cond);
 142
 143     // See https://github.com/boostorg/math/issues/757:
 144     // The polynomial is ((x+1)^2+1)*(x+1) which has roots -1, and two complex
 145     // roots:
 146     roots = cubic_roots<double>(1, 3, 4, 2);
 147     CHECK_ULP_CLOSE(roots[0], -1.0, 3);
 148     CHECK_NAN(roots[1]);
 149     CHECK_NAN(roots[2]);
 150     return;
 151 }
 152
 153 int main() {
 154     test_zero_coefficients<float>();
 155     test_zero_coefficients<double>();
 156 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
 157     test_zero_coefficients<long double>();
 158 #endif
 159     test_ill_conditioned();
 160 #ifdef BOOST_HAS_FLOAT128
 161     // For some reason, the quadmath is way less accurate than the
 162     // float/double/long double:
 163     // test_zero_coefficients<float128>();
 164 #endif
 165
 166     return boost::math::test::report_errors();
 167 }