third_party/boost/include/boost/math/special_functions/detail/ibeta_inv_ab.hpp - webm/webmlive - Git at Google

 //  (C) Copyright John Maddock 2006.
 //  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)

 //
 // This is not a complete header file, it is included by beta.hpp
 // after it has defined it's definitions.  This inverts the incomplete
 // beta functions ibeta and ibetac on the first parameters "a"
 // and "b" using a generic root finding algorithm (TOMS Algorithm 748).
 //

 #ifndef BOOST_MATH_SP_DETAIL_BETA_INV_AB
 #define BOOST_MATH_SP_DETAIL_BETA_INV_AB

 #ifdef _MSC_VER
 #pragma once
 #endif

 #include <boost/math/tools/toms748_solve.hpp>
 #include <boost/cstdint.hpp>

 namespace boost{ namespace math{ namespace detail{

 template <class T, class Policy>
 struct beta_inv_ab_t
 {
    beta_inv_ab_t(T b_, T z_, T p_, bool invert_, bool swap_ab_) : b(b_), z(z_), p(p_), invert(invert_), swap_ab(swap_ab_) {}
    T operator()(T a)
    {
       return invert ?
          p - boost::math::ibetac(swap_ab ? b : a, swap_ab ? a : b, z, Policy())
          : boost::math::ibeta(swap_ab ? b : a, swap_ab ? a : b, z, Policy()) - p;
    }
 private:
    T b, z, p;
    bool invert, swap_ab;
 };

 template <class T, class Policy>
 T inverse_negative_binomial_cornish_fisher(T n, T sf, T sfc, T p, T q, const Policy& pol)
 {
    BOOST_MATH_STD_USING
    // mean:
    T m = n * (sfc) / sf;
    T t = sqrt(n * (sfc));
    // standard deviation:
    T sigma = t / sf;
    // skewness
    T sk = (1 + sfc) / t;
    // kurtosis:
    T k = (6 - sf * (5+sfc)) / (n * (sfc));
    // Get the inverse of a std normal distribution:
    T x = boost::math::erfc_inv(p > q ? 2 * q : 2 * p, pol) * constants::root_two<T>();
    // Set the sign:
    if(p < 0.5)
       x = -x;
    T x2 = x * x;
    // w is correction term due to skewness
    T w = x + sk * (x2 - 1) / 6;
    //
    // Add on correction due to kurtosis.
    //
    if(n >= 10)
       w += k * x * (x2 - 3) / 24 + sk * sk * x * (2 * x2 - 5) / -36;

    w = m + sigma * w;
    if(w < tools::min_value<T>())
       return tools::min_value<T>();
    return w;
 }

 template <class T, class Policy>
 T ibeta_inv_ab_imp(const T& b, const T& z, const T& p, const T& q, bool swap_ab, const Policy& pol)
 {
    BOOST_MATH_STD_USING  // for ADL of std lib math functions
    //
    // Special cases first:
    //
    BOOST_MATH_INSTRUMENT_CODE("b = " << b << " z = " << z << " p = " << p << " q = " << " swap = " << swap_ab);
    if(p == 0)
    {
       return swap_ab ? tools::min_value<T>() : tools::max_value<T>();
    }
    if(q == 0)
    {
       return swap_ab ? tools::max_value<T>() : tools::min_value<T>();
    }
    //
    // Function object, this is the functor whose root
    // we have to solve:
    //
    beta_inv_ab_t<T, Policy> f(b, z, (p < q) ? p : q, (p < q) ? false : true, swap_ab);
    //
    // Tolerance: full precision.
    //
    tools::eps_tolerance<T> tol(policies::digits<T, Policy>());
    //
    // Now figure out a starting guess for what a may be,
    // we'll start out with a value that'll put p or q
    // right bang in the middle of their range, the functions
    // are quite sensitive so we should need too many steps
    // to bracket the root from there:
    //
    T guess = 0;
    T factor = 5;
    //
    // Convert variables to parameters of a negative binomial distribution:
    //
    T n = b;
    T sf = swap_ab ? z : 1-z;
    T sfc = swap_ab ? 1-z : z;
    T u = swap_ab ? p : q;
    T v = swap_ab ? q : p;
    if(u <= pow(sf, n))
    {
       //
       // Result is less than 1, negative binomial approximation
       // is useless....
       //
       if((p < q) != swap_ab)
       {
          guess = (std::min)(T(b * 2), T(1));
       }
       else
       {
          guess = (std::min)(T(b / 2), T(1));
       }
    }
    if(n * n * n * u * sf > 0.005)
       guess = 1 + inverse_negative_binomial_cornish_fisher(n, sf, sfc, u, v, pol);

    if(guess < 10)
    {
       //
       // Negative binomial approximation not accurate in this area:
       //
       if((p < q) != swap_ab)
       {
          guess = (std::min)(T(b * 2), T(10));
       }
       else
       {
          guess = (std::min)(T(b / 2), T(10));
       }
    }
    else
       factor = (v < sqrt(tools::epsilon<T>())) ? 2 : (guess < 20 ? 1.2f : 1.1f);
    BOOST_MATH_INSTRUMENT_CODE("guess = " << guess);
    //
    // Max iterations permitted:
    //
    boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
    std::pair<T, T> r = bracket_and_solve_root(f, guess, factor, swap_ab ? true : false, tol, max_iter, pol);
    if(max_iter >= policies::get_max_root_iterations<Policy>())
       policies::raise_evaluation_error<T>("boost::math::ibeta_invab_imp<%1%>(%1%,%1%,%1%)", "Unable to locate the root within a reasonable number of iterations, closest approximation so far was %1%", r.first, pol);
    return (r.first + r.second) / 2;
 }

 } // namespace detail

 template <class RT1, class RT2, class RT3, class Policy>
 typename tools::promote_args<RT1, RT2, RT3>::type
       ibeta_inva(RT1 b, RT2 x, RT3 p, const Policy& pol)
 {
    typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
    typedef typename policies::evaluation<result_type, Policy>::type value_type;
    typedef typename policies::normalise<
       Policy,
       policies::promote_float<false>,
       policies::promote_double<false>,
       policies::discrete_quantile<>,
       policies::assert_undefined<> >::type forwarding_policy;

    if(p == 0)
    {
       return tools::max_value<result_type>();
    }
    if(p == 1)
    {
       return tools::min_value<result_type>();
    }

    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
       detail::ibeta_inv_ab_imp(
          static_cast<value_type>(b),
          static_cast<value_type>(x),
          static_cast<value_type>(p),
          static_cast<value_type>(1 - static_cast<value_type>(p)),
          false, pol),
       "boost::math::ibeta_inva<%1%>(%1%,%1%,%1%)");
 }

 template <class RT1, class RT2, class RT3, class Policy>
 typename tools::promote_args<RT1, RT2, RT3>::type
       ibetac_inva(RT1 b, RT2 x, RT3 q, const Policy& pol)
 {
    typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
    typedef typename policies::evaluation<result_type, Policy>::type value_type;
    typedef typename policies::normalise<
       Policy,
       policies::promote_float<false>,
       policies::promote_double<false>,
       policies::discrete_quantile<>,
       policies::assert_undefined<> >::type forwarding_policy;

    if(q == 1)
    {
       return tools::max_value<result_type>();
    }
    if(q == 0)
    {
       return tools::min_value<result_type>();
    }

    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
       detail::ibeta_inv_ab_imp(
          static_cast<value_type>(b),
          static_cast<value_type>(x),
          static_cast<value_type>(1 - static_cast<value_type>(q)),
          static_cast<value_type>(q),
          false, pol),
       "boost::math::ibetac_inva<%1%>(%1%,%1%,%1%)");
 }

 template <class RT1, class RT2, class RT3, class Policy>
 typename tools::promote_args<RT1, RT2, RT3>::type
       ibeta_invb(RT1 a, RT2 x, RT3 p, const Policy& pol)
 {
    typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
    typedef typename policies::evaluation<result_type, Policy>::type value_type;
    typedef typename policies::normalise<
       Policy,
       policies::promote_float<false>,
       policies::promote_double<false>,
       policies::discrete_quantile<>,
       policies::assert_undefined<> >::type forwarding_policy;

    if(p == 0)
    {
       return tools::min_value<result_type>();
    }
    if(p == 1)
    {
       return tools::max_value<result_type>();
    }

    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
       detail::ibeta_inv_ab_imp(
          static_cast<value_type>(a),
          static_cast<value_type>(x),
          static_cast<value_type>(p),
          static_cast<value_type>(1 - static_cast<value_type>(p)),
          true, pol),
       "boost::math::ibeta_invb<%1%>(%1%,%1%,%1%)");
 }

 template <class RT1, class RT2, class RT3, class Policy>
 typename tools::promote_args<RT1, RT2, RT3>::type
       ibetac_invb(RT1 a, RT2 x, RT3 q, const Policy& pol)
 {
    typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
    typedef typename policies::evaluation<result_type, Policy>::type value_type;
    typedef typename policies::normalise<
       Policy,
       policies::promote_float<false>,
       policies::promote_double<false>,
       policies::discrete_quantile<>,
       policies::assert_undefined<> >::type forwarding_policy;

    if(q == 1)
    {
       return tools::min_value<result_type>();
    }
    if(q == 0)
    {
       return tools::max_value<result_type>();
    }

    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
       detail::ibeta_inv_ab_imp(
          static_cast<value_type>(a),
          static_cast<value_type>(x),
          static_cast<value_type>(1 - static_cast<value_type>(q)),
          static_cast<value_type>(q),
          true, pol),
          "boost::math::ibetac_invb<%1%>(%1%,%1%,%1%)");
 }

 template <class RT1, class RT2, class RT3>
 inline typename tools::promote_args<RT1, RT2, RT3>::type
          ibeta_inva(RT1 b, RT2 x, RT3 p)
 {
    return boost::math::ibeta_inva(b, x, p, policies::policy<>());
 }

 template <class RT1, class RT2, class RT3>
 inline typename tools::promote_args<RT1, RT2, RT3>::type
          ibetac_inva(RT1 b, RT2 x, RT3 q)
 {
    return boost::math::ibetac_inva(b, x, q, policies::policy<>());
 }

 template <class RT1, class RT2, class RT3>
 inline typename tools::promote_args<RT1, RT2, RT3>::type
          ibeta_invb(RT1 a, RT2 x, RT3 p)
 {
    return boost::math::ibeta_invb(a, x, p, policies::policy<>());
 }

 template <class RT1, class RT2, class RT3>
 inline typename tools::promote_args<RT1, RT2, RT3>::type
          ibetac_invb(RT1 a, RT2 x, RT3 q)
 {
    return boost::math::ibetac_invb(a, x, q, policies::policy<>());
 }

 } // namespace math
 } // namespace boost

 #endif // BOOST_MATH_SP_DETAIL_BETA_INV_AB
	// (C) Copyright John Maddock 2006.
	// 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)

	//
	// This is not a complete header file, it is included by beta.hpp
	// after it has defined it's definitions. This inverts the incomplete
	// beta functions ibeta and ibetac on the first parameters "a"
	// and "b" using a generic root finding algorithm (TOMS Algorithm 748).
	//

	#ifndef BOOST_MATH_SP_DETAIL_BETA_INV_AB
	#define BOOST_MATH_SP_DETAIL_BETA_INV_AB

	#ifdef _MSC_VER
	#pragma once
	#endif

	#include <boost/math/tools/toms748_solve.hpp>
	#include <boost/cstdint.hpp>

	namespace boost{ namespace math{ namespace detail{

	template <class T, class Policy>
	struct beta_inv_ab_t
	{
	beta_inv_ab_t(T b_, T z_, T p_, bool invert_, bool swap_ab_) : b(b_), z(z_), p(p_), invert(invert_), swap_ab(swap_ab_) {}
	T operator()(T a)
	{
	return invert ?
	p - boost::math::ibetac(swap_ab ? b : a, swap_ab ? a : b, z, Policy())
	: boost::math::ibeta(swap_ab ? b : a, swap_ab ? a : b, z, Policy()) - p;
	}
	private:
	T b, z, p;
	bool invert, swap_ab;
	};

	template <class T, class Policy>
	T inverse_negative_binomial_cornish_fisher(T n, T sf, T sfc, T p, T q, const Policy& pol)
	{
	BOOST_MATH_STD_USING
	// mean:
	T m = n * (sfc) / sf;
	T t = sqrt(n * (sfc));
	// standard deviation:
	T sigma = t / sf;
	// skewness
	T sk = (1 + sfc) / t;
	// kurtosis:
	T k = (6 - sf * (5+sfc)) / (n * (sfc));
	// Get the inverse of a std normal distribution:
	T x = boost::math::erfc_inv(p > q ? 2 * q : 2 * p, pol) * constants::root_two<T>();
	// Set the sign:
	if(p < 0.5)
	x = -x;
	T x2 = x * x;
	// w is correction term due to skewness
	T w = x + sk * (x2 - 1) / 6;
	//
	// Add on correction due to kurtosis.
	//
	if(n >= 10)
	w += k * x * (x2 - 3) / 24 + sk * sk * x * (2 * x2 - 5) / -36;

	w = m + sigma * w;
	if(w < tools::min_value<T>())
	return tools::min_value<T>();
	return w;
	}

	template <class T, class Policy>
	T ibeta_inv_ab_imp(const T& b, const T& z, const T& p, const T& q, bool swap_ab, const Policy& pol)
	{
	BOOST_MATH_STD_USING // for ADL of std lib math functions
	//
	// Special cases first:
	//
	BOOST_MATH_INSTRUMENT_CODE("b = " << b << " z = " << z << " p = " << p << " q = " << " swap = " << swap_ab);
	if(p == 0)
	{
	return swap_ab ? tools::min_value<T>() : tools::max_value<T>();
	}
	if(q == 0)
	{
	return swap_ab ? tools::max_value<T>() : tools::min_value<T>();
	}
	//
	// Function object, this is the functor whose root
	// we have to solve:
	//
	beta_inv_ab_t<T, Policy> f(b, z, (p < q) ? p : q, (p < q) ? false : true, swap_ab);
	//
	// Tolerance: full precision.
	//
	tools::eps_tolerance<T> tol(policies::digits<T, Policy>());
	//
	// Now figure out a starting guess for what a may be,
	// we'll start out with a value that'll put p or q
	// right bang in the middle of their range, the functions
	// are quite sensitive so we should need too many steps
	// to bracket the root from there:
	//
	T guess = 0;
	T factor = 5;
	//
	// Convert variables to parameters of a negative binomial distribution:
	//
	T n = b;
	T sf = swap_ab ? z : 1-z;
	T sfc = swap_ab ? 1-z : z;
	T u = swap_ab ? p : q;
	T v = swap_ab ? q : p;
	if(u <= pow(sf, n))
	{
	//
	// Result is less than 1, negative binomial approximation
	// is useless....
	//
	if((p < q) != swap_ab)
	{
	guess = (std::min)(T(b * 2), T(1));
	}
	else
	{
	guess = (std::min)(T(b / 2), T(1));
	}
	}
	if(n * n * n * u * sf > 0.005)
	guess = 1 + inverse_negative_binomial_cornish_fisher(n, sf, sfc, u, v, pol);

	if(guess < 10)
	{
	//
	// Negative binomial approximation not accurate in this area:
	//
	if((p < q) != swap_ab)
	{
	guess = (std::min)(T(b * 2), T(10));
	}
	else
	{
	guess = (std::min)(T(b / 2), T(10));
	}
	}
	else
	factor = (v < sqrt(tools::epsilon<T>())) ? 2 : (guess < 20 ? 1.2f : 1.1f);
	BOOST_MATH_INSTRUMENT_CODE("guess = " << guess);
	//
	// Max iterations permitted:
	//
	boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
	std::pair<T, T> r = bracket_and_solve_root(f, guess, factor, swap_ab ? true : false, tol, max_iter, pol);
	if(max_iter >= policies::get_max_root_iterations<Policy>())
	policies::raise_evaluation_error<T>("boost::math::ibeta_invab_imp<%1%>(%1%,%1%,%1%)", "Unable to locate the root within a reasonable number of iterations, closest approximation so far was %1%", r.first, pol);
	return (r.first + r.second) / 2;
	}

	} // namespace detail

	template <class RT1, class RT2, class RT3, class Policy>
	typename tools::promote_args<RT1, RT2, RT3>::type
	ibeta_inva(RT1 b, RT2 x, RT3 p, const Policy& pol)
	{
	typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
	typedef typename policies::evaluation<result_type, Policy>::type value_type;
	typedef typename policies::normalise<
	Policy,
	policies::promote_float<false>,
	policies::promote_double<false>,
	policies::discrete_quantile<>,
	policies::assert_undefined<> >::type forwarding_policy;

	if(p == 0)
	{
	return tools::max_value<result_type>();
	}
	if(p == 1)
	{
	return tools::min_value<result_type>();
	}

	return policies::checked_narrowing_cast<result_type, forwarding_policy>(
	detail::ibeta_inv_ab_imp(
	static_cast<value_type>(b),
	static_cast<value_type>(x),
	static_cast<value_type>(p),
	static_cast<value_type>(1 - static_cast<value_type>(p)),
	false, pol),
	"boost::math::ibeta_inva<%1%>(%1%,%1%,%1%)");
	}

	template <class RT1, class RT2, class RT3, class Policy>
	typename tools::promote_args<RT1, RT2, RT3>::type
	ibetac_inva(RT1 b, RT2 x, RT3 q, const Policy& pol)
	{
	typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
	typedef typename policies::evaluation<result_type, Policy>::type value_type;
	typedef typename policies::normalise<
	Policy,
	policies::promote_float<false>,
	policies::promote_double<false>,
	policies::discrete_quantile<>,
	policies::assert_undefined<> >::type forwarding_policy;

	if(q == 1)
	{
	return tools::max_value<result_type>();
	}
	if(q == 0)
	{
	return tools::min_value<result_type>();
	}

	return policies::checked_narrowing_cast<result_type, forwarding_policy>(
	detail::ibeta_inv_ab_imp(
	static_cast<value_type>(b),
	static_cast<value_type>(x),
	static_cast<value_type>(1 - static_cast<value_type>(q)),
	static_cast<value_type>(q),
	false, pol),
	"boost::math::ibetac_inva<%1%>(%1%,%1%,%1%)");
	}

	template <class RT1, class RT2, class RT3, class Policy>
	typename tools::promote_args<RT1, RT2, RT3>::type
	ibeta_invb(RT1 a, RT2 x, RT3 p, const Policy& pol)
	{
	typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
	typedef typename policies::evaluation<result_type, Policy>::type value_type;
	typedef typename policies::normalise<
	Policy,
	policies::promote_float<false>,
	policies::promote_double<false>,
	policies::discrete_quantile<>,
	policies::assert_undefined<> >::type forwarding_policy;

	if(p == 0)
	{
	return tools::min_value<result_type>();
	}
	if(p == 1)
	{
	return tools::max_value<result_type>();
	}

	return policies::checked_narrowing_cast<result_type, forwarding_policy>(
	detail::ibeta_inv_ab_imp(
	static_cast<value_type>(a),
	static_cast<value_type>(x),
	static_cast<value_type>(p),
	static_cast<value_type>(1 - static_cast<value_type>(p)),
	true, pol),
	"boost::math::ibeta_invb<%1%>(%1%,%1%,%1%)");
	}

	template <class RT1, class RT2, class RT3, class Policy>
	typename tools::promote_args<RT1, RT2, RT3>::type
	ibetac_invb(RT1 a, RT2 x, RT3 q, const Policy& pol)
	{
	typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
	typedef typename policies::evaluation<result_type, Policy>::type value_type;
	typedef typename policies::normalise<
	Policy,
	policies::promote_float<false>,
	policies::promote_double<false>,
	policies::discrete_quantile<>,
	policies::assert_undefined<> >::type forwarding_policy;

	if(q == 1)
	{
	return tools::min_value<result_type>();
	}
	if(q == 0)
	{
	return tools::max_value<result_type>();
	}

	return policies::checked_narrowing_cast<result_type, forwarding_policy>(
	detail::ibeta_inv_ab_imp(
	static_cast<value_type>(a),
	static_cast<value_type>(x),
	static_cast<value_type>(1 - static_cast<value_type>(q)),
	static_cast<value_type>(q),
	true, pol),
	"boost::math::ibetac_invb<%1%>(%1%,%1%,%1%)");
	}

	template <class RT1, class RT2, class RT3>
	inline typename tools::promote_args<RT1, RT2, RT3>::type
	ibeta_inva(RT1 b, RT2 x, RT3 p)
	{
	return boost::math::ibeta_inva(b, x, p, policies::policy<>());
	}

	template <class RT1, class RT2, class RT3>
	inline typename tools::promote_args<RT1, RT2, RT3>::type
	ibetac_inva(RT1 b, RT2 x, RT3 q)
	{
	return boost::math::ibetac_inva(b, x, q, policies::policy<>());
	}

	template <class RT1, class RT2, class RT3>
	inline typename tools::promote_args<RT1, RT2, RT3>::type
	ibeta_invb(RT1 a, RT2 x, RT3 p)
	{
	return boost::math::ibeta_invb(a, x, p, policies::policy<>());
	}

	template <class RT1, class RT2, class RT3>
	inline typename tools::promote_args<RT1, RT2, RT3>::type
	ibetac_invb(RT1 a, RT2 x, RT3 q)
	{
	return boost::math::ibetac_invb(a, x, q, policies::policy<>());
	}

	} // namespace math
	} // namespace boost

	#endif // BOOST_MATH_SP_DETAIL_BETA_INV_AB