[Issue 22502] New: Potential error where std.internal.math.gammafunction.betaIncompleteInv gives different results from Wolfram Alpha for parameters aa and bb being half integer and yy0 being 0.025.

[d-bugmail at puremagic.com]

https://issues.dlang.org/show_bug.cgi?id=22502

          Issue ID: 22502
           Summary: Potential error where
                    std.internal.math.gammafunction.betaIncompleteInv
                    gives different results from Wolfram Alpha for
                    parameters aa and bb being half integer and yy0  being
                    0.025.
           Product: D
           Version: D2
          Hardware: x86_64
                OS: Linux
            Status: NEW
          Severity: minor
          Priority: P1
         Component: phobos
          Assignee: nobody at puremagic.com
          Reporter: alex at dengenius.com

So I was trying to match a table in this article
(https://www.jstor.org/stable/2676784), in particular Table 5. which basically
corresponds to betaIncompleteInverse(x+1/2,N-x+1/2,0.05/2).  Using Wolfram
Alpha I can match the table, but I am getting numbers which seem surprisingly
off when using betaIncompleteInverse.

For example: 

InverseBetaRegularized[0.05/2, 7+1/2, 1+1/2] == 0.546280678967585... and
betaIncompleteInv(7+1/2,1+1/2,0.05/2) ~= 0.590384

So I did a few sanity checks using the test cases in the unittests in the file
(https://github.com/dlang/phobos/blob/v2.098.0/std/internal/math/gammafunction.d).
For example:

InverseBetaRegularized[0.0109343152340992, 8, 10] == 0.2 just like
betaIncompleteInv(8, 10, 0.010_934_315_234_099_2L) ~= 0.2


Based on that I think I am using the function correctly, it is still possible I
am doing something wrong though.

Using the program:
```
import std.mathspecial;
import std.stdio;

real jefferys(real x, real N, real alpha=0.05){
 return betaIncompleteInverse(x+1/2,N-x+1/2,alpha/2);
}


int main(){
    for(real N =7; N< 10; N++){
        for(real x=1; x < N;x++){
         writef("(%f,%f) %f\n", N,x,jefferys(x,N));
        }
    }
    return 0;
}
```
and Wolfram Alpha like
(https://www.wolframalpha.com/input/?i=InverseBetaRegularized%5B0.05%2F2%2C+1%2B1%2F2%2C+7-1%2B1%2F2%5D),
just varying the parameters and Table 5 from the article I linked I get:



 (N,x)             Wolfram               D             Table 5
 (7,1)             0.0158765...         0.004211..      0.016
 (7,2)             0.0647282...         0.043272..      0.065
 (7,3)             0.1388642...         0.118117..      0.139
 (7,4)             0.2345012...         0.222778..      0.234

And so on.

Sorry if this is just some kind of user error on my part.

--