libs/math/test/test_root_iterations.cpp

   1 //  (C) Copyright John Maddock 2015.
   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 #ifndef BOOST_NO_CXX11_HDR_TUPLE
   9
  10 #define BOOST_TEST_MAIN
  11 #include <boost/test/unit_test.hpp>
  12 #include <boost/test/tools/floating_point_comparison.hpp>
  13 #include <boost/math/tools/roots.hpp>
  14 #include <boost/test/results_collector.hpp>
  15 #include <boost/test/unit_test.hpp>
  16 #include <boost/math/special_functions/cbrt.hpp>
  17 #include <boost/math/special_functions/beta.hpp>
  18 #include <iostream>
  19 #include <iomanip>
  20 #include <tuple>
  21 #include "table_type.hpp"
  22
  23 // No derivatives - using TOMS748 internally.
  24 struct cbrt_functor_noderiv
  25 { //  cube root of x using only function - no derivatives.
  26    cbrt_functor_noderiv(double to_find_root_of) : a(to_find_root_of)
  27    { // Constructor just stores value a to find root of.
  28    }
  29    double operator()(double x)
  30    {
  31       double fx = x*x*x - a; // Difference (estimate x^3 - a).
  32       return fx;
  33    }
  34 private:
  35    double a; // to be 'cube_rooted'.
  36 }; // template <class T> struct cbrt_functor_noderiv
  37
  38 // Using 1st derivative only Newton-Raphson
  39 struct cbrt_functor_deriv
  40 { // Functor also returning 1st derviative.
  41    cbrt_functor_deriv(double const& to_find_root_of) : a(to_find_root_of)
  42    { // Constructor stores value a to find root of,
  43       // for example: calling cbrt_functor_deriv<double>(x) to use to get cube root of x.
  44    }
  45    std::pair<double, double> operator()(double const& x)
  46    { // Return both f(x) and f'(x).
  47       double fx = x*x*x - a; // Difference (estimate x^3 - value).
  48       double dx = 3 * x*x; // 1st derivative = 3x^2.
  49       return std::make_pair(fx, dx); // 'return' both fx and dx.
  50    }
  51 private:
  52    double a; // to be 'cube_rooted'.
  53 };
  54 // Using 1st and 2nd derivatives with Halley algorithm.
  55 struct cbrt_functor_2deriv
  56 { // Functor returning both 1st and 2nd derivatives.
  57    cbrt_functor_2deriv(double const& to_find_root_of) : a(to_find_root_of)
  58    { // Constructor stores value a to find root of, for example:
  59       // calling cbrt_functor_2deriv<double>(x) to get cube root of x,
  60    }
  61    std::tuple<double, double, double> operator()(double const& x)
  62    { // Return both f(x) and f'(x) and f''(x).
  63       double fx = x*x*x - a; // Difference (estimate x^3 - value).
  64       double dx = 3 * x*x; // 1st derivative = 3x^2.
  65       double d2x = 6 * x; // 2nd derivative = 6x.
  66       return std::make_tuple(fx, dx, d2x); // 'return' fx, dx and d2x.
  67    }
  68 private:
  69    double a; // to be 'cube_rooted'.
  70 };
  71
  72 template <class T, class Policy>
  73 struct ibeta_roots_1   // for first order algorithms
  74 {
  75    ibeta_roots_1(T _a, T _b, T t, bool inv = false)
  76       : a(_a), b(_b), target(t), invert(inv) {}
  77
  78    T operator()(const T& x)
  79    {
  80       return boost::math::detail::ibeta_imp(a, b, x, Policy(), invert, true) - target;
  81    }
  82 private:
  83    T a, b, target;
  84    bool invert;
  85 };
  86
  87 template <class T, class Policy>
  88 struct ibeta_roots_2   // for second order algorithms
  89 {
  90    ibeta_roots_2(T _a, T _b, T t, bool inv = false)
  91       : a(_a), b(_b), target(t), invert(inv) {}
  92
  93    boost::math::tuple<T, T> operator()(const T& x)
  94    {
  95       typedef boost::math::lanczos::lanczos<T, Policy> S;
  96       typedef typename S::type L;
  97       T f = boost::math::detail::ibeta_imp(a, b, x, Policy(), invert, true) - target;
  98       T f1 = invert ?
  99          -boost::math::detail::ibeta_power_terms(b, a, 1 - x, x, L(), true, Policy())
 100          : boost::math::detail::ibeta_power_terms(a, b, x, 1 - x, L(), true, Policy());
 101       T y = 1 - x;
 102       if (y == 0)
 103          y = boost::math::tools::min_value<T>() * 8;
 104       f1 /= y * x;
 105
 106       // make sure we don't have a zero derivative:
 107       if (f1 == 0)
 108          f1 = (invert ? -1 : 1) * boost::math::tools::min_value<T>() * 64;
 109
 110       return boost::math::make_tuple(f, f1);
 111    }
 112 private:
 113    T a, b, target;
 114    bool invert;
 115 };
 116
 117 template <class T, class Policy>
 118 struct ibeta_roots_3   // for third order algorithms
 119 {
 120    ibeta_roots_3(T _a, T _b, T t, bool inv = false)
 121       : a(_a), b(_b), target(t), invert(inv) {}
 122
 123    boost::math::tuple<T, T, T> operator()(const T& x)
 124    {
 125       typedef typename boost::math::lanczos::lanczos<T, Policy>::type L;
 126       T f = boost::math::detail::ibeta_imp(a, b, x, Policy(), invert, true) - target;
 127       T f1 = invert ?
 128          -boost::math::detail::ibeta_power_terms(b, a, 1 - x, x, L(), true, Policy())
 129          : boost::math::detail::ibeta_power_terms(a, b, x, 1 - x, L(), true, Policy());
 130       T y = 1 - x;
 131       if (y == 0)
 132          y = boost::math::tools::min_value<T>() * 8;
 133       f1 /= y * x;
 134       T f2 = f1 * (-y * a + (b - 2) * x + 1) / (y * x);
 135       if (invert)
 136          f2 = -f2;
 137
 138       // make sure we don't have a zero derivative:
 139       if (f1 == 0)
 140          f1 = (invert ? -1 : 1) * boost::math::tools::min_value<T>() * 64;
 141
 142       return boost::math::make_tuple(f, f1, f2);
 143    }
 144 private:
 145    T a, b, target;
 146    bool invert;
 147 };
 148
 149
 150 BOOST_AUTO_TEST_CASE( test_main )
 151 {
 152    int newton_limits = static_cast<int>(std::numeric_limits<double>::digits * 0.6);
 153
 154    double arg = 1e-50;
 155    boost::uintmax_t iters;
 156    double guess;
 157    double dr;
 158
 159    while(arg < 1e50)
 160    {
 161       double result = boost::math::cbrt(arg);
 162       //
 163       // Start with a really bad guess 5 times below the result:
 164       //
 165       guess = result / 5;
 166       iters = 1000;
 167       // TOMS algo first:
 168       std::pair<double, double> r = boost::math::tools::bracket_and_solve_root(cbrt_functor_noderiv(arg), guess, 2.0, true, boost::math::tools::eps_tolerance<double>(), iters);
 169       BOOST_CHECK_CLOSE_FRACTION((r.first + r.second) / 2, result, std::numeric_limits<double>::epsilon() * 4);
 170       BOOST_CHECK_LE(iters, 14);
 171       // Newton next:
 172       iters = 1000;
 173       dr = boost::math::tools::newton_raphson_iterate(cbrt_functor_deriv(arg), guess, guess / 2, result * 10, newton_limits, iters);
 174       BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
 175       BOOST_CHECK_LE(iters, 12);
 176       // Halley next:
 177       iters = 1000;
 178       dr = boost::math::tools::halley_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters);
 179       BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
 180       BOOST_CHECK_LE(iters, 7);
 181       // Schroder next:
 182       iters = 1000;
 183       dr = boost::math::tools::schroder_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters);
 184       BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
 185       BOOST_CHECK_LE(iters, 11);
 186       //
 187       // Over again with a bad guess 5 times larger than the result:
 188       //
 189       iters = 1000;
 190       guess = result * 5;
 191       r = boost::math::tools::bracket_and_solve_root(cbrt_functor_noderiv(arg), guess, 2.0, true, boost::math::tools::eps_tolerance<double>(), iters);
 192       BOOST_CHECK_CLOSE_FRACTION((r.first + r.second) / 2, result, std::numeric_limits<double>::epsilon() * 4);
 193       BOOST_CHECK_LE(iters, 14);
 194       // Newton next:
 195       iters = 1000;
 196       dr = boost::math::tools::newton_raphson_iterate(cbrt_functor_deriv(arg), guess, result / 10, result * 10, newton_limits, iters);
 197       BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
 198       BOOST_CHECK_LE(iters, 12);
 199       // Halley next:
 200       iters = 1000;
 201       dr = boost::math::tools::halley_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters);
 202       BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
 203       BOOST_CHECK_LE(iters, 7);
 204       // Schroder next:
 205       iters = 1000;
 206       dr = boost::math::tools::schroder_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters);
 207       BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
 208       BOOST_CHECK_LE(iters, 11);
 209       //
 210       // A much better guess, 1% below result:
 211       //
 212       iters = 1000;
 213       guess = result * 0.9;
 214       r = boost::math::tools::bracket_and_solve_root(cbrt_functor_noderiv(arg), guess, 2.0, true, boost::math::tools::eps_tolerance<double>(), iters);
 215       BOOST_CHECK_CLOSE_FRACTION((r.first + r.second) / 2, result, std::numeric_limits<double>::epsilon() * 4);
 216       BOOST_CHECK_LE(iters, 12);
 217       // Newton next:
 218       iters = 1000;
 219       dr = boost::math::tools::newton_raphson_iterate(cbrt_functor_deriv(arg), guess, result / 10, result * 10, newton_limits, iters);
 220       BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
 221       BOOST_CHECK_LE(iters, 5);
 222       // Halley next:
 223       iters = 1000;
 224       dr = boost::math::tools::halley_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters);
 225       BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
 226       BOOST_CHECK_LE(iters, 3);
 227       // Schroder next:
 228       iters = 1000;
 229       dr = boost::math::tools::schroder_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters);
 230       BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
 231       BOOST_CHECK_LE(iters, 4);
 232       //
 233       // A much better guess, 1% above result:
 234       //
 235       iters = 1000;
 236       guess = result * 1.1;
 237       r = boost::math::tools::bracket_and_solve_root(cbrt_functor_noderiv(arg), guess, 2.0, true, boost::math::tools::eps_tolerance<double>(), iters);
 238       BOOST_CHECK_CLOSE_FRACTION((r.first + r.second) / 2, result, std::numeric_limits<double>::epsilon() * 4);
 239       BOOST_CHECK_LE(iters, 12);
 240       // Newton next:
 241       iters = 1000;
 242       dr = boost::math::tools::newton_raphson_iterate(cbrt_functor_deriv(arg), guess, result / 10, result * 10, newton_limits, iters);
 243       BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
 244       BOOST_CHECK_LE(iters, 5);
 245       // Halley next:
 246       iters = 1000;
 247       dr = boost::math::tools::halley_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters);
 248       BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
 249       BOOST_CHECK_LE(iters, 3);
 250       // Schroder next:
 251       iters = 1000;
 252       dr = boost::math::tools::schroder_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters);
 253       BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
 254       BOOST_CHECK_LE(iters, 4);
 255
 256       arg *= 3.5;
 257    }
 258
 259    //
 260    // Test ibeta as this triggers all the pathological cases!
 261    //
 262 #ifndef SC_
 263 #define SC_(x) x
 264 #endif
 265 #define T double
 266
 267 #  include "ibeta_small_data.ipp"
 268
 269    for (unsigned i = 0; i < ibeta_small_data.size(); ++i)
 270    {
 271       //
 272       // These inverse tests are thrown off if the output of the
 273       // incomplete beta is too close to 1: basically there is insuffient
 274       // information left in the value we're using as input to the inverse
 275       // to be able to get back to the original value.
 276       //
 277       if (ibeta_small_data[i][5] == 0)
 278       {
 279          iters = 1000;
 280          dr = boost::math::tools::newton_raphson_iterate(ibeta_roots_2<double, boost::math::policies::policy<> >(ibeta_small_data[i][0], ibeta_small_data[i][1], ibeta_small_data[i][5]), 0.5, 0.0, 1.0, 53, iters);
 281          BOOST_CHECK_EQUAL(dr, 0.0);
 282          BOOST_CHECK_LE(iters, 27);
 283          iters = 1000;
 284          dr = boost::math::tools::halley_iterate(ibeta_roots_3<double, boost::math::policies::policy<> >(ibeta_small_data[i][0], ibeta_small_data[i][1], ibeta_small_data[i][5]), 0.5, 0.0, 1.0, 53, iters);
 285          BOOST_CHECK_EQUAL(dr, 0.0);
 286          BOOST_CHECK_LE(iters, 10);
 287       }
 288       else if ((1 - ibeta_small_data[i][5] > 0.001)
 289          && (fabs(ibeta_small_data[i][5]) > 2 * boost::math::tools::min_value<double>()))
 290       {
 291          iters = 1000;
 292          double result = ibeta_small_data[i][2];
 293          dr = boost::math::tools::newton_raphson_iterate(ibeta_roots_2<double, boost::math::policies::policy<> >(ibeta_small_data[i][0], ibeta_small_data[i][1], ibeta_small_data[i][5]), 0.5, 0.0, 1.0, 53, iters);
 294          BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 200);
 295 #if defined(BOOST_MSVC) && (BOOST_MSVC == 1600)
 296          BOOST_CHECK_LE(iters, 40);
 297 #else
 298          BOOST_CHECK_LE(iters, 27);
 299 #endif
 300          iters = 1000;
 301          result = ibeta_small_data[i][2];
 302          dr = boost::math::tools::halley_iterate(ibeta_roots_3<double, boost::math::policies::policy<> >(ibeta_small_data[i][0], ibeta_small_data[i][1], ibeta_small_data[i][5]), 0.5, 0.0, 1.0, 53, iters);
 303          BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 200);
 304          BOOST_CHECK_LE(iters, 40);
 305       }
 306       else if (1 == ibeta_small_data[i][5])
 307       {
 308          iters = 1000;
 309          dr = boost::math::tools::newton_raphson_iterate(ibeta_roots_2<double, boost::math::policies::policy<> >(ibeta_small_data[i][0], ibeta_small_data[i][1], ibeta_small_data[i][5]), 0.5, 0.0, 1.0, 53, iters);
 310          BOOST_CHECK_EQUAL(dr, 1.0);
 311          BOOST_CHECK_LE(iters, 27);
 312          iters = 1000;
 313          dr = boost::math::tools::halley_iterate(ibeta_roots_3<double, boost::math::policies::policy<> >(ibeta_small_data[i][0], ibeta_small_data[i][1], ibeta_small_data[i][5]), 0.5, 0.0, 1.0, 53, iters);
 314          BOOST_CHECK_EQUAL(dr, 1.0);
 315          BOOST_CHECK_LE(iters, 10);
 316       }
 317    }
 318
 319 }
 320
 321 #else
 322
 323 int main() { return 0; }
 324
 325 #endif