[ceph.git] / ceph / src / boost / libs / math / test / test_root_iterations.cpp

//  (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 derivative.
   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
Commit	Line	Data
7c673cae FG	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>
92f5a8d4	12	#include <boost/test/tools/floating_point_comparison.hpp>
7c673cae FG	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>
92f5a8d4	17	#include <boost/math/special_functions/beta.hpp>
7c673cae FG	18	#include <iostream>
	19	#include <iomanip>
	20	#include <tuple>
92f5a8d4	21	#include "table_type.hpp"
7c673cae FG	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 = xxx - 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
f67539c2	40	{ // Functor also returning 1st derivative.
7c673cae FG	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 = xxx - 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 = xxx - 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
92f5a8d4 TL	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
7c673cae FG	150	BOOST_AUTO_TEST_CASE( test_main )
	151	{
	152	int newton_limits = static_cast<int>(std::numeric_limits<double>::digits * 0.6);
b32b8144	153
7c673cae	154	double arg = 1e-50;
92f5a8d4 TL	155	boost::uintmax_t iters;
	156	double guess;
	157	double dr;
	158
7c673cae FG	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	//
92f5a8d4 TL	165	guess = result / 5;
92f5a8d4 TL	166	iters = 1000;
7c673cae FG	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;
92f5a8d4	173	dr = boost::math::tools::newton_raphson_iterate(cbrt_functor_deriv(arg), guess, guess / 2, result * 10, newton_limits, iters);
7c673cae FG	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	}
92f5a8d4 TL	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
7c673cae FG	319	}
	320
	321	#else
	322
	323	int main() { return 0; }
	324
	325	#endif