sf_gamma questions

sf_gamma questions

[Jonathan Leto]

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

I was writing some test cases for Math::Gsl, when I saw some
failing, because gamma(19) was off by 1 and gamma(20) is off by
16. So I wrote a little test program to compare absolute difference
between gamma(x) and factorial(x-1) . The output was: 

-------------- Actual error vs. Theoretical Error
Difference (1):0 vs 2.22044604925031308e-16
Difference (2):0 vs 2.22044604925031308e-16
Difference (3):1.77635683940025046e-15 vs 1.4246133846675888e-14
Difference (4):7.99360577730112709e-15 vs 4.71293383088833937e-14
Difference (5):1.06581410364015028e-14 vs 6.92426984602974577e-14
Difference (6):2.84217094304040074e-14 vs 2.96806016881859386e-13
Difference (7):1.1368683772161603e-13 vs 1.45976589371855462e-12
Difference (8):9.09494701772928238e-13 vs 1.02398955184829092e-11
Difference (9):7.2759576141834259e-12 vs 1.01934226616975847e-10
Difference (10):5.82076609134674072e-11 vs 1.00719432793994222e-09
Difference (11):4.65661287307739258e-10 vs 1.08776987417513721e-08
Difference (12):0 vs 1.28517996245136601e-07
Difference (13):5.9604644775390625e-08 vs 1.64857567597209688e-06
Difference (14):0 vs 2.28141601610332145e-05
Difference (15):1.52587890625e-05 vs 0.000338755711482008391
Difference (16):0 vs 0.00537169771064327506
Difference (17):0 vs 0.090592955984902801
Difference (18):0.0625 vs 1.61905872619172486
Difference (19):1 vs 30.5646696115218255
Difference (20):16 vs 607.739360880259596
Difference (21):0 vs 12694.999982832087
Difference (22):0 vs 277939.467709238641
Difference (23):131072 vs 6364246.58713807818
Difference (24):1572864 vs 152117972.34747678
Difference (25):104857600 vs 3788598556.57866812

Why is the error so large for  gamma(x>18) ? Would it be possible
to just return factorial(x-1) if gamma is given an integer argument?

Test program (perl) is attached.

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



[-- Attachment #2: a --]
[-- Type: text/plain, Size: 318 bytes --]

#!/usr/bin/perl -w

use strict;
use Math::Gsl::Sf qw(gamma gamma_e);
my $r = new Math::Gsl::Sf::Result;


for (1 .. 25 ){
	gamma_e( $_ , $r );
	print "Difference ($_):" .  abs( $r->val - fact($_- 1) )  . " vs ";
	print $r->err . "\n";
}

sub fact {
	my $x = shift;
	return 1 if $x == 0;
	return $x * fact( $x - 1 );
}

Re: sf_gamma questions

[Dan, Ho-Jin]

I found that the results from slatec are slightly accurate compare to 
that of slatec.
The relative error of the order O(1.e-15) shows that both routines are 
very accurate.

[17!] error     slatec: 1.875000        [rel:0.000000]
                gsl   : -0.062500       [rel:-0.000000]
[18!] error     slatec: 6.000000        [rel:0.000000]
                gsl   : -1.000000       [rel:-0.000000]
[19!] error     slatec: -80.000000      [rel:-0.000000]
                gsl   : -16.000000      [rel:-0.000000]
[20!] error     slatec: 0.000000        [rel:0.000000]
                gsl   : 0.000000        [rel:0.000000]
[21!] error     slatec: -344064.000000  [rel:-0.000000]
                gsl   : 0.000000        [rel:0.000000]
[22!] error     slatec: 4194304.000000  [rel:0.000000]
                gsl   : -131072.000000  [rel:-0.000000]
[23!] error     slatec: -83886080.000000        [rel:-0.000000]
                gsl   : 0.000000        [rel:0.000000]
[24!] error     slatec: 4697620480.000000       [rel:0.000000]
                gsl   : -134217728.000000       [rel:-0.000000]
[25!] error     slatec: -79456894976.000000     [rel:-0.000000]
                gsl   : -2147483648.000000      [rel:-0.000000]


Jonathan Leto wrote:

>I was writing some test cases for Math::Gsl, when I saw some
>failing, because gamma(19) was off by 1 and gamma(20) is off by
>16. So I wrote a little test program to compare absolute difference
>between gamma(x) and factorial(x-1) . The output was: 
>

Re: sf_gamma questions

[Brian Gough]

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.

Re: sf_gamma questions

[Jonathan Leto]

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

Re: sf_gamma questions

[Brian Gough]

Jonathan Leto writes:
 > --- 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 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);

Thanks. I've added something along these lines but in the function
gamma_xgthalf, which is used in several places and so gives some other
opportunities to gain from a bit of extra accuracy.

