GAMMA.INV silently returned the first bisection midpoint for alpha in
~[143, 171.62] (e.g. GAMMA.INV(0.5, 143, 1) gave 358.0 instead of
142.667) and #NUM! above that, because the Newton-step pdf evaluates
Gamma(a), b**a and value**(a-1) in linear domain, all of which overflow
even though the density itself is a small representable number. The same
pattern breaks the GAMMA.DIST, CHISQ.DIST and F.DIST densities and
GAMMALN, which computed log(Gamma(x)) through Gamma(x).
Evaluate these in log domain via the existing logGamma, and scale the
incomplete-gamma series/continued-fraction iteration cap as O(sqrt(a)),
which both expansions need to converge near x ~ a once the shape is
large (GAMMA.INV drifted from the true quantile above alpha ~5000 and
returned alpha+1 by alpha=10000; same for CHISQ.INV at high df).
While we might never be able to have 100% of our code strict, we can at
the very least do it for all of our tests. This ensures that our tests
are using our API with the types as intended by the test author, and not
silently be cast to what our API requires.
Vincent Gao
Adrien Crivelli
Mark Baker