special_functions/detail/gamma_inva.hpp

0001 //  (C) Copyright John Maddock 2006.
0002 //  (C) Copyright Matt Borland 2024.
0003 //  Use, modification and distribution are subject to the
0004 //  Boost Software License, Version 1.0. (See accompanying file
0005 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0006
0007 //
0008 // This is not a complete header file, it is included by gamma.hpp
0009 // after it has defined it's definitions.  This inverts the incomplete
0010 // gamma functions P and Q on the first parameter "a" using a generic
0011 // root finding algorithm (TOMS Algorithm 748).
0012 //
0013
0014 #ifndef BOOST_MATH_SP_DETAIL_GAMMA_INVA
0015 #define BOOST_MATH_SP_DETAIL_GAMMA_INVA
0016
0017 #ifdef _MSC_VER
0018 #pragma once
0019 #endif
0020
0021 #include <boost/math/tools/config.hpp>
0022 #include <boost/math/tools/toms748_solve.hpp>
0023
0024 namespace boost{ namespace math{
0025
0026 #ifdef BOOST_MATH_HAS_NVRTC
0027 template <typename T, typename Policy>
0028 BOOST_MATH_GPU_ENABLED auto erfc_inv(T x, const Policy&);
0029 #endif
0030
0031 namespace detail{
0032
0033 template <class T, class Policy>
0034 struct gamma_inva_t
0035 {
0036    BOOST_MATH_GPU_ENABLED gamma_inva_t(T z_, T p_, bool invert_) : z(z_), p(p_), invert(invert_) {}
0037    BOOST_MATH_GPU_ENABLED T operator()(T a)
0038    {
0039       return invert ? p - boost::math::gamma_q(a, z, Policy()) : boost::math::gamma_p(a, z, Policy()) - p;
0040    }
0041 private:
0042    T z, p;
0043    bool invert;
0044 };
0045
0046 template <class T, class Policy>
0047 BOOST_MATH_GPU_ENABLED T inverse_poisson_cornish_fisher(T lambda, T p, T q, const Policy& pol)
0048 {
0049    BOOST_MATH_STD_USING
0050    // mean:
0051    T m = lambda;
0052    // standard deviation:
0053    T sigma = sqrt(lambda);
0054    // skewness
0055    T sk = 1 / sigma;
0056    // kurtosis:
0057    // T k = 1/lambda;
0058    // Get the inverse of a std normal distribution:
0059    T x = boost::math::erfc_inv(p > q ? 2 * q : 2 * p, pol) * constants::root_two<T>();
0060    // Set the sign:
0061    if(p < 0.5)
0062       x = -x;
0063    T x2 = x * x;
0064    // w is correction term due to skewness
0065    T w = x + sk * (x2 - 1) / 6;
0066    /*
0067    // Add on correction due to kurtosis.
0068    // Disabled for now, seems to make things worse?
0069    //
0070    if(lambda >= 10)
0071       w += k * x * (x2 - 3) / 24 + sk * sk * x * (2 * x2 - 5) / -36;
0072    */
0073    w = m + sigma * w;
0074    return w > tools::min_value<T>() ? w : tools::min_value<T>();
0075 }
0076
0077 template <class T, class Policy>
0078 BOOST_MATH_GPU_ENABLED T gamma_inva_imp(const T& z, const T& p, const T& q, const Policy& pol)
0079 {
0080    BOOST_MATH_STD_USING  // for ADL of std lib math functions
0081    //
0082    // Special cases first:
0083    //
0084    if(p == 0)
0085    {
0086       return policies::raise_overflow_error<T>("boost::math::gamma_p_inva<%1%>(%1%, %1%)", nullptr, Policy());
0087    }
0088    if(q == 0)
0089    {
0090       return tools::min_value<T>();
0091    }
0092    //
0093    // Function object, this is the functor whose root
0094    // we have to solve:
0095    //
0096    gamma_inva_t<T, Policy> f(z, (p < q) ? p : q, (p < q) ? false : true);
0097    //
0098    // Tolerance: full precision.
0099    //
0100    tools::eps_tolerance<T> tol(policies::digits<T, Policy>());
0101    //
0102    // Now figure out a starting guess for what a may be,
0103    // we'll start out with a value that'll put p or q
0104    // right bang in the middle of their range, the functions
0105    // are quite sensitive so we should need too many steps
0106    // to bracket the root from there:
0107    //
0108    T guess;
0109    T factor = 8;
0110    if(z >= 1)
0111    {
0112       //
0113       // We can use the relationship between the incomplete
0114       // gamma function and the poisson distribution to
0115       // calculate an approximate inverse, for large z
0116       // this is actually pretty accurate, but it fails badly
0117       // when z is very small.  Also set our step-factor according
0118       // to how accurate we think the result is likely to be:
0119       //
0120       guess = 1 + inverse_poisson_cornish_fisher(z, q, p, pol);
0121       if(z > 5)
0122       {
0123          if(z > 1000)
0124             factor = 1.01f;
0125          else if(z > 50)
0126             factor = 1.1f;
0127          else if(guess > 10)
0128             factor = 1.25f;
0129          else
0130             factor = 2;
0131          if(guess < 1.1)
0132             factor = 8;
0133       }
0134    }
0135    else if(z > 0.5)
0136    {
0137       guess = z * 1.2f;
0138    }
0139    else
0140    {
0141       guess = -0.4f / log(z);
0142    }
0143    //
0144    // Max iterations permitted:
0145    //
0146    std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
0147    //
0148    // Use our generic derivative-free root finding procedure.
0149    // We could use Newton steps here, taking the PDF of the
0150    // Poisson distribution as our derivative, but that's
0151    // even worse performance-wise than the generic method :-(
0152    //
0153    std::pair<T, T> r = bracket_and_solve_root(f, guess, factor, false, tol, max_iter, pol);
0154    if(max_iter >= policies::get_max_root_iterations<Policy>())
0155       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);
0156    return (r.first + r.second) / 2;
0157 }
0158
0159 } // namespace detail
0160
0161 template <class T1, class T2, class Policy>
0162 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T1, T2>::type
0163    gamma_p_inva(T1 x, T2 p, const Policy& pol)
0164 {
0165    typedef typename tools::promote_args<T1, T2>::type result_type;
0166    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0167    typedef typename policies::normalise<
0168       Policy,
0169       policies::promote_float<false>,
0170       policies::promote_double<false>,
0171       policies::discrete_quantile<>,
0172       policies::assert_undefined<> >::type forwarding_policy;
0173
0174    if(p == 0)
0175    {
0176       policies::raise_overflow_error<result_type>("boost::math::gamma_p_inva<%1%>(%1%, %1%)", nullptr, Policy());
0177    }
0178    if(p == 1)
0179    {
0180       return tools::min_value<result_type>();
0181    }
0182
0183    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
0184       detail::gamma_inva_imp(
0185          static_cast<value_type>(x),
0186          static_cast<value_type>(p),
0187          static_cast<value_type>(1 - static_cast<value_type>(p)),
0188          pol), "boost::math::gamma_p_inva<%1%>(%1%, %1%)");
0189 }
0190
0191 template <class T1, class T2, class Policy>
0192 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T1, T2>::type
0193    gamma_q_inva(T1 x, T2 q, const Policy& pol)
0194 {
0195    typedef typename tools::promote_args<T1, T2>::type result_type;
0196    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0197    typedef typename policies::normalise<
0198       Policy,
0199       policies::promote_float<false>,
0200       policies::promote_double<false>,
0201       policies::discrete_quantile<>,
0202       policies::assert_undefined<> >::type forwarding_policy;
0203
0204    if(q == 1)
0205    {
0206       policies::raise_overflow_error<result_type>("boost::math::gamma_q_inva<%1%>(%1%, %1%)", nullptr, Policy());
0207    }
0208    if(q == 0)
0209    {
0210       return tools::min_value<result_type>();
0211    }
0212
0213    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
0214       detail::gamma_inva_imp(
0215          static_cast<value_type>(x),
0216          static_cast<value_type>(1 - static_cast<value_type>(q)),
0217          static_cast<value_type>(q),
0218          pol), "boost::math::gamma_q_inva<%1%>(%1%, %1%)");
0219 }
0220
0221 template <class T1, class T2>
0222 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T1, T2>::type
0223    gamma_p_inva(T1 x, T2 p)
0224 {
0225    return boost::math::gamma_p_inva(x, p, policies::policy<>());
0226 }
0227
0228 template <class T1, class T2>
0229 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T1, T2>::type
0230    gamma_q_inva(T1 x, T2 q)
0231 {
0232    return boost::math::gamma_q_inva(x, q, policies::policy<>());
0233 }
0234
0235 } // namespace math
0236 } // namespace boost
0237
0238 #endif // BOOST_MATH_SP_DETAIL_GAMMA_INVA
0239
0240
0241