Brian

Re: sf_gamma questions

[Brian Gough]

Jonathan Leto writes:
 > --- 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 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);

Thanks. I've added something along these lines but in the function
gamma_xgthalf, which is used in several places and so gives some other
opportunities to gain from a bit of extra accuracy.

Brian

sf_gamma questions

[Jonathan Leto]

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

I was writing some test cases for Math::Gsl, when I saw some
failing, because gamma(19) was off by 1 and gamma(20) is off by
16. So I wrote a little test program to compare absolute difference
between gamma(x) and factorial(x-1) . The output was: 

-------------- Actual error vs. Theoretical Error
Difference (1):0 vs 2.22044604925031308e-16
Difference (2):0 vs 2.22044604925031308e-16
Difference (3):1.77635683940025046e-15 vs 1.4246133846675888e-14
Difference (4):7.99360577730112709e-15 vs 4.71293383088833937e-14
Difference (5):1.06581410364015028e-14 vs 6.92426984602974577e-14
Difference (6):2.84217094304040074e-14 vs 2.96806016881859386e-13
Difference (7):1.1368683772161603e-13 vs 1.45976589371855462e-12
Difference (8):9.09494701772928238e-13 vs 1.02398955184829092e-11
Difference (9):7.2759576141834259e-12 vs 1.01934226616975847e-10
Difference (10):5.82076609134674072e-11 vs 1.00719432793994222e-09
Difference (11):4.65661287307739258e-10 vs 1.08776987417513721e-08
Difference (12):0 vs 1.28517996245136601e-07
Difference (13):5.9604644775390625e-08 vs 1.64857567597209688e-06
Difference (14):0 vs 2.28141601610332145e-05
Difference (15):1.52587890625e-05 vs 0.000338755711482008391
Difference (16):0 vs 0.00537169771064327506
Difference (17):0 vs 0.090592955984902801
Difference (18):0.0625 vs 1.61905872619172486
Difference (19):1 vs 30.5646696115218255
Difference (20):16 vs 607.739360880259596
Difference (21):0 vs 12694.999982832087
Difference (22):0 vs 277939.467709238641
Difference (23):131072 vs 6364246.58713807818
Difference (24):1572864 vs 152117972.34747678
Difference (25):104857600 vs 3788598556.57866812

Why is the error so large for  gamma(x>18) ? Would it be possible
to just return factorial(x-1) if gamma is given an integer argument?

Test program (perl) is attached.

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



