ceph/src/boost/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/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 <iostream>
  18 #include <iomanip>
  19 #include <tuple>
  20
  21 // No derivatives - using TOMS748 internally.
  22 struct cbrt_functor_noderiv
  23 { //  cube root of x using only function - no derivatives.
  24    cbrt_functor_noderiv(double to_find_root_of) : a(to_find_root_of)
  25    { // Constructor just stores value a to find root of.
  26    }
  27    double operator()(double x)
  28    {
  29       double fx = x*x*x - a; // Difference (estimate x^3 - a).
  30       return fx;
  31    }
  32 private:
  33    double a; // to be 'cube_rooted'.
  34 }; // template <class T> struct cbrt_functor_noderiv
  35
  36 // Using 1st derivative only Newton-Raphson
  37 struct cbrt_functor_deriv
  38 { // Functor also returning 1st derviative.
  39    cbrt_functor_deriv(double const& to_find_root_of) : a(to_find_root_of)
  40    { // Constructor stores value a to find root of,
  41       // for example: calling cbrt_functor_deriv<double>(x) to use to get cube root of x.
  42    }
  43    std::pair<double, double> operator()(double const& x)
  44    { // Return both f(x) and f'(x).
  45       double fx = x*x*x - a; // Difference (estimate x^3 - value).
  46       double dx = 3 * x*x; // 1st derivative = 3x^2.
  47       return std::make_pair(fx, dx); // 'return' both fx and dx.
  48    }
  49 private:
  50    double a; // to be 'cube_rooted'.
  51 };
  52 // Using 1st and 2nd derivatives with Halley algorithm.
  53 struct cbrt_functor_2deriv
  54 { // Functor returning both 1st and 2nd derivatives.
  55    cbrt_functor_2deriv(double const& to_find_root_of) : a(to_find_root_of)
  56    { // Constructor stores value a to find root of, for example:
  57       // calling cbrt_functor_2deriv<double>(x) to get cube root of x,
  58    }
  59    std::tuple<double, double, double> operator()(double const& x)
  60    { // Return both f(x) and f'(x) and f''(x).
  61       double fx = x*x*x - a; // Difference (estimate x^3 - value).
  62       double dx = 3 * x*x; // 1st derivative = 3x^2.
  63       double d2x = 6 * x; // 2nd derivative = 6x.
  64       return std::make_tuple(fx, dx, d2x); // 'return' fx, dx and d2x.
  65    }
  66 private:
  67    double a; // to be 'cube_rooted'.
  68 };
  69
  70 BOOST_AUTO_TEST_CASE( test_main )
  71 {
  72    int newton_limits = static_cast<int>(std::numeric_limits<double>::digits * 0.6);
  73
  74    double arg = 1e-50;
  75    while(arg < 1e50)
  76    {
  77       double result = boost::math::cbrt(arg);
  78       //
  79       // Start with a really bad guess 5 times below the result:
  80       //
  81       double guess = result / 5;
  82       boost::uintmax_t iters = 1000;
  83       // TOMS algo first:
  84       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);
  85       BOOST_CHECK_CLOSE_FRACTION((r.first + r.second) / 2, result, std::numeric_limits<double>::epsilon() * 4);
  86       BOOST_CHECK_LE(iters, 14);
  87       // Newton next:
  88       iters = 1000;
  89       double dr = boost::math::tools::newton_raphson_iterate(cbrt_functor_deriv(arg), guess, guess / 2, result * 10, newton_limits, iters);
  90       BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
  91       BOOST_CHECK_LE(iters, 12);
  92       // Halley next:
  93       iters = 1000;
  94       dr = boost::math::tools::halley_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters);
  95       BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
  96       BOOST_CHECK_LE(iters, 7);
  97       // Schroder next:
  98       iters = 1000;
  99       dr = boost::math::tools::schroder_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters);
 100       BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
 101       BOOST_CHECK_LE(iters, 11);
 102       //
 103       // Over again with a bad guess 5 times larger than the result:
 104       //
 105       iters = 1000;
 106       guess = result * 5;
 107       r = boost::math::tools::bracket_and_solve_root(cbrt_functor_noderiv(arg), guess, 2.0, true, boost::math::tools::eps_tolerance<double>(), iters);
 108       BOOST_CHECK_CLOSE_FRACTION((r.first + r.second) / 2, result, std::numeric_limits<double>::epsilon() * 4);
 109       BOOST_CHECK_LE(iters, 14);
 110       // Newton next:
 111       iters = 1000;
 112       dr = boost::math::tools::newton_raphson_iterate(cbrt_functor_deriv(arg), guess, result / 10, result * 10, newton_limits, iters);
 113       BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
 114       BOOST_CHECK_LE(iters, 12);
 115       // Halley next:
 116       iters = 1000;
 117       dr = boost::math::tools::halley_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters);
 118       BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
 119       BOOST_CHECK_LE(iters, 7);
 120       // Schroder next:
 121       iters = 1000;
 122       dr = boost::math::tools::schroder_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters);
 123       BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
 124       BOOST_CHECK_LE(iters, 11);
 125       //
 126       // A much better guess, 1% below result:
 127       //
 128       iters = 1000;
 129       guess = result * 0.9;
 130       r = boost::math::tools::bracket_and_solve_root(cbrt_functor_noderiv(arg), guess, 2.0, true, boost::math::tools::eps_tolerance<double>(), iters);
 131       BOOST_CHECK_CLOSE_FRACTION((r.first + r.second) / 2, result, std::numeric_limits<double>::epsilon() * 4);
 132       BOOST_CHECK_LE(iters, 12);
 133       // Newton next:
 134       iters = 1000;
 135       dr = boost::math::tools::newton_raphson_iterate(cbrt_functor_deriv(arg), guess, result / 10, result * 10, newton_limits, iters);
 136       BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
 137       BOOST_CHECK_LE(iters, 5);
 138       // Halley next:
 139       iters = 1000;
 140       dr = boost::math::tools::halley_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters);
 141       BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
 142       BOOST_CHECK_LE(iters, 3);
 143       // Schroder next:
 144       iters = 1000;
 145       dr = boost::math::tools::schroder_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters);
 146       BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
 147       BOOST_CHECK_LE(iters, 4);
 148       //
 149       // A much better guess, 1% above result:
 150       //
 151       iters = 1000;
 152       guess = result * 1.1;
 153       r = boost::math::tools::bracket_and_solve_root(cbrt_functor_noderiv(arg), guess, 2.0, true, boost::math::tools::eps_tolerance<double>(), iters);
 154       BOOST_CHECK_CLOSE_FRACTION((r.first + r.second) / 2, result, std::numeric_limits<double>::epsilon() * 4);
 155       BOOST_CHECK_LE(iters, 12);
 156       // Newton next:
 157       iters = 1000;
 158       dr = boost::math::tools::newton_raphson_iterate(cbrt_functor_deriv(arg), guess, result / 10, result * 10, newton_limits, iters);
 159       BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
 160       BOOST_CHECK_LE(iters, 5);
 161       // Halley next:
 162       iters = 1000;
 163       dr = boost::math::tools::halley_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters);
 164       BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
 165       BOOST_CHECK_LE(iters, 3);
 166       // Schroder next:
 167       iters = 1000;
 168       dr = boost::math::tools::schroder_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters);
 169       BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
 170       BOOST_CHECK_LE(iters, 4);
 171
 172       arg *= 3.5;
 173    }
 174 }
 175
 176 #else
 177
 178 int main() { return 0; }
 179
 180 #endif