libs/math/test/test_roots.cpp

   1 //  (C) Copyright John Maddock 2006.
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 #include <pch.hpp>
   7
   8 #define BOOST_TEST_MAIN
   9 #include <boost/test/unit_test.hpp>
  10 #include <boost/test/tools/floating_point_comparison.hpp>
  11 #include <boost/test/results_collector.hpp>
  12 #include <boost/math/special_functions/beta.hpp>
  13 #include <boost/math/distributions/skew_normal.hpp>
  14 #include <boost/math/tools/polynomial.hpp>
  15 #include <boost/math/tools/roots.hpp>
  16 #include <boost/math/constants/constants.hpp>
  17 #include <boost/test/results_collector.hpp>
  18 #include <boost/test/unit_test.hpp>
  19 #include <boost/array.hpp>
  20 #include <boost/type_index.hpp>
  21 #include "table_type.hpp"
  22 #include <iostream>
  23 #include <iomanip>
  24
  25 #include <boost/multiprecision/cpp_bin_float.hpp>
  26 #include <boost/multiprecision/cpp_complex.hpp>
  27
  28 #define BOOST_CHECK_CLOSE_EX(a, b, prec, i) \
  29    {\
  30       unsigned int failures = boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed;\
  31       BOOST_CHECK_CLOSE(a, b, prec); \
  32       if(failures != boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed)\
  33       {\
  34          std::cerr << "Failure was at row " << i << std::endl;\
  35          std::cerr << std::setprecision(35); \
  36          std::cerr << "{ " << data[i][0] << " , " << data[i][1] << " , " << data[i][2];\
  37          std::cerr << " , " << data[i][3] << " , " << data[i][4] << " , " << data[i][5] << " } " << std::endl;\
  38       }\
  39    }
  40
  41
  42 //
  43 // Implement various versions of inverse of the incomplete beta
  44 // using different root finding algorithms, and deliberately "bad"
  45 // starting conditions: that way we get all the pathological cases
  46 // we could ever wish for!!!
  47 //
  48
  49 template <class T, class Policy>
  50 struct ibeta_roots_1   // for first order algorithms
  51 {
  52    ibeta_roots_1(T _a, T _b, T t, bool inv = false, bool neg = false)
  53       : a(_a), b(_b), target(t), invert(inv), neg(neg) {}
  54
  55    T operator()(const T& x)
  56    {
  57       return boost::math::detail::ibeta_imp(a, b, (neg ? -x : x), Policy(), invert, true) - target;
  58    }
  59 private:
  60    T a, b, target;
  61    bool invert, neg;
  62 };
  63
  64 template <class T, class Policy>
  65 struct ibeta_roots_2   // for second order algorithms
  66 {
  67    ibeta_roots_2(T _a, T _b, T t, bool inv = false, bool neg = false)
  68       : a(_a), b(_b), target(t), invert(inv), neg(neg) {}
  69
  70    boost::math::tuple<T, T> operator()(const T& xx)
  71    {
  72       typedef typename boost::math::lanczos::lanczos<T, Policy>::type L;
  73       T x = neg ? -xx : xx;
  74       T f = boost::math::detail::ibeta_imp(a, b, x, Policy(), invert, true) - target;
  75       T f1 = invert ?
  76          -boost::math::detail::ibeta_power_terms(b, a, 1 - x, x, L(), true, Policy())
  77                : boost::math::detail::ibeta_power_terms(a, b, x, 1 - x, L(), true, Policy());
  78       T y = 1 - x;
  79       if(y == 0)
  80          y = boost::math::tools::min_value<T>() * 8;
  81       f1 /= y * x;
  82
  83       // make sure we don't have a zero derivative:
  84       if(f1 == 0)
  85          f1 = (invert ? -1 : 1) * boost::math::tools::min_value<T>() * 64;
  86
  87       return boost::math::make_tuple(f, neg ? -f1 : f1);
  88    }
  89 private:
  90    T a, b, target;
  91    bool invert, neg;
  92 };
  93
  94 template <class T, class Policy>
  95 struct ibeta_roots_3   // for third order algorithms
  96 {
  97    ibeta_roots_3(T _a, T _b, T t, bool inv = false, bool neg = false)
  98       : a(_a), b(_b), target(t), invert(inv), neg(neg) {}
  99
 100    boost::math::tuple<T, T, T> operator()(const T& xx)
 101    {
 102       typedef typename boost::math::lanczos::lanczos<T, Policy>::type L;
 103       T x = neg ? -xx : xx;
 104       T f = boost::math::detail::ibeta_imp(a, b, x, Policy(), invert, true) - target;
 105       T f1 = invert ?
 106                -boost::math::detail::ibeta_power_terms(b, a, 1 - x, x, L(), true, Policy())
 107                : boost::math::detail::ibeta_power_terms(a, b, x, 1 - x, L(), true, Policy());
 108       T y = 1 - x;
 109       if(y == 0)
 110          y = boost::math::tools::min_value<T>() * 8;
 111       f1 /= y * x;
 112       T f2 = f1 * (-y * a + (b - 2) * x + 1) / (y * x);
 113       if(invert)
 114          f2 = -f2;
 115
 116       // make sure we don't have a zero derivative:
 117       if(f1 == 0)
 118          f1 = (invert ? -1 : 1) * boost::math::tools::min_value<T>() * 64;
 119
 120       if (neg)
 121       {
 122          f1 = -f1;
 123       }
 124
 125       return boost::math::make_tuple(f, f1, f2);
 126    }
 127 private:
 128    T a, b, target;
 129    bool invert, neg;
 130 };
 131
 132 double inverse_ibeta_bisect(double a, double b, double z)
 133 {
 134    typedef boost::math::policies::policy<> pol;
 135    bool invert = false;
 136    int bits = std::numeric_limits<double>::digits;
 137
 138    //
 139    // special cases, we need to have these because there may be other
 140    // possible answers:
 141    //
 142    if(z == 1) return 1;
 143    if(z == 0) return 0;
 144
 145    //
 146    // We need a good estimate of the error in the incomplete beta function
 147    // so that we don't set the desired precision too high.  Assume that 3-bits
 148    // are lost each time the arguments increase by a factor of 10:
 149    //
 150    using namespace std;
 151    int bits_lost = static_cast<int>(ceil(log10((std::max)(a, b)) * 3));
 152    if(bits_lost < 0)
 153       bits_lost = 3;
 154    else
 155       bits_lost += 3;
 156    int precision = bits - bits_lost;
 157
 158    double min = 0;
 159    double max = 1;
 160    boost::math::tools::eps_tolerance<double> tol(precision);
 161    return boost::math::tools::bisect(ibeta_roots_1<double, pol>(a, b, z, invert), min, max, tol).first;
 162 }
 163
 164 double inverse_ibeta_bisect_neg(double a, double b, double z)
 165 {
 166    typedef boost::math::policies::policy<> pol;
 167    bool invert = false;
 168    int bits = std::numeric_limits<double>::digits;
 169
 170    //
 171    // special cases, we need to have these because there may be other
 172    // possible answers:
 173    //
 174    if(z == 1) return 1;
 175    if(z == 0) return 0;
 176
 177    //
 178    // We need a good estimate of the error in the incomplete beta function
 179    // so that we don't set the desired precision too high.  Assume that 3-bits
 180    // are lost each time the arguments increase by a factor of 10:
 181    //
 182    using namespace std;
 183    int bits_lost = static_cast<int>(ceil(log10((std::max)(a, b)) * 3));
 184    if(bits_lost < 0)
 185       bits_lost = 3;
 186    else
 187       bits_lost += 3;
 188    int precision = bits - bits_lost;
 189
 190    double min = -1;
 191    double max = 0;
 192    boost::math::tools::eps_tolerance<double> tol(precision);
 193    return -boost::math::tools::bisect(ibeta_roots_1<double, pol>(a, b, z, invert, true), min, max, tol).first;
 194 }
 195
 196 double inverse_ibeta_newton(double a, double b, double z)
 197 {
 198    double guess = 0.5;
 199    bool invert = false;
 200    int bits = std::numeric_limits<double>::digits;
 201
 202    //
 203    // special cases, we need to have these because there may be other
 204    // possible answers:
 205    //
 206    if(z == 1) return 1;
 207    if(z == 0) return 0;
 208
 209    //
 210    // We need a good estimate of the error in the incomplete beta function
 211    // so that we don't set the desired precision too high.  Assume that 3-bits
 212    // are lost each time the arguments increase by a factor of 10:
 213    //
 214    using namespace std;
 215    int bits_lost = static_cast<int>(ceil(log10((std::max)(a, b)) * 3));
 216    if(bits_lost < 0)
 217       bits_lost = 3;
 218    else
 219       bits_lost += 3;
 220    int precision = bits - bits_lost;
 221
 222    double min = 0;
 223    double max = 1;
 224    return boost::math::tools::newton_raphson_iterate(ibeta_roots_2<double, boost::math::policies::policy<> >(a, b, z, invert), guess, min, max, precision);
 225 }
 226
 227 double inverse_ibeta_newton_neg(double a, double b, double z)
 228 {
 229    double guess = 0.5;
 230    bool invert = false;
 231    int bits = std::numeric_limits<double>::digits;
 232
 233    //
 234    // special cases, we need to have these because there may be other
 235    // possible answers:
 236    //
 237    if(z == 1) return 1;
 238    if(z == 0) return 0;
 239
 240    //
 241    // We need a good estimate of the error in the incomplete beta function
 242    // so that we don't set the desired precision too high.  Assume that 3-bits
 243    // are lost each time the arguments increase by a factor of 10:
 244    //
 245    using namespace std;
 246    int bits_lost = static_cast<int>(ceil(log10((std::max)(a, b)) * 3));
 247    if(bits_lost < 0)
 248       bits_lost = 3;
 249    else
 250       bits_lost += 3;
 251    int precision = bits - bits_lost;
 252
 253    double min = -1;
 254    double max = 0;
 255    return -boost::math::tools::newton_raphson_iterate(ibeta_roots_2<double, boost::math::policies::policy<> >(a, b, z, invert, true), -guess, min, max, precision);
 256 }
 257
 258 double inverse_ibeta_halley(double a, double b, double z)
 259 {
 260    double guess = 0.5;
 261    bool invert = false;
 262    int bits = std::numeric_limits<double>::digits;
 263
 264    //
 265    // special cases, we need to have these because there may be other
 266    // possible answers:
 267    //
 268    if(z == 1) return 1;
 269    if(z == 0) return 0;
 270
 271    //
 272    // We need a good estimate of the error in the incomplete beta function
 273    // so that we don't set the desired precision too high.  Assume that 3-bits
 274    // are lost each time the arguments increase by a factor of 10:
 275    //
 276    using namespace std;
 277    int bits_lost = static_cast<int>(ceil(log10((std::max)(a, b)) * 3));
 278    if(bits_lost < 0)
 279       bits_lost = 3;
 280    else
 281       bits_lost += 3;
 282    int precision = bits - bits_lost;
 283
 284    double min = 0;
 285    double max = 1;
 286    return boost::math::tools::halley_iterate(ibeta_roots_3<double, boost::math::policies::policy<> >(a, b, z, invert), guess, min, max, precision);
 287 }
 288
 289 double inverse_ibeta_halley_neg(double a, double b, double z)
 290 {
 291    double guess = -0.5;
 292    bool invert = false;
 293    int bits = std::numeric_limits<double>::digits;
 294
 295    //
 296    // special cases, we need to have these because there may be other
 297    // possible answers:
 298    //
 299    if(z == 1) return 1;
 300    if(z == 0) return 0;
 301
 302    //
 303    // We need a good estimate of the error in the incomplete beta function
 304    // so that we don't set the desired precision too high.  Assume that 3-bits
 305    // are lost each time the arguments increase by a factor of 10:
 306    //
 307    using namespace std;
 308    int bits_lost = static_cast<int>(ceil(log10((std::max)(a, b)) * 3));
 309    if(bits_lost < 0)
 310       bits_lost = 3;
 311    else
 312       bits_lost += 3;
 313    int precision = bits - bits_lost;
 314
 315    double min = -1;
 316    double max = 0;
 317    return -boost::math::tools::halley_iterate(ibeta_roots_3<double, boost::math::policies::policy<> >(a, b, z, invert, true), guess, min, max, precision);
 318 }
 319
 320 double inverse_ibeta_schroder(double a, double b, double z)
 321 {
 322    double guess = 0.5;
 323    bool invert = false;
 324    int bits = std::numeric_limits<double>::digits;
 325
 326    //
 327    // special cases, we need to have these because there may be other
 328    // possible answers:
 329    //
 330    if(z == 1) return 1;
 331    if(z == 0) return 0;
 332
 333    //
 334    // We need a good estimate of the error in the incomplete beta function
 335    // so that we don't set the desired precision too high.  Assume that 3-bits
 336    // are lost each time the arguments increase by a factor of 10:
 337    //
 338    using namespace std;
 339    int bits_lost = static_cast<int>(ceil(log10((std::max)(a, b)) * 3));
 340    if(bits_lost < 0)
 341       bits_lost = 3;
 342    else
 343       bits_lost += 3;
 344    int precision = bits - bits_lost;
 345
 346    double min = 0;
 347    double max = 1;
 348    return boost::math::tools::schroder_iterate(ibeta_roots_3<double, boost::math::policies::policy<> >(a, b, z, invert), guess, min, max, precision);
 349 }
 350
 351
 352 template <class Real, class T>
 353 void test_inverses(const T& data)
 354 {
 355    using namespace std;
 356    typedef Real                   value_type;
 357
 358    value_type precision = static_cast<value_type>(ldexp(1.0, 1-boost::math::policies::digits<value_type, boost::math::policies::policy<> >()/2)) * 150;
 359    if(boost::math::policies::digits<value_type, boost::math::policies::policy<> >() < 50)
 360       precision = 1;   // 1% or two decimal digits, all we can hope for when the input is truncated
 361
 362    for(unsigned i = 0; i < data.size(); ++i)
 363    {
 364       //
 365       // These inverse tests are thrown off if the output of the
 366       // incomplete beta is too close to 1: basically there is insuffient
 367       // information left in the value we're using as input to the inverse
 368       // to be able to get back to the original value.
 369       //
 370       if(data[i][5] == 0)
 371       {
 372          BOOST_CHECK_EQUAL(inverse_ibeta_halley(Real(data[i][0]), Real(data[i][1]), Real(data[i][5])), value_type(0));
 373          BOOST_CHECK_EQUAL(inverse_ibeta_halley_neg(Real(data[i][0]), Real(data[i][1]), Real(data[i][5])), value_type(0));
 374          BOOST_CHECK_EQUAL(inverse_ibeta_schroder(Real(data[i][0]), Real(data[i][1]), Real(data[i][5])), value_type(0));
 375          BOOST_CHECK_EQUAL(inverse_ibeta_newton(Real(data[i][0]), Real(data[i][1]), Real(data[i][5])), value_type(0));
 376          BOOST_CHECK_EQUAL(inverse_ibeta_newton_neg(Real(data[i][0]), Real(data[i][1]), Real(data[i][5])), value_type(0));
 377          BOOST_CHECK_EQUAL(inverse_ibeta_bisect(Real(data[i][0]), Real(data[i][1]), Real(data[i][5])), value_type(0));
 378          BOOST_CHECK_EQUAL(inverse_ibeta_bisect_neg(Real(data[i][0]), Real(data[i][1]), Real(data[i][5])), value_type(0));
 379       }
 380       else if((1 - data[i][5] > 0.001)
 381          && (fabs(data[i][5]) > 2 * boost::math::tools::min_value<value_type>())
 382          && (fabs(data[i][5]) > 2 * boost::math::tools::min_value<double>()))
 383       {
 384          value_type inv = inverse_ibeta_halley(Real(data[i][0]), Real(data[i][1]), Real(data[i][5]));
 385          BOOST_CHECK_CLOSE_EX(Real(data[i][2]), inv, precision, i);
 386          inv = inverse_ibeta_halley_neg(Real(data[i][0]), Real(data[i][1]), Real(data[i][5]));
 387          BOOST_ASSERT(boost::math::isfinite(inv));
 388          BOOST_CHECK_CLOSE_EX(Real(data[i][2]), inv, precision, i);
 389          inv = inverse_ibeta_schroder(Real(data[i][0]), Real(data[i][1]), Real(data[i][5]));
 390          BOOST_CHECK_CLOSE_EX(Real(data[i][2]), inv, precision, i);
 391          inv = inverse_ibeta_newton(Real(data[i][0]), Real(data[i][1]), Real(data[i][5]));
 392          BOOST_CHECK_CLOSE_EX(Real(data[i][2]), inv, precision, i);
 393          inv = inverse_ibeta_newton_neg(Real(data[i][0]), Real(data[i][1]), Real(data[i][5]));
 394          BOOST_CHECK_CLOSE_EX(Real(data[i][2]), inv, precision, i);
 395          inv = inverse_ibeta_bisect(Real(data[i][0]), Real(data[i][1]), Real(data[i][5]));
 396          BOOST_CHECK_CLOSE_EX(Real(data[i][2]), inv, precision, i);
 397          inv = inverse_ibeta_bisect_neg(Real(data[i][0]), Real(data[i][1]), Real(data[i][5]));
 398          BOOST_CHECK_CLOSE_EX(Real(data[i][2]), inv, precision, i);
 399       }
 400       else if(1 == data[i][5])
 401       {
 402          BOOST_CHECK_EQUAL(inverse_ibeta_halley(Real(data[i][0]), Real(data[i][1]), Real(data[i][5])), value_type(1));
 403          BOOST_CHECK_EQUAL(inverse_ibeta_halley_neg(Real(data[i][0]), Real(data[i][1]), Real(data[i][5])), value_type(1));
 404          BOOST_CHECK_EQUAL(inverse_ibeta_schroder(Real(data[i][0]), Real(data[i][1]), Real(data[i][5])), value_type(1));
 405          BOOST_CHECK_EQUAL(inverse_ibeta_newton(Real(data[i][0]), Real(data[i][1]), Real(data[i][5])), value_type(1));
 406          BOOST_CHECK_EQUAL(inverse_ibeta_newton_neg(Real(data[i][0]), Real(data[i][1]), Real(data[i][5])), value_type(1));
 407          BOOST_CHECK_EQUAL(inverse_ibeta_bisect(Real(data[i][0]), Real(data[i][1]), Real(data[i][5])), value_type(1));
 408          BOOST_CHECK_EQUAL(inverse_ibeta_bisect_neg(Real(data[i][0]), Real(data[i][1]), Real(data[i][5])), value_type(1));
 409       }
 410
 411    }
 412 }
 413
 414 #ifndef SC_
 415 #define SC_(x) static_cast<typename table_type<T>::type>(BOOST_JOIN(x, L))
 416 #endif
 417
 418 template <class T>
 419 void test_beta(T, const char* /* name */)
 420 {
 421    //
 422    // The actual test data is rather verbose, so it's in a separate file
 423    //
 424    // The contents are as follows, each row of data contains
 425    // five items, input value a, input value b, integration limits x, beta(a, b, x) and ibeta(a, b, x):
 426    //
 427 #  include "ibeta_small_data.ipp"
 428
 429    test_inverses<T>(ibeta_small_data);
 430
 431 #  include "ibeta_data.ipp"
 432
 433    test_inverses<T>(ibeta_data);
 434
 435 #  include "ibeta_large_data.ipp"
 436
 437    test_inverses<T>(ibeta_large_data);
 438 }
 439
 440 #if !defined(BOOST_NO_CXX11_AUTO_DECLARATIONS) && !defined(BOOST_NO_CXX11_UNIFIED_INITIALIZATION_SYNTAX) && !defined(BOOST_NO_CXX11_LAMBDAS)
 441 template <class Complex>
 442 void test_complex_newton()
 443 {
 444     typedef typename Complex::value_type Real;
 445     std::cout << "Testing complex Newton's Method on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 446     using std::abs;
 447     using std::sqrt;
 448     using boost::math::tools::complex_newton;
 449     using boost::math::tools::polynomial;
 450     using boost::math::constants::half;
 451
 452     Real tol = std::numeric_limits<Real>::epsilon();
 453     // p(z) = z^2 + 1, roots: \pm i.
 454     polynomial<Complex> p{{1,0}, {0, 0}, {1,0}};
 455     Complex guess{1,1};
 456     polynomial<Complex> p_prime = p.prime();
 457     auto f = [&](Complex z) { return std::make_pair<Complex, Complex>(p(z), p_prime(z)); };
 458     Complex root = complex_newton(f, guess);
 459
 460     BOOST_CHECK(abs(root.real()) <= tol);
 461     BOOST_CHECK_CLOSE(root.imag(), (Real)1, tol);
 462
 463     guess = -guess;
 464     root = complex_newton(f, guess);
 465     BOOST_CHECK(abs(root.real()) <= tol);
 466     BOOST_CHECK_CLOSE(root.imag(), (Real)-1, tol);
 467
 468     // Test that double roots are handled correctly-as correctly as possible.
 469     // Convergence at a double root is not quadratic.
 470     // This sets p = (z-i)^2:
 471     p = polynomial<Complex>({{-1,0}, {0,-2}, {1,0}});
 472     p_prime = p.prime();
 473     guess = -guess;
 474     auto g = [&](Complex z) { return std::make_pair<Complex, Complex>(p(z), p_prime(z)); };
 475     root = complex_newton(g, guess);
 476     BOOST_CHECK(abs(root.real()) < 10*sqrt(tol));
 477     BOOST_CHECK_CLOSE(root.imag(), (Real)1, tol);
 478
 479     // Test that zero derivatives are handled.
 480     // p(z) = z^2 + iz + 1
 481     p = polynomial<Complex>({{1,0}, {0,1}, {1,0}});
 482     // p'(z) = 2z + i
 483     p_prime = p.prime();
 484     guess = Complex(0,-boost::math::constants::half<Real>());
 485     auto g2 = [&](Complex z) { return std::make_pair<Complex, Complex>(p(z), p_prime(z)); };
 486     root = complex_newton(g2, guess);
 487
 488     // Here's the other root, in case code changes cause it to be found:
 489     //Complex expected_root1{0, half<Real>()*(sqrt(static_cast<Real>(5)) - static_cast<Real>(1))};
 490     Complex expected_root2{0, -half<Real>()*(sqrt(static_cast<Real>(5)) + static_cast<Real>(1))};
 491
 492     BOOST_CHECK_CLOSE(expected_root2.imag(),root.imag(), tol);
 493     BOOST_CHECK(abs(root.real()) < tol);
 494
 495     // Does a zero root pass the termination criteria?
 496     p = polynomial<Complex>({{0,0}, {0,0}, {1,0}});
 497     p_prime = p.prime();
 498     guess = Complex(0, -boost::math::constants::half<Real>());
 499     auto g3 = [&](Complex z) { return std::make_pair<Complex, Complex>(p(z), p_prime(z)); };
 500     root = complex_newton(g3, guess);
 501     BOOST_CHECK(abs(root.real()) < tol);
 502
 503     // Does a monstrous root pass?
 504     Real x = -pow(static_cast<Real>(10), 20);
 505     p = polynomial<Complex>({{x, x}, {1,0}});
 506     p_prime = p.prime();
 507     guess = Complex(0, -boost::math::constants::half<Real>());
 508     auto g4 = [&](Complex z) { return std::make_pair<Complex, Complex>(p(z), p_prime(z)); };
 509     root = complex_newton(g4, guess);
 510     BOOST_CHECK(abs(root.real() + x) < tol);
 511     BOOST_CHECK(abs(root.imag() + x) < tol);
 512
 513 }
 514
 515 // Polynomials which didn't factorize using Newton's method at first:
 516 void test_daubechies_fails()
 517 {
 518     std::cout << "Testing failures from Daubechies filter computation.\n";
 519     using std::abs;
 520     using std::sqrt;
 521     using boost::math::tools::complex_newton;
 522     using boost::math::tools::polynomial;
 523     using boost::math::constants::half;
 524
 525     double tol = 500*std::numeric_limits<double>::epsilon();
 526     polynomial<std::complex<double>> p{{-185961388.136908293,141732493.98435241}, {601080390,0}};
 527     std::complex<double> guess{1,1};
 528     polynomial<std::complex<double>> p_prime = p.prime();
 529     auto f = [&](std::complex<double> z) { return std::make_pair<std::complex<double>, std::complex<double>>(p(z), p_prime(z)); };
 530     std::complex<double> root = complex_newton(f, guess);
 531
 532     std::complex<double> expected_root = -p.data()[0]/p.data()[1];
 533     BOOST_CHECK_CLOSE(expected_root.imag(), root.imag(), tol);
 534     BOOST_CHECK_CLOSE(expected_root.real(), root.real(), tol);
 535 }
 536 #endif
 537
 538 #if !defined(BOOST_NO_CXX17_IF_CONSTEXPR)
 539 template<class Real>
 540 void test_solve_real_quadratic()
 541 {
 542     Real tol = std::numeric_limits<Real>::epsilon();
 543     using boost::math::tools::quadratic_roots;
 544     auto [x0, x1] = quadratic_roots<Real>(1, 0, -1);
 545     BOOST_CHECK_CLOSE(x0, Real(-1), tol);
 546     BOOST_CHECK_CLOSE(x1, Real(1), tol);
 547
 548     auto p = quadratic_roots((Real)7, (Real)0, (Real)0);
 549     BOOST_CHECK_SMALL(p.first, tol);
 550     BOOST_CHECK_SMALL(p.second, tol);
 551
 552     // (x-7)^2 = x^2 - 14*x + 49:
 553     p = quadratic_roots((Real)1, (Real)-14, (Real)49);
 554     BOOST_CHECK_CLOSE(p.first, Real(7), tol);
 555     BOOST_CHECK_CLOSE(p.second, Real(7), tol);
 556
 557     // This test does not pass in multiprecision,
 558     // due to the fact it does not have an fma:
 559     if (std::is_floating_point<Real>::value)
 560     {
 561         // (x-1)(x-1-eps) = x^2 + (-eps - 2)x + (1)(1+eps)
 562         Real eps = 2*std::numeric_limits<Real>::epsilon();
 563         Real b = 256 * (-2 - eps);
 564         Real c = 256 * (1 + eps);
 565         p = quadratic_roots((Real)256, b, c);
 566         BOOST_CHECK_CLOSE(p.first, Real(1), tol);
 567         BOOST_CHECK_CLOSE(p.second, Real(1) + eps, tol);
 568     }
 569
 570     if (std::is_same<Real, double>::value)
 571     {
 572         // Kahan's example: This is the test that demonstrates the necessity of the fma instruction.
 573         // https://en.wikipedia.org/wiki/Loss_of_significance#Instability_of_the_quadratic_equation
 574         p = quadratic_roots<Real>((Real)94906265.625, (Real )-189812534, (Real)94906268.375);
 575         BOOST_CHECK_CLOSE_FRACTION(p.first, Real(1), tol);
 576         BOOST_CHECK_CLOSE_FRACTION(p.second, 1.000000028975958, 4*tol);
 577     }
 578 }
 579
 580 template<class Z>
 581 void test_solve_int_quadratic()
 582 {
 583     double tol = std::numeric_limits<double>::epsilon();
 584     using boost::math::tools::quadratic_roots;
 585     auto [x0, x1] = quadratic_roots(1, 0, -1);
 586     BOOST_CHECK_CLOSE(x0, double(-1), tol);
 587     BOOST_CHECK_CLOSE(x1, double(1), tol);
 588
 589     auto p = quadratic_roots(7, 0, 0);
 590     BOOST_CHECK_SMALL(p.first, tol);
 591     BOOST_CHECK_SMALL(p.second, tol);
 592
 593     // (x-7)^2 = x^2 - 14*x + 49:
 594     p = quadratic_roots(1, -14, 49);
 595     BOOST_CHECK_CLOSE(p.first, double(7), tol);
 596     BOOST_CHECK_CLOSE(p.second, double(7), tol);
 597 }
 598
 599 template<class Complex>
 600 void test_solve_complex_quadratic()
 601 {
 602     using Real = typename Complex::value_type;
 603     Real tol = std::numeric_limits<Real>::epsilon();
 604     using boost::math::tools::quadratic_roots;
 605     auto [x0, x1] = quadratic_roots<Complex>({1,0}, {0,0}, {-1,0});
 606     BOOST_CHECK_CLOSE(x0.real(), Real(-1), tol);
 607     BOOST_CHECK_CLOSE(x1.real(), Real(1), tol);
 608     BOOST_CHECK_SMALL(x0.imag(), tol);
 609     BOOST_CHECK_SMALL(x1.imag(), tol);
 610
 611     auto p = quadratic_roots<Complex>({7,0}, {0,0}, {0,0});
 612     BOOST_CHECK_SMALL(p.first.real(), tol);
 613     BOOST_CHECK_SMALL(p.second.real(), tol);
 614
 615     // (x-7)^2 = x^2 - 14*x + 49:
 616     p = quadratic_roots<Complex>({1,0}, {-14,0}, {49,0});
 617     BOOST_CHECK_CLOSE(p.first.real(), Real(7), tol);
 618     BOOST_CHECK_CLOSE(p.second.real(), Real(7), tol);
 619
 620 }
 621
 622 #endif
 623
 624 void test_failures()
 625 {
 626 #if !defined(BOOST_NO_CXX11_LAMBDAS)
 627    // There is no root:
 628    BOOST_CHECK_THROW(boost::math::tools::newton_raphson_iterate([](double x) { return std::make_pair(x * x + 1, 2 * x); }, 10.0, -12.0, 12.0, 52), boost::math::evaluation_error);
 629    BOOST_CHECK_THROW(boost::math::tools::newton_raphson_iterate([](double x) { return std::make_pair(x * x + 1, 2 * x); }, -10.0, -12.0, 12.0, 52), boost::math::evaluation_error);
 630    // There is a root, but a bad guess takes us into a local minima:
 631    BOOST_CHECK_THROW(boost::math::tools::newton_raphson_iterate([](double x) { return std::make_pair(boost::math::pow<6>(x) - 2 * boost::math::pow<4>(x) + x + 0.5, 6 * boost::math::pow<5>(x) - 8 * boost::math::pow<3>(x) + 1); }, 0.75, -20., 20., 52), boost::math::evaluation_error);
 632
 633    // There is no root:
 634    BOOST_CHECK_THROW(boost::math::tools::halley_iterate([](double x) { return std::make_tuple(x * x + 1, 2 * x, 2); }, 10.0, -12.0, 12.0, 52), boost::math::evaluation_error);
 635    BOOST_CHECK_THROW(boost::math::tools::halley_iterate([](double x) { return std::make_tuple(x * x + 1, 2 * x, 2); }, -10.0, -12.0, 12.0, 52), boost::math::evaluation_error);
 636    // There is a root, but a bad guess takes us into a local minima:
 637    BOOST_CHECK_THROW(boost::math::tools::halley_iterate([](double x) { return std::make_tuple(boost::math::pow<6>(x) - 2 * boost::math::pow<4>(x) + x + 0.5, 6 * boost::math::pow<5>(x) - 8 * boost::math::pow<3>(x) + 1, 30 * boost::math::pow<4>(x) - 24 * boost::math::pow<2>(x)); }, 0.75, -20., 20., 52), boost::math::evaluation_error);
 638 #endif
 639 }
 640
 641 BOOST_AUTO_TEST_CASE( test_main )
 642 {
 643
 644    test_beta(0.1, "double");
 645
 646 #if !defined(BOOST_NO_CXX11_AUTO_DECLARATIONS) && !defined(BOOST_NO_CXX11_UNIFIED_INITIALIZATION_SYNTAX) && !defined(BOOST_NO_CXX11_LAMBDAS)
 647    test_complex_newton<std::complex<float>>();
 648    test_complex_newton<std::complex<double>>();
 649    test_complex_newton<boost::multiprecision::cpp_complex_100>();
 650    test_daubechies_fails();
 651 #endif
 652
 653 #if !defined(BOOST_NO_CXX17_IF_CONSTEXPR)
 654     test_solve_real_quadratic<float>();
 655     test_solve_real_quadratic<double>();
 656     test_solve_real_quadratic<long double>();
 657     test_solve_real_quadratic<boost::multiprecision::cpp_bin_float_50>();
 658
 659     test_solve_int_quadratic<int>();
 660     test_solve_complex_quadratic<std::complex<double>>();
 661 #endif
 662     test_failures();
 663 }