...one of the most highly
regarded and expertly designed C++ library projects in the
world.
— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards
#include <boost/math/special_functions/gamma.hpp>
namespace boost{ namespace math{ template <class T1, class T2> calculatedresulttype gamma_q_inv(T1 a, T2 q); template <class T1, class T2, class Policy> calculatedresulttype gamma_q_inv(T1 a, T2 q, const Policy&); template <class T1, class T2> calculatedresulttype gamma_p_inv(T1 a, T2 p); template <class T1, class T2, class Policy> calculatedresulttype gamma_p_inv(T1 a, T2 p, const Policy&); template <class T1, class T2> calculatedresulttype gamma_q_inva(T1 x, T2 q); template <class T1, class T2, class Policy> calculatedresulttype gamma_q_inva(T1 x, T2 q, const Policy&); template <class T1, class T2> calculatedresulttype gamma_p_inva(T1 x, T2 p); template <class T1, class T2, class Policy> calculatedresulttype gamma_p_inva(T1 x, T2 p, const Policy&); }} // namespaces
There are four incomplete gamma function inverses which either compute x given a and p or q, or else compute a given x and either p or q.
The return type of these functions is computed using the result type calculation rules when T1 and T2 are different types, otherwise the return type is simply T1.
The final Policy argument is optional and can be used to control the behaviour of the function: how it handles errors, what level of precision to use etc. Refer to the policy documentation for more details.
Tip  

When people normally talk about the inverse of the incomplete gamma function, they are talking about inverting on parameter x. These are implemented here as gamma_p_inv and gamma_q_inv, and are by far the most efficient of the inverses presented here. The inverse on the a parameter finds use in some statistical applications but has to be computed by rather brute force numerical techniques and is consequently several times slower. These are implemented here as gamma_p_inva and gamma_q_inva. 
template <class T1, class T2> calculatedresulttype gamma_q_inv(T1 a, T2 q); template <class T1, class T2, class Policy> calculatedresulttype gamma_q_inv(T1 a, T2 q, const Policy&);
Returns a value x such that: q
= gamma_q(a,
x);
Requires: a > 0 and 1 >= p,q >= 0.
template <class T1, class T2> calculatedresulttype gamma_p_inv(T1 a, T2 p); template <class T1, class T2, class Policy> calculatedresulttype gamma_p_inv(T1 a, T2 p, const Policy&);
Returns a value x such that: p
= gamma_p(a,
x);
Requires: a > 0 and 1 >= p,q >= 0.
template <class T1, class T2> calculatedresulttype gamma_q_inva(T1 x, T2 q); template <class T1, class T2, class Policy> calculatedresulttype gamma_q_inva(T1 x, T2 q, const Policy&);
Returns a value a such that: q
= gamma_q(a,
x);
Requires: x > 0 and 1 >= p,q >= 0.
template <class T1, class T2> calculatedresulttype gamma_p_inva(T1 x, T2 p); template <class T1, class T2, class Policy> calculatedresulttype gamma_p_inva(T1 x, T2 p, const Policy&);
Returns a value a such that: p
= gamma_p(a,
x);
Requires: x > 0 and 1 >= p,q >= 0.
The accuracy of these functions doesn't vary much by platform or by the type T. Given that these functions are computed by iterative methods, they are deliberately "detuned" so as not to be too accurate: it is in any case impossible for these function to be more accurate than the regular forward incomplete gamma functions. In practice, the accuracy of these functions is very similar to that of gamma_p and gamma_q functions:
Table 6.13. Error rates for gamma_p_inv
GNU C++ version 7.1.0 
GNU C++ version 7.1.0 
Sun compiler version 0x5150 
Microsoft Visual C++ version 14.1 


incomplete gamma inverse(a, z) medium values 
Max = 0.993ε (Mean = 0.15ε) 
Max = 1.8ε (Mean = 0.406ε) 
Max = 1.89ε (Mean = 0.466ε) 
Max = 1.71ε (Mean = 0.34ε) 
incomplete gamma inverse(a, z) large values 
Max = 0ε (Mean = 0ε) 
Max = 0.509ε (Mean = 0.0447ε) 
Max = 0.509ε (Mean = 0.0447ε) 
Max = 0.924ε (Mean = 0.108ε) 
incomplete gamma inverse(a, z) small values 
Max = 441ε (Mean = 53.9ε) 
Max = 9.17e+03ε (Mean = 1.45e+03ε) 
Max = 1.09e+04ε (Mean = 1.3e+03ε) 
Max = 1.1e+03ε (Mean = 131ε) 
Table 6.14. Error rates for gamma_q_inv
GNU C++ version 7.1.0 
GNU C++ version 7.1.0 
Sun compiler version 0x5150 
Microsoft Visual C++ version 14.1 


incomplete gamma inverse(a, z) medium values 
Max = 0.912ε (Mean = 0.154ε) 
Max = 6.2ε (Mean = 0.627ε) 
Max = 6.2ε (Mean = 0.683ε) 
Max = 2.88ε (Mean = 0.469ε) 
incomplete gamma inverse(a, z) large values 
Max = 0.894ε (Mean = 0.0915ε) 
Max = 0ε (Mean = 0ε) 
Max = 0ε (Mean = 0ε) 
Max = 0.814ε (Mean = 0.0856ε) 
incomplete gamma inverse(a, z) small values 
Max = 292ε (Mean = 36.4ε) 
Max = 8.28e+03ε (Mean = 1.09e+03ε) 
Max = 8.98e+03ε (Mean = 877ε) 
Max = 451ε (Mean = 64.7ε) 
Table 6.15. Error rates for gamma_p_inva
GNU C++ version 7.1.0 
GNU C++ version 7.1.0 
Sun compiler version 0x5150 
Microsoft Visual C++ version 14.1 


Incomplete gamma inverses. 
Max = 0ε (Mean = 0ε) 
Max = 7.87ε (Mean = 1.15ε) 
Max = 4.08ε (Mean = 1.12ε) 
Max = 4.92ε (Mean = 1.03ε) 
Table 6.16. Error rates for gamma_q_inva
GNU C++ version 7.1.0 
GNU C++ version 7.1.0 
Sun compiler version 0x5150 
Microsoft Visual C++ version 14.1 


Incomplete gamma inverses. 
Max = 0ε (Mean = 0ε) 
Max = 8.42ε (Mean = 1.3ε) 
Max = 7.86ε (Mean = 1.24ε) 
Max = 5.05ε (Mean = 1.08ε) 
There are two sets of tests:
The functions gamma_p_inv and gamma_q_inv share a common implementation.
First an initial approximation is computed using the methodology described in:
Finally, the last few bits are cleaned up using Halley iteration, the iteration limit is set to 2/3 of the number of bits in T, which by experiment is sufficient to ensure that the inverses are at least as accurate as the normal incomplete gamma functions. In testing, no more than 3 iterations are required to produce a result as accurate as the forward incomplete gamma function, and in many cases only one iteration is required.
The functions gamma_p_inva and gamma_q_inva also share a common implementation but are handled separately from gamma_p_inv and gamma_q_inv.
An initial approximation for a is computed very crudely so that gamma_p(a, x) ~ 0.5, this value is then used as a starting point for a generic derivativefree root finding algorithm. As a consequence, these two functions are rather more expensive to compute than the gamma_p_inv or gamma_q_inv functions. Even so, the root is usually found in fewer than 10 iterations.