libs/math/test/test_root_finding_concepts.cpp

   1 // Copyright John Maddock 2014
   2
   3 // Use, modification and distribution are subject to the
   4 // Boost Software License, Version 1.0.
   5 // (See accompanying file LICENSE_1_0.txt
   6 // or copy at http://www.boost.org/LICENSE_1_0.txt)
   7
   8 #define BOOST_TEST_MAIN
   9 #include <boost/test/unit_test.hpp> // Boost.Test
  10 #include <boost/test/results_collector.hpp>
  11 #include <boost/test/unit_test.hpp>
  12 #include <boost/test/tools/floating_point_comparison.hpp>
  13 #include <boost/tuple/tuple.hpp>
  14 #include <boost/fusion/include/tuple.hpp>
  15 #include <boost/fusion/include/std_pair.hpp>
  16 #ifndef BOOST_NO_CXX11_HDR_TUPLE
  17 #include <tuple>
  18 #endif
  19
  20 #include <iostream>
  21 #include <iomanip>
  22    using std::cout;
  23    using std::endl;
  24    using std::setprecision;
  25
  26 #include <boost/math/tools/roots.hpp>
  27
  28 //
  29 // We'll use cbrt as an example:
  30 //
  31 struct cbtr_functor_1
  32 {
  33    cbtr_functor_1(double x) : m_target(x) {}
  34    double operator()(double x)
  35    {
  36       return x * x * x - m_target;
  37    }
  38 private:
  39    double m_target;
  40 };
  41
  42 struct cbtr_functor_2a
  43 {
  44    cbtr_functor_2a(double x) : m_target(x) {}
  45    std::pair<double, double> operator()(double x)
  46    {
  47       return std::make_pair(x * x * x - m_target, 3 * x * x);
  48    }
  49 private:
  50    double m_target;
  51 };
  52 #ifndef BOOST_NO_CXX11_HDR_TUPLE
  53 struct cbtr_functor_2b
  54 {
  55    cbtr_functor_2b(double x) : m_target(x) {}
  56    std::tuple<double, double> operator()(double x)
  57    {
  58       return std::tuple<double, double>(x * x * x - m_target, 3 * x * x);
  59    }
  60 private:
  61    double m_target;
  62 };
  63 #endif
  64 struct cbtr_functor_2c
  65 {
  66    cbtr_functor_2c(double x) : m_target(x) {}
  67    boost::tuple<double, double> operator()(double x)
  68    {
  69       return boost::tuple<double, double>(x * x * x - m_target, 3 * x * x);
  70    }
  71 private:
  72    double m_target;
  73 };
  74 struct cbtr_functor_2d
  75 {
  76    cbtr_functor_2d(double x) : m_target(x) {}
  77    boost::fusion::tuple<double, double> operator()(double x)
  78    {
  79       return boost::fusion::tuple<double, double>(x * x * x - m_target, 3 * x * x);
  80    }
  81 private:
  82    double m_target;
  83 };
  84
  85 #ifndef BOOST_NO_CXX11_HDR_TUPLE
  86 struct cbtr_functor_3b
  87 {
  88    cbtr_functor_3b(double x) : m_target(x) {}
  89    std::tuple<double, double, double> operator()(double x)
  90    {
  91       return std::tuple<double, double, double>(x * x * x - m_target, 3 * x * x, 6 * x);
  92    }
  93 private:
  94    double m_target;
  95 };
  96 #endif
  97 struct cbtr_functor_3c
  98 {
  99    cbtr_functor_3c(double x) : m_target(x) {}
 100    boost::tuple<double, double, double> operator()(double x)
 101    {
 102       return boost::tuple<double, double, double>(x * x * x - m_target, 3 * x * x, 6 * x);
 103    }
 104 private:
 105    double m_target;
 106 };
 107 struct cbtr_functor_3d
 108 {
 109    cbtr_functor_3d(double x) : m_target(x) {}
 110    boost::fusion::tuple<double, double, double> operator()(double x)
 111    {
 112       return boost::fusion::tuple<double, double, double>(x * x * x - m_target, 3 * x * x, 6 * x);
 113    }
 114 private:
 115    double m_target;
 116 };
 117
 118
 119 BOOST_AUTO_TEST_CASE( test_main )
 120 {
 121    double x = 27;
 122    double expected = 3;
 123    double result;
 124    double tolerance = std::numeric_limits<double>::epsilon() * 5;
 125    std::pair<double, double> p;
 126    //
 127    // Start by trying the unary functors, bisect first:
 128    //
 129    cbtr_functor_1 f1(x);
 130    boost::math::tools::eps_tolerance<double> t(std::numeric_limits<double>::digits - 1);
 131    p = boost::math::tools::bisect(f1, 0.0, x, t);
 132    result = (p.first + p.second) / 2;
 133    BOOST_CHECK_CLOSE_FRACTION(expected, result, tolerance);
 134    //
 135    // bracket_and_solve_root:
 136    //
 137    boost::uintmax_t max_iter = boost::math::policies::get_max_root_iterations<boost::math::policies::policy<> >();
 138    p = boost::math::tools::bracket_and_solve_root(f1, x, 2.0, true, t, max_iter);
 139    result = (p.first + p.second) / 2;
 140    BOOST_CHECK_CLOSE_FRACTION(expected, result, tolerance);
 141    //
 142    // toms748_solve:
 143    //
 144    max_iter = boost::math::policies::get_max_root_iterations<boost::math::policies::policy<> >();
 145    p = boost::math::tools::toms748_solve(f1, 0.0, x, t, max_iter);
 146    result = (p.first + p.second) / 2;
 147    BOOST_CHECK_CLOSE_FRACTION(expected, result, tolerance);
 148
 149 #ifndef BOOST_NO_CXX11_LAMBDAS
 150    //
 151    // Now try again with C++11 lambda's
 152    //
 153    p = boost::math::tools::bisect([x](double z){ return z * z * z - x; }, 0.0, x, t);
 154    result = (p.first + p.second) / 2;
 155    BOOST_CHECK_CLOSE_FRACTION(expected, result, tolerance);
 156    //
 157    // bracket_and_solve_root:
 158    //
 159    max_iter = boost::math::policies::get_max_root_iterations<boost::math::policies::policy<> >();
 160    p = boost::math::tools::bracket_and_solve_root([x](double z){ return z * z * z - x; }, x, 2.0, true, t, max_iter);
 161    result = (p.first + p.second) / 2;
 162    BOOST_CHECK_CLOSE_FRACTION(expected, result, tolerance);
 163    //
 164    // toms748_solve:
 165    //
 166    max_iter = boost::math::policies::get_max_root_iterations<boost::math::policies::policy<> >();
 167    p = boost::math::tools::toms748_solve([x](double z){ return z * z * z - x; }, 0.0, x, t, max_iter);
 168    result = (p.first + p.second) / 2;
 169    BOOST_CHECK_CLOSE_FRACTION(expected, result, tolerance);
 170 #endif
 171
 172    cbtr_functor_2a f2(x);
 173 #ifndef BOOST_NO_CXX11_HDR_TUPLE
 174    cbtr_functor_2b f3(x);
 175 #endif
 176    cbtr_functor_2c f4(x);
 177    cbtr_functor_2d f5(x);
 178
 179    //
 180    // Binary Functors - newton_raphson_iterate - test each possible tuple type:
 181    //
 182 #ifndef BOOST_NO_CXX11_HDR_TUPLE
 183    result = boost::math::tools::newton_raphson_iterate(f2, x, 0.0, x, std::numeric_limits<double>::digits - 1);
 184    BOOST_CHECK_CLOSE_FRACTION(expected, result, tolerance);
 185    result = boost::math::tools::newton_raphson_iterate(f3, x, 0.0, x, std::numeric_limits<double>::digits - 1);
 186    BOOST_CHECK_CLOSE_FRACTION(expected, result, tolerance);
 187 #endif
 188    result = boost::math::tools::newton_raphson_iterate(f4, x, 0.0, x, std::numeric_limits<double>::digits - 1);
 189    BOOST_CHECK_CLOSE_FRACTION(expected, result, tolerance);
 190    result = boost::math::tools::newton_raphson_iterate(f5, x, 0.0, x, std::numeric_limits<double>::digits - 1);
 191    BOOST_CHECK_CLOSE_FRACTION(expected, result, tolerance);
 192    //
 193    // And again but with lambdas:
 194    //
 195 #ifndef BOOST_NO_CXX11_LAMBDAS
 196    result = boost::math::tools::newton_raphson_iterate([x](double z){ return std::make_pair(z * z * z - x, 3 * z * z); }, x, 0.0, x, std::numeric_limits<double>::digits - 1);
 197    BOOST_CHECK_CLOSE_FRACTION(expected, result, tolerance);
 198 #ifndef BOOST_NO_CXX11_HDR_TUPLE
 199    result = boost::math::tools::newton_raphson_iterate([x](double z){ return std::make_tuple(z * z * z - x, 3 * z * z); }, x, 0.0, x, std::numeric_limits<double>::digits - 1);
 200    BOOST_CHECK_CLOSE_FRACTION(expected, result, tolerance);
 201 #endif
 202    result = boost::math::tools::newton_raphson_iterate([x](double z){ return boost::tuple<double, double>(z * z * z - x, 3 * z * z); }, x, 0.0, x, std::numeric_limits<double>::digits - 1);
 203    BOOST_CHECK_CLOSE_FRACTION(expected, result, tolerance);
 204    result = boost::math::tools::newton_raphson_iterate([x](double z){ return boost::fusion::tuple<double, double>(z * z * z - x, 3 * z * z); }, x, 0.0, x, std::numeric_limits<double>::digits - 1);
 205    BOOST_CHECK_CLOSE_FRACTION(expected, result, tolerance);
 206 #endif
 207
 208 #ifndef BOOST_NO_CXX11_HDR_TUPLE
 209    cbtr_functor_3b f6(x);
 210 #endif
 211    cbtr_functor_3c f7(x);
 212    cbtr_functor_3d f8(x);
 213
 214    //
 215    // Ternary functors:
 216    //
 217 #ifndef BOOST_NO_CXX11_HDR_TUPLE
 218    result = boost::math::tools::halley_iterate(f6, x, 0.0, x, std::numeric_limits<double>::digits - 1);
 219    BOOST_CHECK_CLOSE_FRACTION(expected, result, tolerance);
 220 #endif
 221    result = boost::math::tools::halley_iterate(f7, x, 0.0, x, std::numeric_limits<double>::digits - 1);
 222    BOOST_CHECK_CLOSE_FRACTION(expected, result, tolerance);
 223    result = boost::math::tools::halley_iterate(f8, x, 0.0, x, std::numeric_limits<double>::digits - 1);
 224    BOOST_CHECK_CLOSE_FRACTION(expected, result, tolerance);
 225 #ifndef BOOST_NO_CXX11_LAMBDAS
 226 #ifndef BOOST_NO_CXX11_HDR_TUPLE
 227    result = boost::math::tools::halley_iterate([x](double z){ return std::make_tuple(z * z * z - x, 3 * z * z, 6 * z); }, x, 0.0, x, std::numeric_limits<double>::digits - 1);
 228    BOOST_CHECK_CLOSE_FRACTION(expected, result, tolerance);
 229 #endif
 230    result = boost::math::tools::halley_iterate([x](double z){ return boost::tuple<double, double, double>(z * z * z - x, 3 * z * z, 6 * z); }, x, 0.0, x, std::numeric_limits<double>::digits - 1);
 231    BOOST_CHECK_CLOSE_FRACTION(expected, result, tolerance);
 232    result = boost::math::tools::halley_iterate([x](double z){ return boost::fusion::tuple<double, double, double>(z * z * z - x, 3 * z * z, 6 * z); }, x, 0.0, x, std::numeric_limits<double>::digits - 1);
 233    BOOST_CHECK_CLOSE_FRACTION(expected, result, tolerance);
 234 #endif
 235 #ifndef BOOST_NO_CXX11_HDR_TUPLE
 236    result = boost::math::tools::schroder_iterate(f6, x, 0.0, x, std::numeric_limits<double>::digits - 1);
 237    BOOST_CHECK_CLOSE_FRACTION(expected, result, tolerance);
 238 #endif
 239    result = boost::math::tools::schroder_iterate(f7, x, 0.0, x, std::numeric_limits<double>::digits - 1);
 240    BOOST_CHECK_CLOSE_FRACTION(expected, result, tolerance);
 241    result = boost::math::tools::schroder_iterate(f8, x, 0.0, x, std::numeric_limits<double>::digits - 1);
 242    BOOST_CHECK_CLOSE_FRACTION(expected, result, tolerance);
 243 #ifndef BOOST_NO_CXX11_LAMBDAS
 244 #ifndef BOOST_NO_CXX11_HDR_TUPLE
 245    result = boost::math::tools::schroder_iterate([x](double z){ return std::make_tuple(z * z * z - x, 3 * z * z, 6 * z); }, x, 0.0, x, std::numeric_limits<double>::digits - 1);
 246    BOOST_CHECK_CLOSE_FRACTION(expected, result, tolerance);
 247 #endif
 248    result = boost::math::tools::schroder_iterate([x](double z){ return boost::tuple<double, double, double>(z * z * z - x, 3 * z * z, 6 * z); }, x, 0.0, x, std::numeric_limits<double>::digits - 1);
 249    BOOST_CHECK_CLOSE_FRACTION(expected, result, tolerance);
 250    result = boost::math::tools::schroder_iterate([x](double z){ return boost::fusion::tuple<double, double, double>(z * z * z - x, 3 * z * z, 6 * z); }, x, 0.0, x, std::numeric_limits<double>::digits - 1);
 251    BOOST_CHECK_CLOSE_FRACTION(expected, result, tolerance);
 252 #endif
 253 } // BOOST_AUTO_TEST_CASE( test_main )
 254