gsl_cdf_binomial_P seems to choke on some, but not all, large values for k and n. gsl-2.6 built by evensen@shin-niryoko:~/src/gsl/gsl-2.6$ ./configure; make; make install evensen@shin-niryoko:~/src/gsl/gsl-2.6$ gcc --version gcc (Ubuntu 9.2.1-9ubuntu2) 9.2.1 20191008
# The snippet below exposes the bug: evensen@shin-niryoko:~/src/nasam$ cat binomcdf-error.c // Short program to expose NaN-bug in gsl_cdf_binomial_P: #include <math.h> #include <stdio.h> #include <gsl/gsl_cdf.h> #include <gsl/gsl_version.h> int main(int argc, char* argv[]) { printf("GSL version: %s\n", GSL_VERSION); unsigned int k = 4193750; unsigned int n = 8388600; double p = 0.5; int scale = 1; double bp; do { bp = gsl_cdf_binomial_P(k / scale, p, n / scale); printf("Pr(X <= %u) for %u trials with p = %f: %f\n", k / scale, n / scale, p, bp); scale *= 2; } while(isnan(bp)); return 0; } # Compiling the snippet: evensen@shin-niryoko:~/src/nasam$ gcc -Wall binomcdf-error.c -lgsl -lm -lblas evensen@shin-niryoko:~/src/nasam$ ./a.out GSL version: 2.6 Pr(X <= 8387760) for 16777216 trials with p = 0.500000: nan <=== NaN Pr(X <= 4193880) for 8388608 trials with p = 0.500000: nan <=== NaN Pr(X <= 2096940) for 4194304 trials with p = 0.500000: nan <=== NaN Pr(X <= 1048470) for 2097152 trials with p = 0.500000: 0.442078 Regards, Pelle Evensen