incomplete_gamma_inv : Unexpected NaNs for specific values

[tdegeus]

I tried incomplete_gamma_inv to help me convert a flat random distribution to a gamma-distribution. When I try a large sample most inputs work well, but I get NaNs with some specific values as listed below:

#include <gcem.hpp>
#include <iostream>
#include <vector>

int main()
{
    std::vector<double> r = {2.26397533e-05, 9.99999672e-01, 7.41391908e-04, 3.88840912e-04,
                             7.33291963e-04, 2.62747984e-04, 2.54816143e-04, 1.69432024e-04,
                             1.75788999e-04, 2.59617111e-04, 4.94140433e-04, 6.44463347e-04};

    for (auto& i : r) {
        std::cout << gcem::incomplete_gamma_inv(2.0, i) << std::endl;
    }

    return 0;
}

which gives:

nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan

Note that most are quite small, but not even all of them. Do you know what is going on here?

[tdegeus]

Note that instead, for example scipy does compute (presumably) correct values:

from scipy.special import gammaincinv
import numpy as np

r = np.array([2.26397533e-05, 9.99999672e-01, 7.41391908e-04, 3.88840912e-04, 
              7.33291963e-04, 2.62747984e-04, 2.54816143e-04, 1.69432024e-04, 
              1.75788999e-04, 2.59617111e-04, 4.94140433e-04, 6.44463347e-04])

for i in r:
    print(gammaincinv(2.0, i))

which gives

0.006744144756417914
17.867703900499794
0.03901009547759127
0.028149536712293647
0.03879362103603407
0.0231007268352203
0.022746692866731894
0.01852217725342957
0.01886862303518837
0.022961618953584943
0.03177118744354654
0.0363384721865337

Likewise Boost does the job:

#include <iostream>
#include <vector>
#include <boost/math/special_functions/gamma.hpp>

int main()
{
    std::vector<double> r = {2.26397533e-05, 9.99999672e-01, 7.41391908e-04, 3.88840912e-04,
                             7.33291963e-04, 2.62747984e-04, 2.54816143e-04, 1.69432024e-04,
                             1.75788999e-04, 2.59617111e-04, 4.94140433e-04, 6.44463347e-04};

    for (auto& i : r) {
        std::cout << boost::math::gamma_p_inv(2.0, i) << std::endl;
    }

    return 0;
}

as it gives

0.00674414
17.8677
0.0390101
0.0281495
0.0387936
0.0231007
0.0227467
0.0185222
0.0188686
0.0229616
0.0317712
0.0363385