Re: sf_gamma questions

[Jonathan Leto <jonathan@leto.net>]

[-- Attachment #1: Type: text/plain, Size: 1412 bytes --]

Patch is attached to implement my idea. I am not a C programmer
by trade, so let me know if I am tickling Nyarlathotep or
something. 

New output of my difference program:

-------------- Actual error vs. Theoretical Error
Difference (1):0 vs 0
Difference (2):0 vs 0
Difference (3):0 vs 0
Difference (4):0 vs 0
Difference (5):0 vs 0
Difference (6):0 vs 0
Difference (7):0 vs 0
Difference (8):0 vs 0
Difference (9):0 vs 0
Difference (10):0 vs 0
Difference (11):0 vs 0
Difference (12):0 vs 0
Difference (13):0 vs 0
Difference (14):0 vs 0
Difference (15):0 vs 0
Difference (16):0 vs 0
Difference (17):0 vs 0
Difference (18):0 vs 0
Difference (19):0 vs 2.84322508014156
Difference (20):0 vs 54.0212765226897
Difference (21):0 vs 1080.42553045379
Difference (22):0 vs 22688.9361395297
Difference (23):0 vs 499156.595069653
Difference (24):0 vs 11480601.686602
Difference (25):0 vs 275534440.478449

 Brian Gough (bjg@network-theory.co.uk) was saying:

> Jonathan Leto writes:
>  > Why is the error so large for gamma(x>18) ?
> 
> It's computed with an approximation so the relative error is the
> appropriate measure. Note that 1/Gamma(19) is the first value which is
> less than DBL_EPSILON.
> 
>  > Would it be possible to just return factorial(x-1) if gamma is
>  > given an integer argument?
> 
> Sounds reasonable to me.

-- 
jonathan@leto.net 
"Wir muessen wissen. Wir werden wissen."
	-- David Hilbert, 1931



[-- Attachment #2: gsl-1.0_gamma.patch --]
[-- Type: text/plain, Size: 565 bytes --]

--- 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 );