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.

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.

Basic sanity checks attempt to "round-trip" from a
and x to p or q
and back again. These tests have quite generous tolerances: in general
both the incomplete gamma, and its inverses, change so rapidly that
round tripping to more than a couple of significant digits isn't possible.
This is especially true when p or q
is very near one: in this case there isn't enough "information
content" in the input to the inverse function to get back where
you started.

Accuracy checks using high precision test values. These measure the
accuracy of the result, given exact input values.

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 derivative-free 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.