homepage

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.