[-- Attachment #2: a --]
[-- Type: text/plain, Size: 318 bytes --]

#!/usr/bin/perl -w

use strict;
use Math::Gsl::Sf qw(gamma gamma_e);
my $r = new Math::Gsl::Sf::Result;


for (1 .. 25 ){
	gamma_e( $_ , $r );
	print "Difference ($_):" .  abs( $r->val - fact($_- 1) )  . " vs ";
	print $r->err . "\n";
}

sub fact {
	my $x = shift;
	return 1 if $x == 0;
	return $x * fact( $x - 1 );
}

Re: sf_gamma questions

[Dan, Ho-Jin]

I found that the results from slatec are slightly accurate compare to 
that of slatec.
The relative error of the order O(1.e-15) shows that both routines are 
very accurate.

[17!] error     slatec: 1.875000        [rel:0.000000]
                gsl   : -0.062500       [rel:-0.000000]
[18!] error     slatec: 6.000000        [rel:0.000000]
                gsl   : -1.000000       [rel:-0.000000]
[19!] error     slatec: -80.000000      [rel:-0.000000]
                gsl   : -16.000000      [rel:-0.000000]
[20!] error     slatec: 0.000000        [rel:0.000000]
                gsl   : 0.000000        [rel:0.000000]
[21!] error     slatec: -344064.000000  [rel:-0.000000]
                gsl   : 0.000000        [rel:0.000000]
[22!] error     slatec: 4194304.000000  [rel:0.000000]
                gsl   : -131072.000000  [rel:-0.000000]
[23!] error     slatec: -83886080.000000        [rel:-0.000000]
                gsl   : 0.000000        [rel:0.000000]
[24!] error     slatec: 4697620480.000000       [rel:0.000000]
                gsl   : -134217728.000000       [rel:-0.000000]
[25!] error     slatec: -79456894976.000000     [rel:-0.000000]
                gsl   : -2147483648.000000      [rel:-0.000000]


Jonathan Leto wrote:

>I was writing some test cases for Math::Gsl, when I saw some
>failing, because gamma(19) was off by 1 and gamma(20) is off by
>16. So I wrote a little test program to compare absolute difference
>between gamma(x) and factorial(x-1) . The output was: 
>

Re: sf_gamma questions

[Brian Gough]

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.

Re: sf_gamma questions

[Jonathan Leto]

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

sf_gamma questions

[Jonathan Leto]

I was writing some test cases for Math::Gsl, when I saw some
failing, because gamma(19) was off by 1 and gamma(20) is off by
16. So I wrote a little test program to compare absolute difference
between gamma(x) and factorial(x-1) . The output was: 

-------------- Actual error vs. Theoretical Error
Difference (1):0 vs 2.22044604925031308e-16
Difference (2):0 vs 2.22044604925031308e-16
Difference (3):1.77635683940025046e-15 vs 1.4246133846675888e-14
Difference (4):7.99360577730112709e-15 vs 4.71293383088833937e-14
Difference (5):1.06581410364015028e-14 vs 6.92426984602974577e-14
Difference (6):2.84217094304040074e-14 vs 2.96806016881859386e-13
Difference (7):1.1368683772161603e-13 vs 1.45976589371855462e-12
Difference (8):9.09494701772928238e-13 vs 1.02398955184829092e-11
Difference (9):7.2759576141834259e-12 vs 1.01934226616975847e-10
Difference (10):5.82076609134674072e-11 vs 1.00719432793994222e-09
Difference (11):4.65661287307739258e-10 vs 1.08776987417513721e-08
Difference (12):0 vs 1.28517996245136601e-07
Difference (13):5.9604644775390625e-08 vs 1.64857567597209688e-06
Difference (14):0 vs 2.28141601610332145e-05
Difference (15):1.52587890625e-05 vs 0.000338755711482008391
Difference (16):0 vs 0.00537169771064327506
Difference (17):0 vs 0.090592955984902801
Difference (18):0.0625 vs 1.61905872619172486
Difference (19):1 vs 30.5646696115218255
Difference (20):16 vs 607.739360880259596
Difference (21):0 vs 12694.999982832087
Difference (22):0 vs 277939.467709238641
Difference (23):131072 vs 6364246.58713807818
Difference (24):1572864 vs 152117972.34747678
Difference (25):104857600 vs 3788598556.57866812

Why is the error so large for  gamma(x>18) ? Would it be possible
to just return factorial(x-1) if gamma is given an integer argument?

Test program (perl) is attached.

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

Re: sf_gamma questions

[Jonathan Leto]

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

Re: sf_gamma questions

[Brian Gough]

Dherminder Kainth writes:
 > I have been using the Levin U transform routines within a C++ program and 
 > have been getting memory leaks - for example the program crashes
 > when I attempt to deallocate the U transform workspace. A simpler
 > example program written in C seems to run fine. Is there any information
 > as to why this may be the case ? I can send both programs if somebody would 
 > like to have a look.
 > 
 > The coding is done on a RH 7.2 Linux box.

Hi,

I recommend recompiling your program and GSL with
checkerg++/checkergcc.  

That will trap any stack or malloc errors at the point where they
occur.

regards
Brian Gough

Re: sf_gamma questions

[Dherminder Kainth]

Happy new year !

I have been using the Levin U transform routines within a C++ program and 
have been getting memory leaks - for example the program crashes
when I attempt to deallocate the U transform workspace. A simpler
example program written in C seems to run fine. Is there any information
as to why this may be the case ? I can send both programs if somebody would 
like to have a look.

The coding is done on a RH 7.2 Linux box.

Thanks,

Dherminder Kainth

_________________________________________________________________
Send and receive Hotmail on your mobile device: http://mobile.msn.com

Re: sf_gamma questions

[Brian Gough]

Dherminder Kainth writes:
 > I have been using the Levin U transform routines within a C++ program and 
 > have been getting memory leaks - for example the program crashes
 > when I attempt to deallocate the U transform workspace. A simpler
 > example program written in C seems to run fine. Is there any information
 > as to why this may be the case ? I can send both programs if somebody would 
 > like to have a look.
 > 
 > The coding is done on a RH 7.2 Linux box.

Hi,

I recommend recompiling your program and GSL with
checkerg++/checkergcc.  

That will trap any stack or malloc errors at the point where they
occur.

regards
Brian Gough

Re: sf_gamma questions

[Dherminder Kainth]

Happy new year !

I have been using the Levin U transform routines within a C++ program and 
have been getting memory leaks - for example the program crashes
when I attempt to deallocate the U transform workspace. A simpler
example program written in C seems to run fine. Is there any information
as to why this may be the case ? I can send both programs if somebody would 
like to have a look.

The coding is done on a RH 7.2 Linux box.

Thanks,

Dherminder Kainth

_________________________________________________________________
Send and receive Hotmail on your mobile device: http://mobile.msn.com