--- gamma.c.orig	Mon Dec  3 23:21:36 2001
+++ gamma.c	Mon Dec  3 23:53:27 2001
@@ -1247,7 +1247,14 @@
 int
 gsl_sf_gamma_e(const double x, gsl_sf_result * result)
 {
-  if(x < 0.5) {
+  int y;
+  y = abs((int)x - x );
+
+  /* if x is an integer as far as we can tell, return factorial(x-1) 
+     this gives us much better precision when x>19 or so */
+  if(y < GSL_DBL_EPSILON ){
+	gsl_sf_fact_e((int) x - 1,result);
+  } else if(x < 0.5) {
     int rint_x = (int)floor(x+0.5);
     double f_x = x - rint_x;
     double sgn = ( GSL_IS_EVEN(rint_x) ? 1.0 : -1.0 );