special_functions/detail/igamma_inverse.hpp

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