Fix gamma_inc for larga a and small x (and very large x)

[andreasnoack]

Fixes #311 as explained in my comment there.

[codecov]

Codecov Report

Merging #314 (ef42fcc) into master (2e90d1b) will increase coverage by 0.35%.
The diff coverage is 93.13%.

@@            Coverage Diff             @@
##           master     #314      +/-   ##
==========================================
+ Coverage   88.44%   88.79%   +0.35%     
==========================================
  Files          12       12              
  Lines        2622     2624       +2     
==========================================
+ Hits         2319     2330      +11     
+ Misses        303      294       -9     

Flag	Coverage Δ
unittests	88.79% <93.13%> (+0.35%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

Impacted Files	Coverage Δ
src/gamma_inc.jl	90.22% <93.13%> (+1.61%)	⬆️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 2e90d1b...ef42fcc. Read the comment docs.

[andreasnoack]

I've added some tests to increase the test coverage. While doing that I noticed a bug in rgammax which I've fixed. I'm now folloing the double precision instead the single precision in NSWC. @devmotion Would you be able to review this one? (It will help us towards getting rid of Rmath in StatsFuns).

Update: I've also renamed the stirling function to stirling_error since it returns the difference between the Stirling approximation of loggamma and loggamma so I think it's confusing to call it stirling although that is what's called in IncgamFI module that the some of the implementations are following.

[devmotion]

I did not check the correctness of the parts where only the formatting was changed. I am sure that rgammax matches the Fortran implementation and that _gamma_inc_choose_algorithm also corresponds to the Fortran code (just had a question about the assert statement). It was a bit tricky to follow the Fortran code for _gamma_inc with all the go-to statements. I couldn't find one of the branches in the Fortran code but otherwise it seems to be correct, in particular the additional branch for _gamma_inc_choose_algorithm if a >= big1[iop] and !(abs(s) <= 0.4).

[devmotion]

In the Fortran code (end of block 20) I can only find a check for abs(s) <= 0.4 which appears after this branch here in SpecialFunctions.

[andreasnoack]

I just looked and it seems to be part of the single precision version GRATIO. It's also block 20 and line 13303. However, we shouldn't be using the single precision version (although it seems that they had larger words back when Morris wrote all this code). I'll update this to follow the double precision routine.

[andreasnoack]

Actually, I think it's better to leave that for a separate PR. The single precision version has the extra IOP argument that controls the precision and we expose that in gamma_inc. As mentioned, the word size was larger on the machines that Morris wrote this for and he writes in the comment to the single precision version GRATIO that

IF IND = 0 THEN THE USER IS REQUESTING AS MUCH ACCURACY AS POSSIBLE (UP TO 14 SIGNIFICANT DIGITS).

so I think it's still okay for Float64. Notice that his double precision literals have 30 digits.

[devmotion]

Yes, this seems reasonable. Thanks for looking into this, I didn't check the single precision version.

[devmotion]

It seems in the Fortran code the cutoff is 2 instead of 1.1 (in the section "SELECT THE APPROPRIATE ALGORITHM" above 10).

[andreasnoack]

Indeed this is similar to the one above. I'll go through the while function.

[andreasnoack]

I think this one is ready.

[devmotion]

Did you see my suggestions in #314 (comment) and #314 (comment) and similar lines? It seems they were marked as resolved but maybe that happened automatically when you force pushed the branch.

[devmotion]

Looks good to me 👍

[andreasnoack]

I could have sworn that I'd already accepted all of them. Maybe it's a problem if you don't allow enough time between acceptance of different suggestions.