Bug in GSL gsl_cdf_binomial_P

[Pelle Evensen]

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