diff options
Diffstat (limited to 'src/third_party/boost-1.69.0/boost/math/special_functions/detail/gamma_inva.hpp')
| -rw-r--r-- | src/third_party/boost-1.69.0/boost/math/special_functions/detail/gamma_inva.hpp | 233 |
1 files changed, 0 insertions, 233 deletions
diff --git a/src/third_party/boost-1.69.0/boost/math/special_functions/detail/gamma_inva.hpp b/src/third_party/boost-1.69.0/boost/math/special_functions/detail/gamma_inva.hpp deleted file mode 100644 index 7c32d2946c0..00000000000 --- a/src/third_party/boost-1.69.0/boost/math/special_functions/detail/gamma_inva.hpp +++ /dev/null @@ -1,233 +0,0 @@ -// (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 gamma.hpp -// after it has defined it's definitions. This inverts the incomplete -// gamma functions P and Q on the first parameter "a" using a generic -// root finding algorithm (TOMS Algorithm 748). -// - -#ifndef BOOST_MATH_SP_DETAIL_GAMMA_INVA -#define BOOST_MATH_SP_DETAIL_GAMMA_INVA - -#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 gamma_inva_t -{ - gamma_inva_t(T z_, T p_, bool invert_) : z(z_), p(p_), invert(invert_) {} - T operator()(T a) - { - return invert ? p - boost::math::gamma_q(a, z, Policy()) : boost::math::gamma_p(a, z, Policy()) - p; - } -private: - T z, p; - bool invert; -}; - -template <class T, class Policy> -T inverse_poisson_cornish_fisher(T lambda, T p, T q, const Policy& pol) -{ - BOOST_MATH_STD_USING - // mean: - T m = lambda; - // standard deviation: - T sigma = sqrt(lambda); - // skewness - T sk = 1 / sigma; - // kurtosis: - // T k = 1/lambda; - // 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. - // Disabled for now, seems to make things worse? - // - if(lambda >= 10) - w += k * x * (x2 - 3) / 24 + sk * sk * x * (2 * x2 - 5) / -36; - */ - w = m + sigma * w; - return w > tools::min_value<T>() ? w : tools::min_value<T>(); -} - -template <class T, class Policy> -T gamma_inva_imp(const T& z, const T& p, const T& q, const Policy& pol) -{ - BOOST_MATH_STD_USING // for ADL of std lib math functions - // - // Special cases first: - // - if(p == 0) - { - return policies::raise_overflow_error<T>("boost::math::gamma_p_inva<%1%>(%1%, %1%)", 0, Policy()); - } - if(q == 0) - { - return tools::min_value<T>(); - } - // - // Function object, this is the functor whose root - // we have to solve: - // - gamma_inva_t<T, Policy> f(z, (p < q) ? p : q, (p < q) ? false : true); - // - // 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; - T factor = 8; - if(z >= 1) - { - // - // We can use the relationship between the incomplete - // gamma function and the poisson distribution to - // calculate an approximate inverse, for large z - // this is actually pretty accurate, but it fails badly - // when z is very small. Also set our step-factor according - // to how accurate we think the result is likely to be: - // - guess = 1 + inverse_poisson_cornish_fisher(z, q, p, pol); - if(z > 5) - { - if(z > 1000) - factor = 1.01f; - else if(z > 50) - factor = 1.1f; - else if(guess > 10) - factor = 1.25f; - else - factor = 2; - if(guess < 1.1) - factor = 8; - } - } - else if(z > 0.5) - { - guess = z * 1.2f; - } - else - { - guess = -0.4f / log(z); - } - // - // Max iterations permitted: - // - boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); - // - // Use our generic derivative-free root finding procedure. - // We could use Newton steps here, taking the PDF of the - // Poisson distribution as our derivative, but that's - // even worse performance-wise than the generic method :-( - // - std::pair<T, T> r = bracket_and_solve_root(f, guess, factor, false, tol, max_iter, pol); - if(max_iter >= policies::get_max_root_iterations<Policy>()) - return policies::raise_evaluation_error<T>("boost::math::gamma_p_inva<%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 T1, class T2, class Policy> -inline typename tools::promote_args<T1, T2>::type - gamma_p_inva(T1 x, T2 p, const Policy& pol) -{ - typedef typename tools::promote_args<T1, T2>::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) - { - policies::raise_overflow_error<result_type>("boost::math::gamma_p_inva<%1%>(%1%, %1%)", 0, Policy()); - } - if(p == 1) - { - return tools::min_value<result_type>(); - } - - return policies::checked_narrowing_cast<result_type, forwarding_policy>( - detail::gamma_inva_imp( - static_cast<value_type>(x), - static_cast<value_type>(p), - static_cast<value_type>(1 - static_cast<value_type>(p)), - pol), "boost::math::gamma_p_inva<%1%>(%1%, %1%)"); -} - -template <class T1, class T2, class Policy> -inline typename tools::promote_args<T1, T2>::type - gamma_q_inva(T1 x, T2 q, const Policy& pol) -{ - typedef typename tools::promote_args<T1, T2>::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) - { - policies::raise_overflow_error<result_type>("boost::math::gamma_q_inva<%1%>(%1%, %1%)", 0, Policy()); - } - if(q == 0) - { - return tools::min_value<result_type>(); - } - - return policies::checked_narrowing_cast<result_type, forwarding_policy>( - detail::gamma_inva_imp( - static_cast<value_type>(x), - static_cast<value_type>(1 - static_cast<value_type>(q)), - static_cast<value_type>(q), - pol), "boost::math::gamma_q_inva<%1%>(%1%, %1%)"); -} - -template <class T1, class T2> -inline typename tools::promote_args<T1, T2>::type - gamma_p_inva(T1 x, T2 p) -{ - return boost::math::gamma_p_inva(x, p, policies::policy<>()); -} - -template <class T1, class T2> -inline typename tools::promote_args<T1, T2>::type - gamma_q_inva(T1 x, T2 q) -{ - return boost::math::gamma_q_inva(x, q, policies::policy<>()); -} - -} // namespace math -} // namespace boost - -#endif // BOOST_MATH_SP_DETAIL_GAMMA_INVA - - - |