Chris@16: //  (C) Copyright John Maddock 2006.
Chris@16: //  Use, modification and distribution are subject to the
Chris@16: //  Boost Software License, Version 1.0. (See accompanying file
Chris@16: //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: #ifndef BOOST_MATH_SPECIAL_FUNCTIONS_IGAMMA_INVERSE_HPP
Chris@16: #define BOOST_MATH_SPECIAL_FUNCTIONS_IGAMMA_INVERSE_HPP
Chris@16: 
Chris@16: #ifdef _MSC_VER
Chris@16: #pragma once
Chris@16: #endif
Chris@16: 
Chris@16: #include <boost/math/tools/tuple.hpp>
Chris@16: #include <boost/math/special_functions/gamma.hpp>
Chris@16: #include <boost/math/special_functions/sign.hpp>
Chris@16: #include <boost/math/tools/roots.hpp>
Chris@16: #include <boost/math/policies/error_handling.hpp>
Chris@16: 
Chris@16: namespace boost{ namespace math{
Chris@16: 
Chris@16: namespace detail{
Chris@16: 
Chris@16: template <class T>
Chris@16: T find_inverse_s(T p, T q)
Chris@16: {
Chris@16:    //
Chris@16:    // Computation of the Incomplete Gamma Function Ratios and their Inverse
Chris@16:    // ARMIDO R. DIDONATO and ALFRED H. MORRIS, JR.
Chris@16:    // ACM Transactions on Mathematical Software, Vol. 12, No. 4,
Chris@16:    // December 1986, Pages 377-393.
Chris@16:    //
Chris@16:    // See equation 32.
Chris@16:    //
Chris@16:    BOOST_MATH_STD_USING
Chris@16:    T t;
Chris@16:    if(p < 0.5)
Chris@16:    {
Chris@16:       t = sqrt(-2 * log(p));
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:       t = sqrt(-2 * log(q));
Chris@16:    }
Chris@16:    static const double a[4] = { 3.31125922108741, 11.6616720288968, 4.28342155967104, 0.213623493715853 };
Chris@16:    static const double b[5] = { 1, 6.61053765625462, 6.40691597760039, 1.27364489782223, 0.3611708101884203e-1 };
Chris@16:    T s = t - tools::evaluate_polynomial(a, t) / tools::evaluate_polynomial(b, t);
Chris@16:    if(p < 0.5)
Chris@16:       s = -s;
Chris@16:    return s;
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: T didonato_SN(T a, T x, unsigned N, T tolerance = 0)
Chris@16: {
Chris@16:    //
Chris@16:    // Computation of the Incomplete Gamma Function Ratios and their Inverse
Chris@16:    // ARMIDO R. DIDONATO and ALFRED H. MORRIS, JR.
Chris@16:    // ACM Transactions on Mathematical Software, Vol. 12, No. 4,
Chris@16:    // December 1986, Pages 377-393.
Chris@16:    //
Chris@16:    // See equation 34.
Chris@16:    //
Chris@16:    T sum = 1;
Chris@16:    if(N >= 1)
Chris@16:    {
Chris@16:       T partial = x / (a + 1);
Chris@16:       sum += partial;
Chris@16:       for(unsigned i = 2; i <= N; ++i)
Chris@16:       {
Chris@16:          partial *= x / (a + i);
Chris@16:          sum += partial;
Chris@16:          if(partial < tolerance)
Chris@16:             break;
Chris@16:       }
Chris@16:    }
Chris@16:    return sum;
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline T didonato_FN(T p, T a, T x, unsigned N, T tolerance, const Policy& pol)
Chris@16: {
Chris@16:    //
Chris@16:    // Computation of the Incomplete Gamma Function Ratios and their Inverse
Chris@16:    // ARMIDO R. DIDONATO and ALFRED H. MORRIS, JR.
Chris@16:    // ACM Transactions on Mathematical Software, Vol. 12, No. 4,
Chris@16:    // December 1986, Pages 377-393.
Chris@16:    //
Chris@16:    // See equation 34.
Chris@16:    //
Chris@16:    BOOST_MATH_STD_USING
Chris@16:    T u = log(p) + boost::math::lgamma(a + 1, pol);
Chris@16:    return exp((u + x - log(didonato_SN(a, x, N, tolerance))) / a);
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: T find_inverse_gamma(T a, T p, T q, const Policy& pol, bool* p_has_10_digits)
Chris@16: {
Chris@16:    //
Chris@16:    // In order to understand what's going on here, you will
Chris@16:    // need to refer to:
Chris@16:    //
Chris@16:    // Computation of the Incomplete Gamma Function Ratios and their Inverse
Chris@16:    // ARMIDO R. DIDONATO and ALFRED H. MORRIS, JR.
Chris@16:    // ACM Transactions on Mathematical Software, Vol. 12, No. 4,
Chris@16:    // December 1986, Pages 377-393.
Chris@16:    //
Chris@16:    BOOST_MATH_STD_USING
Chris@16: 
Chris@16:    T result;
Chris@16:    *p_has_10_digits = false;
Chris@16: 
Chris@16:    if(a == 1)
Chris@16:    {
Chris@16:       result = -log(q);
Chris@16:       BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16:    }
Chris@16:    else if(a < 1)
Chris@16:    {
Chris@16:       T g = boost::math::tgamma(a, pol);
Chris@16:       T b = q * g;
Chris@16:       BOOST_MATH_INSTRUMENT_VARIABLE(g);
Chris@16:       BOOST_MATH_INSTRUMENT_VARIABLE(b);
Chris@16:       if((b > 0.6) || ((b >= 0.45) && (a >= 0.3)))
Chris@16:       {
Chris@16:          // DiDonato & Morris Eq 21: 
Chris@16:          //
Chris@16:          // There is a slight variation from DiDonato and Morris here:
Chris@16:          // the first form given here is unstable when p is close to 1,
Chris@16:          // making it impossible to compute the inverse of Q(a,x) for small
Chris@16:          // q.  Fortunately the second form works perfectly well in this case.
Chris@16:          //
Chris@16:          T u;
Chris@16:          if((b * q > 1e-8) && (q > 1e-5))
Chris@16:          {
Chris@16:             u = pow(p * g * a, 1 / a);
Chris@16:             BOOST_MATH_INSTRUMENT_VARIABLE(u);
Chris@16:          }
Chris@16:          else
Chris@16:          {
Chris@16:             u = exp((-q / a) - constants::euler<T>());
Chris@16:             BOOST_MATH_INSTRUMENT_VARIABLE(u);
Chris@16:          }
Chris@16:          result = u / (1 - (u / (a + 1)));
Chris@16:          BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16:       }
Chris@16:       else if((a < 0.3) && (b >= 0.35))
Chris@16:       {
Chris@16:          // DiDonato & Morris Eq 22:
Chris@16:          T t = exp(-constants::euler<T>() - b);
Chris@16:          T u = t * exp(t);
Chris@16:          result = t * exp(u);
Chris@16:          BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16:       }
Chris@16:       else if((b > 0.15) || (a >= 0.3))
Chris@16:       {
Chris@16:          // DiDonato & Morris Eq 23:
Chris@16:          T y = -log(b);
Chris@16:          T u = y - (1 - a) * log(y);
Chris@16:          result = y - (1 - a) * log(u) - log(1 + (1 - a) / (1 + u));
Chris@16:          BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16:       }
Chris@16:       else if (b > 0.1)
Chris@16:       {
Chris@16:          // DiDonato & Morris Eq 24:
Chris@16:          T y = -log(b);
Chris@16:          T u = y - (1 - a) * log(y);
Chris@16:          result = y - (1 - a) * log(u) - log((u * u + 2 * (3 - a) * u + (2 - a) * (3 - a)) / (u * u + (5 - a) * u + 2));
Chris@16:          BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16:       }
Chris@16:       else
Chris@16:       {
Chris@16:          // DiDonato & Morris Eq 25:
Chris@16:          T y = -log(b);
Chris@16:          T c1 = (a - 1) * log(y);
Chris@16:          T c1_2 = c1 * c1;
Chris@16:          T c1_3 = c1_2 * c1;
Chris@16:          T c1_4 = c1_2 * c1_2;
Chris@16:          T a_2 = a * a;
Chris@16:          T a_3 = a_2 * a;
Chris@16: 
Chris@16:          T c2 = (a - 1) * (1 + c1);
Chris@16:          T c3 = (a - 1) * (-(c1_2 / 2) + (a - 2) * c1 + (3 * a - 5) / 2);
Chris@16:          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);
Chris@16:          T c5 = (a - 1) * (-(c1_4 / 4)
Chris@16:                            + (11 * a - 17) * c1_3 / 6
Chris@16:                            + (-3 * a_2 + 13 * a -13) * c1_2
Chris@16:                            + (2 * a_3 - 25 * a_2 + 72 * a - 61) * c1 / 2
Chris@16:                            + (25 * a_3 - 195 * a_2 + 477 * a - 379) / 12);
Chris@16: 
Chris@16:          T y_2 = y * y;
Chris@16:          T y_3 = y_2 * y;
Chris@16:          T y_4 = y_2 * y_2;
Chris@16:          result = y + c1 + (c2 / y) + (c3 / y_2) + (c4 / y_3) + (c5 / y_4);
Chris@16:          BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16:          if(b < 1e-28f)
Chris@16:             *p_has_10_digits = true;
Chris@16:       }
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:       // DiDonato and Morris Eq 31:
Chris@16:       T s = find_inverse_s(p, q);
Chris@16: 
Chris@16:       BOOST_MATH_INSTRUMENT_VARIABLE(s);
Chris@16: 
Chris@16:       T s_2 = s * s;
Chris@16:       T s_3 = s_2 * s;
Chris@16:       T s_4 = s_2 * s_2;
Chris@16:       T s_5 = s_4 * s;
Chris@16:       T ra = sqrt(a);
Chris@16: 
Chris@16:       BOOST_MATH_INSTRUMENT_VARIABLE(ra);
Chris@16: 
Chris@16:       T w = a + s * ra + (s * s -1) / 3;
Chris@16:       w += (s_3 - 7 * s) / (36 * ra);
Chris@16:       w -= (3 * s_4 + 7 * s_2 - 16) / (810 * a);
Chris@16:       w += (9 * s_5 + 256 * s_3 - 433 * s) / (38880 * a * ra);
Chris@16: 
Chris@16:       BOOST_MATH_INSTRUMENT_VARIABLE(w);
Chris@16: 
Chris@16:       if((a >= 500) && (fabs(1 - w / a) < 1e-6))
Chris@16:       {
Chris@16:          result = w;
Chris@16:          *p_has_10_digits = true;
Chris@16:          BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16:       }
Chris@16:       else if (p > 0.5)
Chris@16:       {
Chris@16:          if(w < 3 * a)
Chris@16:          {
Chris@16:             result = w;
Chris@16:             BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16:          }
Chris@16:          else
Chris@16:          {
Chris@16:             T D = (std::max)(T(2), T(a * (a - 1)));
Chris@16:             T lg = boost::math::lgamma(a, pol);
Chris@16:             T lb = log(q) + lg;
Chris@16:             if(lb < -D * 2.3)
Chris@16:             {
Chris@16:                // DiDonato and Morris Eq 25:
Chris@16:                T y = -lb;
Chris@16:                T c1 = (a - 1) * log(y);
Chris@16:                T c1_2 = c1 * c1;
Chris@16:                T c1_3 = c1_2 * c1;
Chris@16:                T c1_4 = c1_2 * c1_2;
Chris@16:                T a_2 = a * a;
Chris@16:                T a_3 = a_2 * a;
Chris@16: 
Chris@16:                T c2 = (a - 1) * (1 + c1);
Chris@16:                T c3 = (a - 1) * (-(c1_2 / 2) + (a - 2) * c1 + (3 * a - 5) / 2);
Chris@16:                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);
Chris@16:                T c5 = (a - 1) * (-(c1_4 / 4)
Chris@16:                                  + (11 * a - 17) * c1_3 / 6
Chris@16:                                  + (-3 * a_2 + 13 * a -13) * c1_2
Chris@16:                                  + (2 * a_3 - 25 * a_2 + 72 * a - 61) * c1 / 2
Chris@16:                                  + (25 * a_3 - 195 * a_2 + 477 * a - 379) / 12);
Chris@16: 
Chris@16:                T y_2 = y * y;
Chris@16:                T y_3 = y_2 * y;
Chris@16:                T y_4 = y_2 * y_2;
Chris@16:                result = y + c1 + (c2 / y) + (c3 / y_2) + (c4 / y_3) + (c5 / y_4);
Chris@16:                BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16:             }
Chris@16:             else
Chris@16:             {
Chris@16:                // DiDonato and Morris Eq 33:
Chris@16:                T u = -lb + (a - 1) * log(w) - log(1 + (1 - a) / (1 + w));
Chris@16:                result = -lb + (a - 1) * log(u) - log(1 + (1 - a) / (1 + u));
Chris@16:                BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16:             }
Chris@16:          }
Chris@16:       }
Chris@16:       else
Chris@16:       {
Chris@16:          T z = w;
Chris@16:          T ap1 = a + 1;
Chris@16:          T ap2 = a + 2;
Chris@16:          if(w < 0.15f * ap1)
Chris@16:          {
Chris@16:             // DiDonato and Morris Eq 35:
Chris@16:             T v = log(p) + boost::math::lgamma(ap1, pol);
Chris@16:             z = exp((v + w) / a);
Chris@101:             s = boost::math::log1p(z / ap1 * (1 + z / ap2), pol);
Chris@16:             z = exp((v + z - s) / a);
Chris@101:             s = boost::math::log1p(z / ap1 * (1 + z / ap2), pol);
Chris@16:             z = exp((v + z - s) / a);
Chris@101:             s = boost::math::log1p(z / ap1 * (1 + z / ap2 * (1 + z / (a + 3))), pol);
Chris@16:             z = exp((v + z - s) / a);
Chris@16:             BOOST_MATH_INSTRUMENT_VARIABLE(z);
Chris@16:          }
Chris@16: 
Chris@16:          if((z <= 0.01 * ap1) || (z > 0.7 * ap1))
Chris@16:          {
Chris@16:             result = z;
Chris@16:             if(z <= 0.002 * ap1)
Chris@16:                *p_has_10_digits = true;
Chris@16:             BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16:          }
Chris@16:          else
Chris@16:          {
Chris@16:             // DiDonato and Morris Eq 36:
Chris@16:             T ls = log(didonato_SN(a, z, 100, T(1e-4)));
Chris@16:             T v = log(p) + boost::math::lgamma(ap1, pol);
Chris@16:             z = exp((v + z - ls) / a);
Chris@16:             result = z * (1 - (a * log(z) - z - v + ls) / (a - z));
Chris@16: 
Chris@16:             BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16:          }
Chris@16:       }
Chris@16:    }
Chris@16:    return result;
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: struct gamma_p_inverse_func
Chris@16: {
Chris@16:    gamma_p_inverse_func(T a_, T p_, bool inv) : a(a_), p(p_), invert(inv)
Chris@16:    {
Chris@16:       //
Chris@16:       // If p is too near 1 then P(x) - p suffers from cancellation
Chris@16:       // errors causing our root-finding algorithms to "thrash", better
Chris@16:       // to invert in this case and calculate Q(x) - (1-p) instead.
Chris@16:       //
Chris@16:       // Of course if p is *very* close to 1, then the answer we get will
Chris@16:       // be inaccurate anyway (because there's not enough information in p)
Chris@16:       // but at least we will converge on the (inaccurate) answer quickly.
Chris@16:       //
Chris@16:       if(p > 0.9)
Chris@16:       {
Chris@16:          p = 1 - p;
Chris@16:          invert = !invert;
Chris@16:       }
Chris@16:    }
Chris@16: 
Chris@16:    boost::math::tuple<T, T, T> operator()(const T& x)const
Chris@16:    {
Chris@16:       BOOST_FPU_EXCEPTION_GUARD
Chris@16:       //
Chris@16:       // Calculate P(x) - p and the first two derivates, or if the invert
Chris@16:       // flag is set, then Q(x) - q and it's derivatives.
Chris@16:       //
Chris@16:       typedef typename policies::evaluation<T, Policy>::type value_type;
Chris@16:       // typedef typename lanczos::lanczos<T, Policy>::type evaluation_type;
Chris@16:       typedef typename policies::normalise<
Chris@16:          Policy, 
Chris@16:          policies::promote_float<false>, 
Chris@16:          policies::promote_double<false>, 
Chris@16:          policies::discrete_quantile<>,
Chris@16:          policies::assert_undefined<> >::type forwarding_policy;
Chris@16: 
Chris@16:       BOOST_MATH_STD_USING  // For ADL of std functions.
Chris@16: 
Chris@16:       T f, f1;
Chris@16:       value_type ft;
Chris@16:       f = static_cast<T>(boost::math::detail::gamma_incomplete_imp(
Chris@16:                static_cast<value_type>(a), 
Chris@16:                static_cast<value_type>(x), 
Chris@16:                true, invert,
Chris@16:                forwarding_policy(), &ft));
Chris@16:       f1 = static_cast<T>(ft);
Chris@16:       T f2;
Chris@16:       T div = (a - x - 1) / x;
Chris@16:       f2 = f1;
Chris@16:       if((fabs(div) > 1) && (tools::max_value<T>() / fabs(div) < f2))
Chris@16:       {
Chris@16:          // overflow:
Chris@16:          f2 = -tools::max_value<T>() / 2;
Chris@16:       }
Chris@16:       else
Chris@16:       {
Chris@16:          f2 *= div;
Chris@16:       }
Chris@16: 
Chris@16:       if(invert)
Chris@16:       {
Chris@16:          f1 = -f1;
Chris@16:          f2 = -f2;
Chris@16:       }
Chris@16: 
Chris@16:       return boost::math::make_tuple(static_cast<T>(f - p), f1, f2);
Chris@16:    }
Chris@16: private:
Chris@16:    T a, p;
Chris@16:    bool invert;
Chris@16: };
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: T gamma_p_inv_imp(T a, T p, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING  // ADL of std functions.
Chris@16: 
Chris@16:    static const char* function = "boost::math::gamma_p_inv<%1%>(%1%, %1%)";
Chris@16: 
Chris@16:    BOOST_MATH_INSTRUMENT_VARIABLE(a);
Chris@16:    BOOST_MATH_INSTRUMENT_VARIABLE(p);
Chris@16: 
Chris@16:    if(a <= 0)
Chris@101:       return policies::raise_domain_error<T>(function, "Argument a in the incomplete gamma function inverse must be >= 0 (got a=%1%).", a, pol);
Chris@16:    if((p < 0) || (p > 1))
Chris@101:       return policies::raise_domain_error<T>(function, "Probabilty must be in the range [0,1] in the incomplete gamma function inverse (got p=%1%).", p, pol);
Chris@16:    if(p == 1)
Chris@101:       return policies::raise_overflow_error<T>(function, 0, Policy());
Chris@16:    if(p == 0)
Chris@16:       return 0;
Chris@16:    bool has_10_digits;
Chris@16:    T guess = detail::find_inverse_gamma<T>(a, p, 1 - p, pol, &has_10_digits);
Chris@16:    if((policies::digits<T, Policy>() <= 36) && has_10_digits)
Chris@16:       return guess;
Chris@16:    T lower = tools::min_value<T>();
Chris@16:    if(guess <= lower)
Chris@16:       guess = tools::min_value<T>();
Chris@16:    BOOST_MATH_INSTRUMENT_VARIABLE(guess);
Chris@16:    //
Chris@16:    // Work out how many digits to converge to, normally this is
Chris@16:    // 2/3 of the digits in T, but if the first derivative is very
Chris@16:    // large convergence is slow, so we'll bump it up to full 
Chris@16:    // precision to prevent premature termination of the root-finding routine.
Chris@16:    //
Chris@16:    unsigned digits = policies::digits<T, Policy>();
Chris@16:    if(digits < 30)
Chris@16:    {
Chris@16:       digits *= 2;
Chris@16:       digits /= 3;
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:       digits /= 2;
Chris@16:       digits -= 1;
Chris@16:    }
Chris@16:    if((a < 0.125) && (fabs(gamma_p_derivative(a, guess, pol)) > 1 / sqrt(tools::epsilon<T>())))
Chris@16:       digits = policies::digits<T, Policy>() - 2;
Chris@16:    //
Chris@16:    // Go ahead and iterate:
Chris@16:    //
Chris@16:    boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
Chris@16:    guess = tools::halley_iterate(
Chris@16:       detail::gamma_p_inverse_func<T, Policy>(a, p, false),
Chris@16:       guess,
Chris@16:       lower,
Chris@16:       tools::max_value<T>(),
Chris@16:       digits,
Chris@16:       max_iter);
Chris@16:    policies::check_root_iterations<T>(function, max_iter, pol);
Chris@16:    BOOST_MATH_INSTRUMENT_VARIABLE(guess);
Chris@16:    if(guess == lower)
Chris@16:       guess = policies::raise_underflow_error<T>(function, "Expected result known to be non-zero, but is smaller than the smallest available number.", pol);
Chris@16:    return guess;
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: T gamma_q_inv_imp(T a, T q, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING  // ADL of std functions.
Chris@16: 
Chris@16:    static const char* function = "boost::math::gamma_q_inv<%1%>(%1%, %1%)";
Chris@16: 
Chris@16:    if(a <= 0)
Chris@101:       return policies::raise_domain_error<T>(function, "Argument a in the incomplete gamma function inverse must be >= 0 (got a=%1%).", a, pol);
Chris@16:    if((q < 0) || (q > 1))
Chris@101:       return policies::raise_domain_error<T>(function, "Probabilty must be in the range [0,1] in the incomplete gamma function inverse (got q=%1%).", q, pol);
Chris@16:    if(q == 0)
Chris@101:       return policies::raise_overflow_error<T>(function, 0, Policy());
Chris@16:    if(q == 1)
Chris@16:       return 0;
Chris@16:    bool has_10_digits;
Chris@16:    T guess = detail::find_inverse_gamma<T>(a, 1 - q, q, pol, &has_10_digits);
Chris@16:    if((policies::digits<T, Policy>() <= 36) && has_10_digits)
Chris@16:       return guess;
Chris@16:    T lower = tools::min_value<T>();
Chris@16:    if(guess <= lower)
Chris@16:       guess = tools::min_value<T>();
Chris@16:    //
Chris@16:    // Work out how many digits to converge to, normally this is
Chris@16:    // 2/3 of the digits in T, but if the first derivative is very
Chris@16:    // large convergence is slow, so we'll bump it up to full 
Chris@16:    // precision to prevent premature termination of the root-finding routine.
Chris@16:    //
Chris@16:    unsigned digits = policies::digits<T, Policy>();
Chris@16:    if(digits < 30)
Chris@16:    {
Chris@16:       digits *= 2;
Chris@16:       digits /= 3;
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:       digits /= 2;
Chris@16:       digits -= 1;
Chris@16:    }
Chris@16:    if((a < 0.125) && (fabs(gamma_p_derivative(a, guess, pol)) > 1 / sqrt(tools::epsilon<T>())))
Chris@16:       digits = policies::digits<T, Policy>();
Chris@16:    //
Chris@16:    // Go ahead and iterate:
Chris@16:    //
Chris@16:    boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
Chris@16:    guess = tools::halley_iterate(
Chris@16:       detail::gamma_p_inverse_func<T, Policy>(a, q, true),
Chris@16:       guess,
Chris@16:       lower,
Chris@16:       tools::max_value<T>(),
Chris@16:       digits,
Chris@16:       max_iter);
Chris@16:    policies::check_root_iterations<T>(function, max_iter, pol);
Chris@16:    if(guess == lower)
Chris@16:       guess = policies::raise_underflow_error<T>(function, "Expected result known to be non-zero, but is smaller than the smallest available number.", pol);
Chris@16:    return guess;
Chris@16: }
Chris@16: 
Chris@16: } // namespace detail
Chris@16: 
Chris@16: template <class T1, class T2, class Policy>
Chris@16: inline typename tools::promote_args<T1, T2>::type 
Chris@16:    gamma_p_inv(T1 a, T2 p, const Policy& pol)
Chris@16: {
Chris@16:    typedef typename tools::promote_args<T1, T2>::type result_type;
Chris@16:    return detail::gamma_p_inv_imp(
Chris@16:       static_cast<result_type>(a),
Chris@16:       static_cast<result_type>(p), pol);
Chris@16: }
Chris@16: 
Chris@16: template <class T1, class T2, class Policy>
Chris@16: inline typename tools::promote_args<T1, T2>::type 
Chris@16:    gamma_q_inv(T1 a, T2 p, const Policy& pol)
Chris@16: {
Chris@16:    typedef typename tools::promote_args<T1, T2>::type result_type;
Chris@16:    return detail::gamma_q_inv_imp(
Chris@16:       static_cast<result_type>(a),
Chris@16:       static_cast<result_type>(p), pol);
Chris@16: }
Chris@16: 
Chris@16: template <class T1, class T2>
Chris@16: inline typename tools::promote_args<T1, T2>::type 
Chris@16:    gamma_p_inv(T1 a, T2 p)
Chris@16: {
Chris@16:    return gamma_p_inv(a, p, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: template <class T1, class T2>
Chris@16: inline typename tools::promote_args<T1, T2>::type 
Chris@16:    gamma_q_inv(T1 a, T2 p)
Chris@16: {
Chris@16:    return gamma_q_inv(a, p, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: } // namespace math
Chris@16: } // namespace boost
Chris@16: 
Chris@16: #endif // BOOST_MATH_SPECIAL_FUNCTIONS_IGAMMA_INVERSE_HPP
Chris@16: 
Chris@16: 
Chris@16: