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