--- gsl-1.15-vanilla/randist/nbinomial.c	2010-12-26 10:57:08.000000000 -0700
+++ gsl-1.15/randist/nbinomial.c	2011-06-02 13:16:33.871584330 -0600
@@ -48,7 +48,7 @@ gsl_ran_negative_binomial_pdf (const uns
   double a = gsl_sf_lngamma (n) ;
   double b = gsl_sf_lngamma (k + 1.0) ;
 
-  P = exp(f-a-b) * pow (p, n) * pow (1 - p, (double)k);
-  
+  P = exp(f - a - b + n * log(p) + k * log(1 - p));
+
   return P;
 }