[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"

#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;
   std::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);
#if defined(__PPC__) || defined(__aarch64__)
         BOOST_CHECK_LE(iters, 55);
#else
         BOOST_CHECK_LE(iters, 40);
#endif
      }
      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);
      }
   }

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