gsl_ran_negative_binomial_pdf gives -NaN if p == 1.0, despite p=1 being a valid negative binomial distribution. I believe this is because of the call to log1p(-p) which evaluates to log(0) and gives NaN. There appear to be guards against this in the positive version (gsl_ran_binomial_pdf) and it works as expected, however these edge cases aren't covered in the negative version. I think that the negative version should have similar guards to return appropriate values for edge cases like p==1.
Thanks for your consideration, Eric Redekopp University of Saskatchewan
