123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325 |
- // (C) Copyright John Maddock 2015.
- // Use, modification and distribution are subject to the
- // Boost Software License, Version 1.0. (See accompanying file
- // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
- #include <pch.hpp>
- #ifndef BOOST_NO_CXX11_HDR_TUPLE
- #define BOOST_TEST_MAIN
- #include <boost/test/unit_test.hpp>
- #include <boost/test/tools/floating_point_comparison.hpp>
- #include <boost/math/tools/roots.hpp>
- #include <boost/test/results_collector.hpp>
- #include <boost/test/unit_test.hpp>
- #include <boost/math/special_functions/cbrt.hpp>
- #include <boost/math/special_functions/beta.hpp>
- #include <iostream>
- #include <iomanip>
- #include <tuple>
- #include "table_type.hpp"
- // No derivatives - using TOMS748 internally.
- struct cbrt_functor_noderiv
- { // cube root of x using only function - no derivatives.
- cbrt_functor_noderiv(double to_find_root_of) : a(to_find_root_of)
- { // Constructor just stores value a to find root of.
- }
- double operator()(double x)
- {
- double fx = x*x*x - a; // Difference (estimate x^3 - a).
- return fx;
- }
- private:
- double a; // to be 'cube_rooted'.
- }; // template <class T> struct cbrt_functor_noderiv
- // Using 1st derivative only Newton-Raphson
- struct cbrt_functor_deriv
- { // Functor also returning 1st derviative.
- cbrt_functor_deriv(double const& to_find_root_of) : a(to_find_root_of)
- { // Constructor stores value a to find root of,
- // for example: calling cbrt_functor_deriv<double>(x) to use to get cube root of x.
- }
- std::pair<double, double> operator()(double const& x)
- { // Return both f(x) and f'(x).
- double fx = x*x*x - a; // Difference (estimate x^3 - value).
- double dx = 3 * x*x; // 1st derivative = 3x^2.
- return std::make_pair(fx, dx); // 'return' both fx and dx.
- }
- private:
- double a; // to be 'cube_rooted'.
- };
- // Using 1st and 2nd derivatives with Halley algorithm.
- struct cbrt_functor_2deriv
- { // Functor returning both 1st and 2nd derivatives.
- cbrt_functor_2deriv(double const& to_find_root_of) : a(to_find_root_of)
- { // Constructor stores value a to find root of, for example:
- // calling cbrt_functor_2deriv<double>(x) to get cube root of x,
- }
- std::tuple<double, double, double> operator()(double const& x)
- { // Return both f(x) and f'(x) and f''(x).
- double fx = x*x*x - a; // Difference (estimate x^3 - value).
- double dx = 3 * x*x; // 1st derivative = 3x^2.
- double d2x = 6 * x; // 2nd derivative = 6x.
- return std::make_tuple(fx, dx, d2x); // 'return' fx, dx and d2x.
- }
- private:
- double a; // to be 'cube_rooted'.
- };
- template <class T, class Policy>
- struct ibeta_roots_1 // for first order algorithms
- {
- ibeta_roots_1(T _a, T _b, T t, bool inv = false)
- : a(_a), b(_b), target(t), invert(inv) {}
- T operator()(const T& x)
- {
- return boost::math::detail::ibeta_imp(a, b, x, Policy(), invert, true) - target;
- }
- private:
- T a, b, target;
- bool invert;
- };
- template <class T, class Policy>
- struct ibeta_roots_2 // for second order algorithms
- {
- ibeta_roots_2(T _a, T _b, T t, bool inv = false)
- : a(_a), b(_b), target(t), invert(inv) {}
- boost::math::tuple<T, T> operator()(const T& x)
- {
- typedef boost::math::lanczos::lanczos<T, Policy> S;
- typedef typename S::type L;
- T f = boost::math::detail::ibeta_imp(a, b, x, Policy(), invert, true) - target;
- T f1 = invert ?
- -boost::math::detail::ibeta_power_terms(b, a, 1 - x, x, L(), true, Policy())
- : boost::math::detail::ibeta_power_terms(a, b, x, 1 - x, L(), true, Policy());
- T y = 1 - x;
- if (y == 0)
- y = boost::math::tools::min_value<T>() * 8;
- f1 /= y * x;
- // make sure we don't have a zero derivative:
- if (f1 == 0)
- f1 = (invert ? -1 : 1) * boost::math::tools::min_value<T>() * 64;
- return boost::math::make_tuple(f, f1);
- }
- private:
- T a, b, target;
- bool invert;
- };
- template <class T, class Policy>
- struct ibeta_roots_3 // for third order algorithms
- {
- ibeta_roots_3(T _a, T _b, T t, bool inv = false)
- : a(_a), b(_b), target(t), invert(inv) {}
- boost::math::tuple<T, T, T> operator()(const T& x)
- {
- typedef typename boost::math::lanczos::lanczos<T, Policy>::type L;
- T f = boost::math::detail::ibeta_imp(a, b, x, Policy(), invert, true) - target;
- T f1 = invert ?
- -boost::math::detail::ibeta_power_terms(b, a, 1 - x, x, L(), true, Policy())
- : boost::math::detail::ibeta_power_terms(a, b, x, 1 - x, L(), true, Policy());
- T y = 1 - x;
- if (y == 0)
- y = boost::math::tools::min_value<T>() * 8;
- f1 /= y * x;
- T f2 = f1 * (-y * a + (b - 2) * x + 1) / (y * x);
- if (invert)
- f2 = -f2;
- // make sure we don't have a zero derivative:
- if (f1 == 0)
- f1 = (invert ? -1 : 1) * boost::math::tools::min_value<T>() * 64;
- return boost::math::make_tuple(f, f1, f2);
- }
- private:
- T a, b, target;
- bool invert;
- };
- BOOST_AUTO_TEST_CASE( test_main )
- {
- int newton_limits = static_cast<int>(std::numeric_limits<double>::digits * 0.6);
- double arg = 1e-50;
- boost::uintmax_t iters;
- double guess;
- double dr;
- while(arg < 1e50)
- {
- double result = boost::math::cbrt(arg);
- //
- // Start with a really bad guess 5 times below the result:
- //
- guess = result / 5;
- iters = 1000;
- // TOMS algo first:
- 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);
- BOOST_CHECK_CLOSE_FRACTION((r.first + r.second) / 2, result, std::numeric_limits<double>::epsilon() * 4);
- BOOST_CHECK_LE(iters, 14);
- // Newton next:
- iters = 1000;
- dr = boost::math::tools::newton_raphson_iterate(cbrt_functor_deriv(arg), guess, guess / 2, result * 10, newton_limits, iters);
- BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
- BOOST_CHECK_LE(iters, 12);
- // Halley next:
- iters = 1000;
- dr = boost::math::tools::halley_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters);
- BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
- BOOST_CHECK_LE(iters, 7);
- // Schroder next:
- iters = 1000;
- dr = boost::math::tools::schroder_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters);
- BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
- BOOST_CHECK_LE(iters, 11);
- //
- // Over again with a bad guess 5 times larger than the result:
- //
- iters = 1000;
- guess = result * 5;
- r = boost::math::tools::bracket_and_solve_root(cbrt_functor_noderiv(arg), guess, 2.0, true, boost::math::tools::eps_tolerance<double>(), iters);
- BOOST_CHECK_CLOSE_FRACTION((r.first + r.second) / 2, result, std::numeric_limits<double>::epsilon() * 4);
- BOOST_CHECK_LE(iters, 14);
- // Newton next:
- iters = 1000;
- dr = boost::math::tools::newton_raphson_iterate(cbrt_functor_deriv(arg), guess, result / 10, result * 10, newton_limits, iters);
- BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
- BOOST_CHECK_LE(iters, 12);
- // Halley next:
- iters = 1000;
- dr = boost::math::tools::halley_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters);
- BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
- BOOST_CHECK_LE(iters, 7);
- // Schroder next:
- iters = 1000;
- dr = boost::math::tools::schroder_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters);
- BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
- BOOST_CHECK_LE(iters, 11);
- //
- // A much better guess, 1% below result:
- //
- iters = 1000;
- guess = result * 0.9;
- r = boost::math::tools::bracket_and_solve_root(cbrt_functor_noderiv(arg), guess, 2.0, true, boost::math::tools::eps_tolerance<double>(), iters);
- BOOST_CHECK_CLOSE_FRACTION((r.first + r.second) / 2, result, std::numeric_limits<double>::epsilon() * 4);
- BOOST_CHECK_LE(iters, 12);
- // Newton next:
- iters = 1000;
- dr = boost::math::tools::newton_raphson_iterate(cbrt_functor_deriv(arg), guess, result / 10, result * 10, newton_limits, iters);
- BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
- BOOST_CHECK_LE(iters, 5);
- // Halley next:
- iters = 1000;
- dr = boost::math::tools::halley_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters);
- BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
- BOOST_CHECK_LE(iters, 3);
- // Schroder next:
- iters = 1000;
- dr = boost::math::tools::schroder_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters);
- BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
- BOOST_CHECK_LE(iters, 4);
- //
- // A much better guess, 1% above result:
- //
- iters = 1000;
- guess = result * 1.1;
- r = boost::math::tools::bracket_and_solve_root(cbrt_functor_noderiv(arg), guess, 2.0, true, boost::math::tools::eps_tolerance<double>(), iters);
- BOOST_CHECK_CLOSE_FRACTION((r.first + r.second) / 2, result, std::numeric_limits<double>::epsilon() * 4);
- BOOST_CHECK_LE(iters, 12);
- // Newton next:
- iters = 1000;
- dr = boost::math::tools::newton_raphson_iterate(cbrt_functor_deriv(arg), guess, result / 10, result * 10, newton_limits, iters);
- BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
- BOOST_CHECK_LE(iters, 5);
- // Halley next:
- iters = 1000;
- dr = boost::math::tools::halley_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters);
- BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
- BOOST_CHECK_LE(iters, 3);
- // Schroder next:
- iters = 1000;
- dr = boost::math::tools::schroder_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters);
- BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2);
- BOOST_CHECK_LE(iters, 4);
- arg *= 3.5;
- }
- //
- // Test ibeta as this triggers all the pathological cases!
- //
- #ifndef SC_
- #define SC_(x) x
- #endif
- #define T double
- # include "ibeta_small_data.ipp"
- for (unsigned i = 0; i < ibeta_small_data.size(); ++i)
- {
- //
- // These inverse tests are thrown off if the output of the
- // incomplete beta is too close to 1: basically there is insuffient
- // information left in the value we're using as input to the inverse
- // to be able to get back to the original value.
- //
- if (ibeta_small_data[i][5] == 0)
- {
- iters = 1000;
- 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);
- BOOST_CHECK_EQUAL(dr, 0.0);
- BOOST_CHECK_LE(iters, 27);
- iters = 1000;
- 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);
- BOOST_CHECK_EQUAL(dr, 0.0);
- BOOST_CHECK_LE(iters, 10);
- }
- else if ((1 - ibeta_small_data[i][5] > 0.001)
- && (fabs(ibeta_small_data[i][5]) > 2 * boost::math::tools::min_value<double>()))
- {
- iters = 1000;
- double result = ibeta_small_data[i][2];
- 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);
- BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 200);
- #if defined(BOOST_MSVC) && (BOOST_MSVC == 1600)
- BOOST_CHECK_LE(iters, 40);
- #else
- BOOST_CHECK_LE(iters, 27);
- #endif
- iters = 1000;
- result = ibeta_small_data[i][2];
- 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);
- BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 200);
- BOOST_CHECK_LE(iters, 40);
- }
- else if (1 == ibeta_small_data[i][5])
- {
- iters = 1000;
- 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);
- BOOST_CHECK_EQUAL(dr, 1.0);
- BOOST_CHECK_LE(iters, 27);
- iters = 1000;
- 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);
- BOOST_CHECK_EQUAL(dr, 1.0);
- BOOST_CHECK_LE(iters, 10);
- }
- }
- }
- #else
- int main() { return 0; }
- #endif
|