New patch for gamma , with some tweaks:

 - return exact values for integral arguments: gamma(1) through gamma(23)
 - apply a cheap correction to improve accuracy of exp and pow
   computations
 - use a different form of the reflection formula:

     gamma(x) = -pi/sinpi(x)/x/gamma(x)

   (the usual reflection formula has accuracy problems for
    x close to a power of 2;  e.g., x in (-64,-63) or x
    in (-128, -127))

 - avoid duplication formula for large negative arguments
 - add a few extra tests

On my machine, testing with approx. 10**7 random samples, this version 
achieves an accuracy of <= 10 ulps across the domain (comparing with 
correctly-rounded results generated by MPFR).  Limiting the test to 
arguments in the range (-256.0, 1/256.0] + [1/256.0, 256.0) (with each 
float in that range equally likely), the error in ulps from 10**6 samples 
has mean -0.104 and standard deviation 1.230.

I plan to check this in in a week or two.  Feedback welcome!  It would be 
especially useful to know whether this patch compiles correctly on 
Windows.