special_functions/detail/igamma_inverse.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 #ifndef BOOST_MATH_SPECIAL_FUNCTIONS_IGAMMA_INVERSE_HPP
0008 #define BOOST_MATH_SPECIAL_FUNCTIONS_IGAMMA_INVERSE_HPP
0009
0010 #ifdef _MSC_VER
0011 #pragma once
0012 #endif
0013
0014 #include <boost/math/tools/config.hpp>
0015 #include <boost/math/tools/cstdint.hpp>
0016 #include <boost/math/tools/tuple.hpp>
0017 #include <boost/math/special_functions/gamma.hpp>
0018 #include <boost/math/special_functions/sign.hpp>
0019 #include <boost/math/tools/roots.hpp>
0020 #include <boost/math/policies/error_handling.hpp>
0021
0022 namespace boost{ namespace math{
0023
0024 namespace detail{
0025
0026 template <class T>
0027 BOOST_MATH_GPU_ENABLED T find_inverse_s(T p, T q)
0028 {
0029    //
0030    // Computation of the Incomplete Gamma Function Ratios and their Inverse
0031    // ARMIDO R. DIDONATO and ALFRED H. MORRIS, JR.
0032    // ACM Transactions on Mathematical Software, Vol. 12, No. 4,
0033    // December 1986, Pages 377-393.
0034    //
0035    // See equation 32.
0036    //
0037    BOOST_MATH_STD_USING
0038    T t;
0039    if(p < T(0.5))
0040    {
0041       t = sqrt(-2 * log(p));
0042    }
0043    else
0044    {
0045       t = sqrt(-2 * log(q));
0046    }
0047    BOOST_MATH_STATIC const double a[4] = { 3.31125922108741, 11.6616720288968, 4.28342155967104, 0.213623493715853 };
0048    BOOST_MATH_STATIC const double b[5] = { 1, 6.61053765625462, 6.40691597760039, 1.27364489782223, 0.3611708101884203e-1 };
0049    T s = t - tools::evaluate_polynomial(a, t) / tools::evaluate_polynomial(b, t);
0050    if(p < T(0.5))
0051       s = -s;
0052    return s;
0053 }
0054
0055 template <class T>
0056 BOOST_MATH_GPU_ENABLED T didonato_SN(T a, T x, unsigned N, T tolerance = 0)
0057 {
0058    //
0059    // Computation of the Incomplete Gamma Function Ratios and their Inverse
0060    // ARMIDO R. DIDONATO and ALFRED H. MORRIS, JR.
0061    // ACM Transactions on Mathematical Software, Vol. 12, No. 4,
0062    // December 1986, Pages 377-393.
0063    //
0064    // See equation 34.
0065    //
0066    T sum = 1;
0067    if(N >= 1)
0068    {
0069       T partial = x / (a + 1);
0070       sum += partial;
0071       for(unsigned i = 2; i <= N; ++i)
0072       {
0073          partial *= x / (a + i);
0074          sum += partial;
0075          if(partial < tolerance)
0076             break;
0077       }
0078    }
0079    return sum;
0080 }
0081
0082 template <class T, class Policy>
0083 BOOST_MATH_GPU_ENABLED inline T didonato_FN(T p, T a, T x, unsigned N, T tolerance, const Policy& pol)
0084 {
0085    //
0086    // Computation of the Incomplete Gamma Function Ratios and their Inverse
0087    // ARMIDO R. DIDONATO and ALFRED H. MORRIS, JR.
0088    // ACM Transactions on Mathematical Software, Vol. 12, No. 4,
0089    // December 1986, Pages 377-393.
0090    //
0091    // See equation 34.
0092    //
0093    BOOST_MATH_STD_USING
0094    T u = log(p) + boost::math::lgamma(a + 1, pol);
0095    return exp((u + x - log(didonato_SN(a, x, N, tolerance))) / a);
0096 }
0097
0098 template <class T, class Policy>
0099 BOOST_MATH_GPU_ENABLED T find_inverse_gamma(T a, T p, T q, const Policy& pol, bool* p_has_10_digits)
0100 {
0101    //
0102    // In order to understand what's going on here, you will
0103    // need to refer to:
0104    //
0105    // Computation of the Incomplete Gamma Function Ratios and their Inverse
0106    // ARMIDO R. DIDONATO and ALFRED H. MORRIS, JR.
0107    // ACM Transactions on Mathematical Software, Vol. 12, No. 4,
0108    // December 1986, Pages 377-393.
0109    //
0110    BOOST_MATH_STD_USING
0111
0112    T result;
0113    *p_has_10_digits = false;
0114
0115    if(a == 1)
0116    {
0117       result = -log(q);
0118       BOOST_MATH_INSTRUMENT_VARIABLE(result);
0119    }
0120    else if(a < 1)
0121    {
0122       T g = boost::math::tgamma(a, pol);
0123       T b = q * g;
0124       BOOST_MATH_INSTRUMENT_VARIABLE(g);
0125       BOOST_MATH_INSTRUMENT_VARIABLE(b);
0126       if((b >T(0.6)) || ((b >= T(0.45)) && (a >= T(0.3))))
0127       {
0128          // DiDonato & Morris Eq 21:
0129          //
0130          // There is a slight variation from DiDonato and Morris here:
0131          // the first form given here is unstable when p is close to 1,
0132          // making it impossible to compute the inverse of Q(a,x) for small
0133          // q.  Fortunately the second form works perfectly well in this case.
0134          //
0135          T u;
0136          if((b * q > T(1e-8)) && (q > T(1e-5)))
0137          {
0138             u = pow(p * g * a, 1 / a);
0139             BOOST_MATH_INSTRUMENT_VARIABLE(u);
0140          }
0141          else
0142          {
0143             u = exp((-q / a) - constants::euler<T>());
0144             BOOST_MATH_INSTRUMENT_VARIABLE(u);
0145          }
0146          result = u / (1 - (u / (a + 1)));
0147          BOOST_MATH_INSTRUMENT_VARIABLE(result);
0148       }
0149       else if((a < 0.3) && (b >= 0.35))
0150       {
0151          // DiDonato & Morris Eq 22:
0152          T t = exp(-constants::euler<T>() - b);
0153          T u = t * exp(t);
0154          result = t * exp(u);
0155          BOOST_MATH_INSTRUMENT_VARIABLE(result);
0156       }
0157       else if((b > 0.15) || (a >= 0.3))
0158       {
0159          // DiDonato & Morris Eq 23:
0160          T y = -log(b);
0161          T u = y - (1 - a) * log(y);
0162          result = y - (1 - a) * log(u) - log(1 + (1 - a) / (1 + u));
0163          BOOST_MATH_INSTRUMENT_VARIABLE(result);
0164       }
0165       else if (b > 0.1)
0166       {
0167          // DiDonato & Morris Eq 24:
0168          T y = -log(b);
0169          T u = y - (1 - a) * log(y);
0170          result = y - (1 - a) * log(u) - log((u * u + 2 * (3 - a) * u + (2 - a) * (3 - a)) / (u * u + (5 - a) * u + 2));
0171          BOOST_MATH_INSTRUMENT_VARIABLE(result);
0172       }
0173       else
0174       {
0175          // DiDonato & Morris Eq 25:
0176          T y = -log(b);
0177          T c1 = (a - 1) * log(y);
0178          T c1_2 = c1 * c1;
0179          T c1_3 = c1_2 * c1;
0180          T c1_4 = c1_2 * c1_2;
0181          T a_2 = a * a;
0182          T a_3 = a_2 * a;
0183
0184          T c2 = (a - 1) * (1 + c1);
0185          T c3 = (a - 1) * (-(c1_2 / 2) + (a - 2) * c1 + (3 * a - 5) / 2);
0186          T c4 = (a - 1) * ((c1_3 / 3) - (3 * a - 5) * c1_2 / 2 + (a_2 - 6 * a + 7) * c1 + (11 * a_2 - 46 * a + 47) / 6);
0187          T c5 = (a - 1) * (-(c1_4 / 4)
0188                            + (11 * a - 17) * c1_3 / 6
0189                            + (-3 * a_2 + 13 * a -13) * c1_2
0190                            + (2 * a_3 - 25 * a_2 + 72 * a - 61) * c1 / 2
0191                            + (25 * a_3 - 195 * a_2 + 477 * a - 379) / 12);
0192
0193          T y_2 = y * y;
0194          T y_3 = y_2 * y;
0195          T y_4 = y_2 * y_2;
0196          result = y + c1 + (c2 / y) + (c3 / y_2) + (c4 / y_3) + (c5 / y_4);
0197          BOOST_MATH_INSTRUMENT_VARIABLE(result);
0198          if(b < 1e-28f)
0199             *p_has_10_digits = true;
0200       }
0201    }
0202    else
0203    {
0204       // DiDonato and Morris Eq 31:
0205       T s = find_inverse_s(p, q);
0206
0207       BOOST_MATH_INSTRUMENT_VARIABLE(s);
0208
0209       T s_2 = s * s;
0210       T s_3 = s_2 * s;
0211       T s_4 = s_2 * s_2;
0212       T s_5 = s_4 * s;
0213       T ra = sqrt(a);
0214
0215       BOOST_MATH_INSTRUMENT_VARIABLE(ra);
0216
0217       T w = a + s * ra + (s * s -1) / 3;
0218       w += (s_3 - 7 * s) / (36 * ra);
0219       w -= (3 * s_4 + 7 * s_2 - 16) / (810 * a);
0220       w += (9 * s_5 + 256 * s_3 - 433 * s) / (38880 * a * ra);
0221
0222       BOOST_MATH_INSTRUMENT_VARIABLE(w);
0223
0224       if((a >= 500) && (fabs(1 - w / a) < 1e-6))
0225       {
0226          result = w;
0227          *p_has_10_digits = true;
0228          BOOST_MATH_INSTRUMENT_VARIABLE(result);
0229       }
0230       else if (p > 0.5)
0231       {
0232          if(w < 3 * a)
0233          {
0234             result = w;
0235             BOOST_MATH_INSTRUMENT_VARIABLE(result);
0236          }
0237          else
0238          {
0239             T D = BOOST_MATH_GPU_SAFE_MAX(T(2), T(a * (a - 1)));
0240             T lg = boost::math::lgamma(a, pol);
0241             T lb = log(q) + lg;
0242             if(lb < -D * T(2.3))
0243             {
0244                // DiDonato and Morris Eq 25:
0245                T y = -lb;
0246                T c1 = (a - 1) * log(y);
0247                T c1_2 = c1 * c1;
0248                T c1_3 = c1_2 * c1;
0249                T c1_4 = c1_2 * c1_2;
0250                T a_2 = a * a;
0251                T a_3 = a_2 * a;
0252
0253                T c2 = (a - 1) * (1 + c1);
0254                T c3 = (a - 1) * (-(c1_2 / 2) + (a - 2) * c1 + (3 * a - 5) / 2);
0255                T c4 = (a - 1) * ((c1_3 / 3) - (3 * a - 5) * c1_2 / 2 + (a_2 - 6 * a + 7) * c1 + (11 * a_2 - 46 * a + 47) / 6);
0256                T c5 = (a - 1) * (-(c1_4 / 4)
0257                                  + (11 * a - 17) * c1_3 / 6
0258                                  + (-3 * a_2 + 13 * a -13) * c1_2
0259                                  + (2 * a_3 - 25 * a_2 + 72 * a - 61) * c1 / 2
0260                                  + (25 * a_3 - 195 * a_2 + 477 * a - 379) / 12);
0261
0262                T y_2 = y * y;
0263                T y_3 = y_2 * y;
0264                T y_4 = y_2 * y_2;
0265                result = y + c1 + (c2 / y) + (c3 / y_2) + (c4 / y_3) + (c5 / y_4);
0266                BOOST_MATH_INSTRUMENT_VARIABLE(result);
0267             }
0268             else
0269             {
0270                // DiDonato and Morris Eq 33:
0271                T u = -lb + (a - 1) * log(w) - log(1 + (1 - a) / (1 + w));
0272                result = -lb + (a - 1) * log(u) - log(1 + (1 - a) / (1 + u));
0273                BOOST_MATH_INSTRUMENT_VARIABLE(result);
0274             }
0275          }
0276       }
0277       else
0278       {
0279          T z = w;
0280          T ap1 = a + 1;
0281          T ap2 = a + 2;
0282          if(w < 0.15f * ap1)
0283          {
0284             // DiDonato and Morris Eq 35:
0285             T v = log(p) + boost::math::lgamma(ap1, pol);
0286             z = exp((v + w) / a);
0287             s = boost::math::log1p(z / ap1 * (1 + z / ap2), pol);
0288             z = exp((v + z - s) / a);
0289             s = boost::math::log1p(z / ap1 * (1 + z / ap2), pol);
0290             z = exp((v + z - s) / a);
0291             s = boost::math::log1p(z / ap1 * (1 + z / ap2 * (1 + z / (a + 3))), pol);
0292             z = exp((v + z - s) / a);
0293             BOOST_MATH_INSTRUMENT_VARIABLE(z);
0294          }
0295
0296          if((z <= 0.01 * ap1) || (z > 0.7 * ap1))
0297          {
0298             result = z;
0299             if(z <= T(0.002) * ap1)
0300                *p_has_10_digits = true;
0301             BOOST_MATH_INSTRUMENT_VARIABLE(result);
0302          }
0303          else
0304          {
0305             // DiDonato and Morris Eq 36:
0306             T ls = log(didonato_SN(a, z, 100, T(1e-4)));
0307             T v = log(p) + boost::math::lgamma(ap1, pol);
0308             z = exp((v + z - ls) / a);
0309             result = z * (1 - (a * log(z) - z - v + ls) / (a - z));
0310
0311             BOOST_MATH_INSTRUMENT_VARIABLE(result);
0312          }
0313       }
0314    }
0315    return result;
0316 }
0317
0318 template <class T, class Policy>
0319 struct gamma_p_inverse_func
0320 {
0321    BOOST_MATH_GPU_ENABLED gamma_p_inverse_func(T a_, T p_, bool inv) : a(a_), p(p_), invert(inv)
0322    {
0323       //
0324       // If p is too near 1 then P(x) - p suffers from cancellation
0325       // errors causing our root-finding algorithms to "thrash", better
0326       // to invert in this case and calculate Q(x) - (1-p) instead.
0327       //
0328       // Of course if p is *very* close to 1, then the answer we get will
0329       // be inaccurate anyway (because there's not enough information in p)
0330       // but at least we will converge on the (inaccurate) answer quickly.
0331       //
0332       if(p > T(0.9))
0333       {
0334          p = 1 - p;
0335          invert = !invert;
0336       }
0337    }
0338
0339    BOOST_MATH_GPU_ENABLED boost::math::tuple<T, T, T> operator()(const T& x)const
0340    {
0341       BOOST_FPU_EXCEPTION_GUARD
0342       //
0343       // Calculate P(x) - p and the first two derivates, or if the invert
0344       // flag is set, then Q(x) - q and it's derivatives.
0345       //
0346       typedef typename policies::evaluation<T, Policy>::type value_type;
0347       // typedef typename lanczos::lanczos<T, Policy>::type evaluation_type;
0348       typedef typename policies::normalise<
0349          Policy,
0350          policies::promote_float<false>,
0351          policies::promote_double<false>,
0352          policies::discrete_quantile<>,
0353          policies::assert_undefined<> >::type forwarding_policy;
0354
0355       BOOST_MATH_STD_USING  // For ADL of std functions.
0356
0357       T f, f1;
0358       value_type ft;
0359       f = static_cast<T>(boost::math::detail::gamma_incomplete_imp(
0360                static_cast<value_type>(a),
0361                static_cast<value_type>(x),
0362                true, invert,
0363                forwarding_policy(), &ft));
0364       f1 = static_cast<T>(ft);
0365       T f2;
0366       T div = (a - x - 1) / x;
0367       f2 = f1;
0368       if(fabs(div) > 1)
0369       {
0370          // split if statement to address M1 mac clang bug;
0371          // see issue 826
0372          if (tools::max_value<T>() / fabs(div) < f2)
0373          {
0374             // overflow:
0375             f2 = -tools::max_value<T>() / 2;
0376          }
0377          else
0378          {
0379             f2 *= div;
0380          }
0381       }
0382       else
0383       {
0384          f2 *= div;
0385       }
0386
0387       if(invert)
0388       {
0389          f1 = -f1;
0390          f2 = -f2;
0391       }
0392
0393       return boost::math::make_tuple(static_cast<T>(f - p), f1, f2);
0394    }
0395 private:
0396    T a, p;
0397    bool invert;
0398 };
0399
0400 template <class T, class Policy>
0401 BOOST_MATH_GPU_ENABLED T gamma_p_inv_imp(T a, T p, const Policy& pol)
0402 {
0403    BOOST_MATH_STD_USING  // ADL of std functions.
0404
0405    constexpr auto function = "boost::math::gamma_p_inv<%1%>(%1%, %1%)";
0406
0407    BOOST_MATH_INSTRUMENT_VARIABLE(a);
0408    BOOST_MATH_INSTRUMENT_VARIABLE(p);
0409
0410    if(a <= 0)
0411       return policies::raise_domain_error<T>(function, "Argument a in the incomplete gamma function inverse must be >= 0 (got a=%1%).", a, pol);
0412    if((p < 0) || (p > 1))
0413       return policies::raise_domain_error<T>(function, "Probability must be in the range [0,1] in the incomplete gamma function inverse (got p=%1%).", p, pol);
0414    if(p == 1)
0415       return policies::raise_overflow_error<T>(function, nullptr, Policy());
0416    if(p == 0)
0417       return 0;
0418    bool has_10_digits;
0419    T guess = detail::find_inverse_gamma<T>(a, p, 1 - p, pol, &has_10_digits);
0420    if((policies::digits<T, Policy>() <= 36) && has_10_digits)
0421       return guess;
0422    T lower = tools::min_value<T>();
0423    if(guess <= lower)
0424       guess = tools::min_value<T>();
0425    BOOST_MATH_INSTRUMENT_VARIABLE(guess);
0426    //
0427    // Work out how many digits to converge to, normally this is
0428    // 2/3 of the digits in T, but if the first derivative is very
0429    // large convergence is slow, so we'll bump it up to full
0430    // precision to prevent premature termination of the root-finding routine.
0431    //
0432    unsigned digits = policies::digits<T, Policy>();
0433    if(digits < 30)
0434    {
0435       digits *= 2;
0436       digits /= 3;
0437    }
0438    else
0439    {
0440       digits /= 2;
0441       digits -= 1;
0442    }
0443    if((a < T(0.125)) && (fabs(gamma_p_derivative(a, guess, pol)) > 1 / sqrt(tools::epsilon<T>())))
0444       digits = policies::digits<T, Policy>() - 2;
0445    //
0446    // Go ahead and iterate:
0447    //
0448    boost::math::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
0449
0450    #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0451    guess = tools::halley_iterate(
0452       detail::gamma_p_inverse_func<T, Policy>(a, p, false),
0453       guess,
0454       lower,
0455       tools::max_value<T>(),
0456       digits,
0457       max_iter);
0458    #else
0459    guess = tools::newton_raphson_iterate(
0460       detail::gamma_p_inverse_func<T, Policy>(a, p, false),
0461       guess,
0462       lower,
0463       tools::max_value<T>(),
0464       digits,
0465       max_iter);
0466    #endif
0467
0468    policies::check_root_iterations<T>(function, max_iter, pol);
0469    BOOST_MATH_INSTRUMENT_VARIABLE(guess);
0470    if(guess == lower)
0471       guess = policies::raise_underflow_error<T>(function, "Expected result known to be non-zero, but is smaller than the smallest available number.", pol);
0472    return guess;
0473 }
0474
0475 template <class T, class Policy>
0476 BOOST_MATH_GPU_ENABLED T gamma_q_inv_imp(T a, T q, const Policy& pol)
0477 {
0478    BOOST_MATH_STD_USING  // ADL of std functions.
0479
0480    constexpr auto function = "boost::math::gamma_q_inv<%1%>(%1%, %1%)";
0481
0482    if(a <= 0)
0483       return policies::raise_domain_error<T>(function, "Argument a in the incomplete gamma function inverse must be >= 0 (got a=%1%).", a, pol);
0484    if((q < 0) || (q > 1))
0485       return policies::raise_domain_error<T>(function, "Probability must be in the range [0,1] in the incomplete gamma function inverse (got q=%1%).", q, pol);
0486    if(q == 0)
0487       return policies::raise_overflow_error<T>(function, nullptr, Policy());
0488    if(q == 1)
0489       return 0;
0490    bool has_10_digits;
0491    T guess = detail::find_inverse_gamma<T>(a, 1 - q, q, pol, &has_10_digits);
0492    if((policies::digits<T, Policy>() <= 36) && has_10_digits)
0493       return guess;
0494    T lower = tools::min_value<T>();
0495    if(guess <= lower)
0496       guess = tools::min_value<T>();
0497    //
0498    // Work out how many digits to converge to, normally this is
0499    // 2/3 of the digits in T, but if the first derivative is very
0500    // large convergence is slow, so we'll bump it up to full
0501    // precision to prevent premature termination of the root-finding routine.
0502    //
0503    unsigned digits = policies::digits<T, Policy>();
0504    if(digits < 30)
0505    {
0506       digits *= 2;
0507       digits /= 3;
0508    }
0509    else
0510    {
0511       digits /= 2;
0512       digits -= 1;
0513    }
0514    if((a < 0.125) && (fabs(gamma_p_derivative(a, guess, pol)) > 1 / sqrt(tools::epsilon<T>())))
0515       digits = policies::digits<T, Policy>();
0516    //
0517    // Go ahead and iterate:
0518    //
0519    boost::math::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
0520
0521    #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0522    guess = tools::halley_iterate(
0523       detail::gamma_p_inverse_func<T, Policy>(a, q, true),
0524       guess,
0525       lower,
0526       tools::max_value<T>(),
0527       digits,
0528       max_iter);
0529    #else
0530    guess = tools::newton_raphson_iterate(
0531       detail::gamma_p_inverse_func<T, Policy>(a, q, true),
0532       guess,
0533       lower,
0534       tools::max_value<T>(),
0535       digits,
0536       max_iter);
0537    #endif
0538
0539    policies::check_root_iterations<T>(function, max_iter, pol);
0540    if(guess == lower)
0541       guess = policies::raise_underflow_error<T>(function, "Expected result known to be non-zero, but is smaller than the smallest available number.", pol);
0542    return guess;
0543 }
0544
0545 } // namespace detail
0546
0547 template <class T1, class T2, class Policy>
0548 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T1, T2>::type
0549    gamma_p_inv(T1 a, T2 p, const Policy& pol)
0550 {
0551    typedef typename tools::promote_args<T1, T2>::type result_type;
0552    return detail::gamma_p_inv_imp(
0553       static_cast<result_type>(a),
0554       static_cast<result_type>(p), pol);
0555 }
0556
0557 template <class T1, class T2, class Policy>
0558 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T1, T2>::type
0559    gamma_q_inv(T1 a, T2 p, const Policy& pol)
0560 {
0561    typedef typename tools::promote_args<T1, T2>::type result_type;
0562    return detail::gamma_q_inv_imp(
0563       static_cast<result_type>(a),
0564       static_cast<result_type>(p), pol);
0565 }
0566
0567 template <class T1, class T2>
0568 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T1, T2>::type
0569    gamma_p_inv(T1 a, T2 p)
0570 {
0571    return gamma_p_inv(a, p, policies::policy<>());
0572 }
0573
0574 template <class T1, class T2>
0575 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T1, T2>::type
0576    gamma_q_inv(T1 a, T2 p)
0577 {
0578    return gamma_q_inv(a, p, policies::policy<>());
0579 }
0580
0581 } // namespace math
0582 } // namespace boost
0583
0584 #endif // BOOST_MATH_SPECIAL_FUNCTIONS_IGAMMA_INVERSE_HPP
0585
0586
0587