sf_gamma questions

[Jonathan Leto <jonathan@leto.net>]

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