This is the mail archive of the libc-alpha@sourceware.org mailing list for the glibc project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

Re: Improve tgamma accuracy (bugs 2546, 2560, 5159, 15426)


On 05/02/2013 06:35 PM, Joseph S. Myers wrote:
> This patch fixes various bugs relating to inaccuracy of results from
> tgamma.  Bugs 2546, 2560 and 5159 are different cases of inaccurate
> finite results, arising from the use of exp (lgamma) which inevitably
> results in large errors for results with large (positive or negative)
> exponent.  Bug 15426, found by the new testcases prepared for this
> patch, is wrong sign of infinity on overflow for negative arguments
> close to 0.  (Bugs 15425 and 15427 were also found by testing unfixed
> glibc with the testcases included in this patch; this patch does
> nothing about those bugs in lgammaf and exp, although they no longer
> affect the results of tgamma after the tgamma implementation changes
> in this patch.)
> 
> Because of the inaccuracy resulting from exp (lgamma) when the result
> of lgamma is large, after this patch lgamma is only called by tgamma
> on arguments in the range 0.5 to 1.5.  For larger positive arguments,
> the recurrence gamma(z+1) = z * gamma(z) is used to reduce either to a
> value in that range where exp (lgamma) does not introduce large
> errors, or to a value large enough to apply Stirling's approximation.
> When using Stirling's approximation, pow is then used to avoid large
> errors from the use of exp on large values, and care is taken that
> errors in arguments of exp and exp2 are small (measured in ulps of 1,
> not ulps of the argument of exp, since ulps of 1 are what's relevant
> to avoid large errors in the final result); a power of 2 is split out
> in the calculation to avoid internal overflows when the final result
> is finite.  Negative arguments (other than those very close to 0) are
> handled through gamma(z) = pi/(-z)sin(pi*z)gamma(-z).
> 
> To reduce accumulation of errors in the reduction, __gamma_product
> keeps an estimate of how much the exact product differs from the
> computed value.  (The estimation in the ldbl-128ibm version is less
> thorough; it deals with errors arising from an adjusted version of the
> argument being used in the product because the result of adding a
> small integer to the unadjusted argument may not be exactly
> representable, but not with errors from each successive multiplication
> step.  The former kind of errors would all be in the same direction,
> whereas the latter would generally be in random directions, but this
> is still expected to make the ldbl-128ibm version slightly less
> accurate - although such reduced accuracy doesn't show in the
> testsuite.)
> 
> Before this patch, on x86_64, the testcases added by this patch showed
> errors of up to 203 ulp for float, 1713 ulp for double and 12714 ulp
> for long double; the patch reduces the maximum errors for its
> testcases to 3 ulp, 4 ulp and 4 ulp respectively.
> 
> Tested x86_64 and x86 and ulps updated accordingly.  Also did spot
> tests on mips64 to sanity test the ldbl-128 version, and powerpc32 to
> sanity test the ldbl-128ibm version.  The former showed errors of at
> most 4 ulp and the latter errors of at most 3 ulp.  There are spurious
> underflows for ldbl-128ibm testing inputs 0x1p-1023, 0x1p-1024,
> 0x1p-1074, which I think (like some other underflows seen on tests of
> other functions) just reflect problems with the underlying long double
> arithmetic giving spurious underflow exceptions when dealing with
> norma long doubles whose lower part is a subnormal double (because of
> the discontiguous mantissa bits), not actual problems with libm itself
> (such issues would need fixing in libgcc's multiplication / division
> code).
> 
> (There are other open bugs about inaccuracy of lgamma.  Those have an
> independent cause - large cancellation when the gamma function has
> value close to +/- 1 - and will need their own fix.  The IA64 version
> of lgammal addresses those by using 35 different degree-19 polynomials
> for values near roots of lgamma (less needed for float and double
> versions, of course); I haven't looked in detail at how feasible it is
> to be accurate near roots of lgamma without so much precomputed data.)
> 
> 2013-05-02  Joseph Myers  <joseph@codesourcery.com>
> 
> 	[BZ #2546]
> 	[BZ #2560]
> 	[BZ #5159]
> 	[BZ #15426]
> 	* sysdeps/ieee754/k_standard.c (__kernel_standard): Copy sign of
> 	input to result for tgamma overflow.
> 	* sysdeps/ieee754/dbl-64/e_gamma_r.c: Include <float.h>.
> 	(gamma_coeff): New variable.
> 	(NCOEFF): New macro.
> 	(gamma_positive): New function.
> 	(__ieee754_gamma_r): Handle positive infinity, NaN, overflow and
> 	underflow here.  Use gamma_positive instead of exp (lgamma) for
> 	other arguments.
> 	* sysdeps/ieee754/flt-32/e_gammaf_r.c: Include <float.h>.
> 	(gamma_coeff): New variable.
> 	(NCOEFF): New macro.
> 	(gammaf_positive): New function.
> 	(__ieee754_gammaf_r): Handle positive infinity, NaN, overflow and
> 	underflow here.  Use gamma_positive instead of exp (lgamma) for
> 	other arguments.
> 	* sysdeps/ieee754/ldbl-128/e_gammal_r.c: Include <float.h>.
> 	(gamma_coeff): New variable.
> 	(NCOEFF): New macro.
> 	(gammal_positive): New function.
> 	(__ieee754_gammal_r): Handle positive infinity, NaN, overflow and
> 	underflow here.  Use gamma_positive instead of exp (lgamma) for
> 	other arguments.
> 	* sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c: Include <float.h>.
> 	(gamma_coeff): New variable.
> 	(NCOEFF): New macro.
> 	(gammal_positive): New function.
> 	(__ieee754_gammal_r): Handle positive infinity, overflow and
> 	underflow here.  Handle NaN the same as positive infinity.  Remove
> 	check x < 0xffffffff for negative integers.  Use gamma_positive
> 	instead of exp (lgamma) for other arguments.
> 	* sysdeps/ieee754/ldbl-96/e_gammal_r.c: Include <float.h>.
> 	(gamma_coeff): New variable.
> 	(NCOEFF): New macro.
> 	(gammal_positive): New function.
> 	* sysdeps/ieee754/dbl-64/gamma_product.c: New file.
> 	* sysdeps/ieee754/dbl-64/gamma_productf.c: Likewise.
> 	* sysdeps/ieee754/ldbl-128/gamma_productl.c: Likewise.
> 	* sysdeps/ieee754/ldbl-128ibm/gamma_productl.c: Likewise.
> 	* sysdeps/ieee754/ldbl-96/gamma_product.c: Likewise.
> 	* sysdeps/ieee754/ldbl-96/gamma_productl.c: Likewise.
> 	* sysdeps/generic/math_private.h (__gamma_productf): New
> 	prototype.
> 	(__gamma_product): Likewise.
> 	(__gamma_productl): Likewise.
> 	* math/Makefile (libm-calls): Add gamma_product.
> 	* math/libm-test.inc (tgamma_test): Add more tests.
> 	* sysdeps/i386/fpu/libm-test-ulps: Update.
> 	* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.

This is a very very large patch, which I have not reviewed completely,
but will attempt to provide some review and comments for.


> diff --git a/math/Makefile b/math/Makefile
> index 9f0bf72..5bbf9d3 100644
> --- a/math/Makefile
> +++ b/math/Makefile
> @@ -58,7 +58,8 @@ libm-calls = e_acos e_acosh e_asin e_atan2 e_atanh e_cosh e_exp e_fmod	\
>  	     s_catan s_casin s_ccos s_csin s_ctan s_ctanh s_cacos	\
>  	     s_casinh s_cacosh s_catanh s_csqrt s_cpow s_cproj s_clog10 \
>  	     s_fma s_lrint s_llrint s_lround s_llround e_exp10 w_log2	\
> -	     s_isinf_ns s_issignaling $(calls:s_%=m_%) x2y2m1 k_casinh
> +	     s_isinf_ns s_issignaling $(calls:s_%=m_%) x2y2m1 k_casinh	\
> +	     gamma_product
>  
>  include ../Makeconfig
>  
> diff --git a/math/libm-test.inc b/math/libm-test.inc
> index 8feb0a6..46f42b6 100644
> --- a/math/libm-test.inc
> +++ b/math/libm-test.inc
> @@ -13482,11 +13482,511 @@ tgamma_test (void)
>    TEST_f_f (tgamma, -0.5, -M_2_SQRT_PIl);
>  
>    TEST_f_f (tgamma, 1, 1);
> +  TEST_f_f (tgamma, 2, 1);
> +  TEST_f_f (tgamma, 3, 2);
>    TEST_f_f (tgamma, 4, 6);
> +  TEST_f_f (tgamma, 5, 24);
> +  TEST_f_f (tgamma, 6, 120);
> +  TEST_f_f (tgamma, 7, 720);
> +  TEST_f_f (tgamma, 8, 5040);
> +  TEST_f_f (tgamma, 9,  40320);
> +  TEST_f_f (tgamma, 10, 362880);
>  
>    TEST_f_f (tgamma, 0.7L, 1.29805533264755778568117117915281162L);
>    TEST_f_f (tgamma, 1.2L, 0.918168742399760610640951655185830401L);
>  
> +  TEST_f_f (tgamma, 1.5L, 8.8622692545275801364908374167057259139877e-01L);
> +  TEST_f_f (tgamma, 2.5L, 1.3293403881791370204736256125058588870982e+00L);
> +  TEST_f_f (tgamma, 3.5L, 3.3233509704478425511840640312646472177454e+00L);
> +  TEST_f_f (tgamma, 4.5L, 1.1631728396567448929144224109426265262109e+01L);
> +  TEST_f_f (tgamma, 5.5L, 5.2342777784553520181149008492418193679490e+01L);
> +  TEST_f_f (tgamma, 6.5L, 2.8788527781504436099631954670830006523720e+02L);
> +  TEST_f_f (tgamma, 7.5L, 1.8712543057977883464760770536039504240418e+03L);
> +  TEST_f_f (tgamma, 8.5L, 1.4034407293483412598570577902029628180313e+04L);
> +  TEST_f_f (tgamma, 9.5L, 1.1929246199460900708784991216725183953266e+05L);
> +  TEST_f_f (tgamma, -1.5L, 2.3632718012073547030642233111215269103967e+00L);
> +  TEST_f_f (tgamma, -2.5L, -9.4530872048294188122568932444861076415869e-01L);
> +  TEST_f_f (tgamma, -3.5L, 2.7008820585226910892162552127103164690248e-01L);
> +  TEST_f_f (tgamma, -4.5L, -6.0019601300504246427027893615784810422774e-02L);
> +  TEST_f_f (tgamma, -5.5L, 1.0912654781909862986732344293779056440504e-02L);
> +  TEST_f_f (tgamma, -6.5L, -1.6788699664476712287280529682737009908468e-03L);
> +  TEST_f_f (tgamma, -7.5L, 2.2384932885968949716374039576982679877958e-04L);
> +  TEST_f_f (tgamma, -8.5L, -2.6335215159963470254557693619979623385833e-05L);
> +  TEST_f_f (tgamma, -9.5L, 2.7721279115751021320587045915768024616666e-06L);

How did you choose the number of digits in the answer?

How did you compute the answer?

A quick check for long double on x86-64 shows that a correct answer rounded
to long double would required this many digits (more than you have):

2.772127911575102132076012628788120976952580321039931732229888439178466796875e-06

or

0xb.a08cb7d65f3679ep-22

Which is 35 more digits than your answer.

It is also different than your answer by ~2e-26.

Computed using Wolfram's Alpha.

http://www.wolframalpha.com/input/?i=gamma%28-9.5%29

What am I misunderstanding about the way libm-test.inc works?

> +  TEST_f_f (tgamma, 0x1p-24L, 1.6777215422784394050795179874582764575261e+07L);
> +  TEST_f_f (tgamma, -0x1p-24L, -1.6777216577215723853867349114260580375249e+07L);
> +  TEST_f_f (tgamma, 0x1p-53L, 9.0071992547409914227843350984672492007618e+15L);
> +  TEST_f_f (tgamma, -0x1p-53L, -9.0071992547409925772156649015329704137860e+15L);
> +  TEST_f_f (tgamma, 0x1p-64L, 1.8446744073709551615422784335098467139447e+19L);
> +  TEST_f_f (tgamma, -0x1p-64L, -1.8446744073709551616577215664901532860660e+19L);
> +  TEST_f_f (tgamma, 0x1p-106L, 8.1129638414606681695789005144063422784335e+31L);
> +  TEST_f_f (tgamma, -0x1p-106L, -8.1129638414606681695789005144064577215665e+31L);
> +  TEST_f_f (tgamma, 0x1p-113L, 1.0384593717069655257060992658440191422784e+34L);
> +  TEST_f_f (tgamma, -0x1p-113L, -1.0384593717069655257060992658440192577216e+34L);
> +  TEST_f_f (tgamma, 0x1p-127L, 1.7014118346046923173168730371588410572742e+38L);
> +  TEST_f_f (tgamma, -0x1p-127L, -1.7014118346046923173168730371588410572858e+38L);
> +#ifdef TEST_FLOAT
> +  TEST_f_f (tgamma, 0x1p-128L, plus_infty, OVERFLOW_EXCEPTION);
> +  TEST_f_f (tgamma, -0x1p-128L, minus_infty, OVERFLOW_EXCEPTION);
> +  TEST_f_f (tgamma, 0x1p-149L, plus_infty, OVERFLOW_EXCEPTION);
> +  TEST_f_f (tgamma, -0x1p-149L, minus_infty, OVERFLOW_EXCEPTION);
> +#else
> +  TEST_f_f (tgamma, 0x1p-128L, 3.4028236692093846346337460743176821145542e+38L);
> +  TEST_f_f (tgamma, -0x1p-128L, -3.4028236692093846346337460743176821145658e+38L);
> +  TEST_f_f (tgamma, 0x1p-149L, 7.1362384635297994052914298472474756819137e+44L);
> +  TEST_f_f (tgamma, -0x1p-149L, -7.1362384635297994052914298472474756819137e+44L);
> +#endif
> +#ifndef TEST_FLOAT
> +  TEST_f_f (tgamma, 0x1p-1023L, 8.9884656743115795386465259539451236680899e+307L);
> +  TEST_f_f (tgamma, -0x1p-1023L, -8.9884656743115795386465259539451236680899e+307L);
> +# if !defined TEST_LDOUBLE || LDBL_MAX_EXP <= 1024
> +  TEST_f_f (tgamma, 0x1p-1024L, plus_infty, OVERFLOW_EXCEPTION);
> +  TEST_f_f (tgamma, -0x1p-1024L, minus_infty, OVERFLOW_EXCEPTION);
> +  TEST_f_f (tgamma, 0x1p-1074L, plus_infty, OVERFLOW_EXCEPTION);
> +  TEST_f_f (tgamma, -0x1p-1074L, minus_infty, OVERFLOW_EXCEPTION);
> +# else
> +  TEST_f_f (tgamma, 0x1p-1024L, 1.7976931348623159077293051907890247336180e+308L);
> +  TEST_f_f (tgamma, -0x1p-1024L, -1.7976931348623159077293051907890247336180e+308L);
> +  TEST_f_f (tgamma, 0x1p-1074L, 2.0240225330731061835249534671891730704956e+323L);
> +  TEST_f_f (tgamma, -0x1p-1074L, -2.0240225330731061835249534671891730704956e+323L);
> +# endif
> +#endif
> +#if defined TEST_LDOUBLE && LDBL_MIN_EXP <= -16381
> +  TEST_f_f (tgamma, 0x1p-16383L, 5.9486574767861588254287966331400356538172e+4931L);
> +  TEST_f_f (tgamma, -0x1p-16383L, -5.9486574767861588254287966331400356538172e+4931L);
> +  TEST_f_f (tgamma, 0x1p-16384L, plus_infty, OVERFLOW_EXCEPTION);
> +  TEST_f_f (tgamma, -0x1p-16384L, minus_infty, OVERFLOW_EXCEPTION);
> +  TEST_f_f (tgamma, 0x1p-16445L, plus_infty, OVERFLOW_EXCEPTION);
> +  TEST_f_f (tgamma, -0x1p-16445L, minus_infty, OVERFLOW_EXCEPTION);
> +# if LDBL_MANT_DIG >= 113
> +  TEST_f_f (tgamma, 0x1p-16494L, plus_infty, OVERFLOW_EXCEPTION);
> +  TEST_f_f (tgamma, -0x1p-16494L, minus_infty, OVERFLOW_EXCEPTION);
> +# endif
> +#endif
> +  TEST_f_f (tgamma, 0x8.00001p0L, 5.0400096882277802019946778420223050233915e+03L);
> +  TEST_f_f (tgamma, 0x7.fffff8p0L, 5.0399951558933225045148935487583089307135e+03L);
> +  TEST_f_f (tgamma, 0x7.000008p0L, 7.2000064296977505705636258629805621178587e+02L);
> +  TEST_f_f (tgamma, 0x6.fffff8p0L, 7.1999935703082425988147448928288557689866e+02L);
> +  TEST_f_f (tgamma, 0x6.000008p0L, 1.2000009762487825358530770343720418162783e+02L);
> +  TEST_f_f (tgamma, 0x5.fffff8p0L, 1.1999990237520611552119807476573441975106e+02L);
> +  TEST_f_f (tgamma, 0x5.000008p0L, 2.4000017236155647574166073485628713443799e+01L);
> +  TEST_f_f (tgamma, 0x4.fffff8p0L, 2.3999982763857938712639837029547357501709e+01L);
> +  TEST_f_f (tgamma, 0x4.000008p0L, 6.0000035937827461765660468073471093546129e+00L);
> +  TEST_f_f (tgamma, 0x3.fffffcp0L, 5.9999982031095793171233994481968816873643e+00L);
> +  TEST_f_f (tgamma, 0x3.000004p0L, 2.0000004400179308360529417942462250547999e+00L);
> +  TEST_f_f (tgamma, 0x2.fffffcp0L, 1.9999995599822108706107786027549565954046e+00L);
> +  TEST_f_f (tgamma, 0x2.000004p0L, 1.0000001007996638509889062631687945799175e+00L);
> +  TEST_f_f (tgamma, 0x1.fffffep0L, 9.9999994960018563231526611134590489120697e-01L);
> +  TEST_f_f (tgamma, 0x1.000002p0L, 9.9999993119054472483596471908942669644327e-01L);
> +  TEST_f_f (tgamma, 0x0.ffffffp0L, 1.0000000344047381790797322460568297132998e+00L);
> +  TEST_f_f (tgamma, -0x0.ffffffp0L, -1.6777216422784419250710305882992376932423e+07L);
> +  TEST_f_f (tgamma, -0x1.000002p0L, 8.3886075772158332060084424806449513922858e+06L);
> +  TEST_f_f (tgamma, -0x1.fffffep0L, 4.1943044613922792026014320172298377770578e+06L);
> +  TEST_f_f (tgamma, -0x2.000004p0L, -2.0971515386080557574407223895988378776747e+06L);
> +  TEST_f_f (tgamma, -0x2.fffffcp0L, -6.9905087601970247876992248591045142913324e+05L);
> +  TEST_f_f (tgamma, -0x3.000004p0L, 6.9905045731381300146131914617735687322025e+05L);
> +  TEST_f_f (tgamma, -0x3.fffffcp0L, 1.7476272942159602684441970627092458855771e+05L);
> +  TEST_f_f (tgamma, -0x4.000008p0L, -8.7381270578483499672965708923121931082305e+04L);
> +  TEST_f_f (tgamma, -0x4.fffff8p0L, -1.7476280884325863043793087474680780379554e+04L);
> +  TEST_f_f (tgamma, -0x5.000008p0L, 1.7476252449031389167286893378510439443844e+04L);
> +  TEST_f_f (tgamma, -0x5.fffff8p0L, 2.9127137122026653716311560165769071985443e+03L);
> +  TEST_f_f (tgamma, -0x6.000008p0L, -2.9127085100239567622341538102130981196910e+03L);
> +  TEST_f_f (tgamma, -0x6.fffff8p0L, -4.1610198723079349791939054365613377035519e+02L);
> +  TEST_f_f (tgamma, -0x7.000008p0L, 4.1610118737306415004517215226199741948733e+02L);
> +  TEST_f_f (tgamma, -0x7.fffff8p0L, 5.2012751504050764429534086402871289946986e+01L);
> +  TEST_f_f (tgamma, -0x8.00001p0L, -2.6006296115134418896533598545925084576702e+01L);
> +  TEST_f_f (tgamma, -0x9.fffffp0L, 2.8896008370721717567612135720915723136310e-01L);
> +  TEST_f_f (tgamma, -0xa.00001p0L, -2.8895878754728051776830454190076999107021e-01L);
> +  TEST_f_f (tgamma, -0x13.ffffep0L, 2.1550026214525536756224040483579183652119e-13L);
> +  TEST_f_f (tgamma, -0x14.00002p0L, -2.1549777908265594916405421768142757507179e-13L);
> +  TEST_f_f (tgamma, -0x1d.ffffep0L, 1.9765721589464867957912772592816027583176e-27L);
> +  TEST_f_f (tgamma, -0x1e.00002p0L, -1.9765463890341964384070157599286498212650e-27L);
> +  TEST_f_f (tgamma, -0x27.ffffcp0L, 3.2129279441390812141195076945616975790225e-43L, UNDERFLOW_EXCEPTION_FLOAT);
> +  TEST_f_f (tgamma, -0x28.00004p0L, -3.2128372159115252365699015758097981155793e-43L, UNDERFLOW_EXCEPTION_FLOAT);
> +  TEST_f_f (tgamma, -0x28.ffffcp0L, -7.8364103489619817539676737414096652170685e-45L, UNDERFLOW_EXCEPTION_FLOAT);
> +  TEST_f_f (tgamma, -0x29.00004p0L, 7.8361876024016854597745353972619195760515e-45L, UNDERFLOW_EXCEPTION_FLOAT);
> +#ifdef TEST_FLOAT
> +  TEST_f_f (tgamma, -0x29.ffffcp0L, plus_zero, UNDERFLOW_EXCEPTION_FLOAT);
> +  TEST_f_f (tgamma, -0x2a.00004p0L, minus_zero, UNDERFLOW_EXCEPTION_FLOAT);
> +#else
> +  TEST_f_f (tgamma, -0x29.ffffcp0L, 1.8658121573125798145204120066590953505132e-46L);
> +  TEST_f_f (tgamma, -0x2a.00004p0L, -1.8657587834931410688246126853566488626385e-46L);
> +#endif
> +#ifndef TEST_FLOAT
> +  TEST_f_f (tgamma, 0x8.0000000000008p0L, 5.0400000000000180457125667322294144477136e+03L);
> +  TEST_f_f (tgamma, 0x7.ffffffffffffcp0L, 5.0399999999999909771437166339103165198442e+03L);
> +  TEST_f_f (tgamma, 0x7.0000000000004p0L, 7.2000000000000119762397445457359071259652e+02L);
> +  TEST_f_f (tgamma, 0x6.ffffffffffffcp0L, 7.1999999999999880237602554542848858572672e+02L);
> +  TEST_f_f (tgamma, 0x6.0000000000004p0L, 1.2000000000000018184042734842640022086408e+02L);
> +  TEST_f_f (tgamma, 0x5.ffffffffffffcp0L, 1.1999999999999981815957265157389249327533e+02L);
> +  TEST_f_f (tgamma, 0x5.0000000000004p0L, 2.4000000000000032104829055124673225982803e+01L);
> +  TEST_f_f (tgamma, 0x4.ffffffffffffcp0L, 2.3999999999999967895170944875373910918544e+01L);
> +  TEST_f_f (tgamma, 0x4.0000000000004p0L, 6.0000000000000066939396342309789716341613e+00L);
> +  TEST_f_f (tgamma, 0x3.ffffffffffffep0L, 5.9999999999999966530301828845138185025345e+00L);
> +  TEST_f_f (tgamma, 0x3.0000000000002p0L, 2.0000000000000008195971324717875960213536e+00L);
> +  TEST_f_f (tgamma, 0x2.ffffffffffffep0L, 1.9999999999999991804028675282128956223990e+00L);
> +  TEST_f_f (tgamma, 0x2.0000000000002p0L, 1.0000000000000001877539613108624482361963e+00L);
> +  TEST_f_f (tgamma, 0x1.fffffffffffffp0L, 9.9999999999999990612301934456883679778984e-01L);
> +  TEST_f_f (tgamma, 0x1.0000000000001p0L, 9.9999999999999987183237573039992418700531e-01L);
> +  TEST_f_f (tgamma, 0x0.fffffffffffff8p0L, 1.0000000000000000640838121348000744796665e+00L);
> +  TEST_f_f (tgamma, -0x0.fffffffffffff8p0L, -9.0071992547409924227843350984672961392521e+15L);
> +  TEST_f_f (tgamma, -0x1.0000000000001p0L, 4.5035996273704955772156649015331740980405e+15L);
> +  TEST_f_f (tgamma, -0x1.fffffffffffffp0L, 2.2517998136852484613921675492337776673289e+15L);
> +  TEST_f_f (tgamma, -0x2.0000000000002p0L, -1.1258999068426235386078324507668462444260e+15L);
> +  TEST_f_f (tgamma, -0x2.ffffffffffffep0L, -3.7529996894754154268627807196691509198813e+14L);
> +  TEST_f_f (tgamma, -0x3.0000000000002p0L, 3.7529996894754112398038859470009084971438e+14L);
> +  TEST_f_f (tgamma, -0x3.ffffffffffffep0L, 9.3824992236885396088236184658402406857503e+13L);
> +  TEST_f_f (tgamma, -0x4.0000000000004p0L, -4.6912496118442603911763815341745722862351e+13L);
> +  TEST_f_f (tgamma, -0x4.ffffffffffffcp0L, -9.3824992236885475509805702650262155809819e+12L);
> +  TEST_f_f (tgamma, -0x5.0000000000004p0L, 9.3824992236885191156860964016850034672946e+12L);
> +  TEST_f_f (tgamma, -0x5.ffffffffffffcp0L, 1.5637498706147581566449098589862357835505e+12L);
> +  TEST_f_f (tgamma, -0x6.0000000000004p0L, -1.5637498706147529544662012521330708016396e+12L);
> +  TEST_f_f (tgamma, -0x6.ffffffffffffcp0L, -2.2339283865925119357965832452642909859289e+11L);
> +  TEST_f_f (tgamma, -0x7.0000000000004p0L, 2.2339283865925039372192897706214475877342e+11L);
> +  TEST_f_f (tgamma, -0x7.ffffffffffffcp0L, 2.7924104832406402297655703264222230055898e+10L);
> +  TEST_f_f (tgamma, -0x8.0000000000008p0L, -1.3962052416203121511868106259843527348026e+10L);
> +  TEST_f_f (tgamma, -0x9.ffffffffffff8p0L, 1.5513391573559147700413058496716749249803e+08L);
> +  TEST_f_f (tgamma, -0xa.0000000000008p0L, -1.5513391573559018084419393002828541166901e+08L);
> +  TEST_f_f (tgamma, -0x13.ffffffffffffp0L, 1.1569515572952029402736625857313236848570e-04L);
> +  TEST_f_f (tgamma, -0x14.000000000001p0L, -1.1569515572951781096476686854873801225397e-04L);
> +  TEST_f_f (tgamma, -0x1d.ffffffffffffp0L, 1.0611571800204311628217068863959963842891e-18L);
> +  TEST_f_f (tgamma, -0x1e.000000000001p0L, -1.0611571800204053929094168642022073530425e-18L);
> +  TEST_f_f (tgamma, -0x27.fffffffffffep0L, 1.7249032006742266376460389310340465554361e-34L);
> +  TEST_f_f (tgamma, -0x28.000000000002p0L, -1.7249032006741359094184881234822934593822e-34L);
> +  TEST_f_f (tgamma, -0x28.fffffffffffep0L, -4.2070809772542120404320040128839297118648e-36L);
> +  TEST_f_f (tgamma, -0x29.000000000002p0L, 4.2070809772539892938717205103652583609422e-36L);
> +  TEST_f_f (tgamma, -0x29.fffffffffffep0L, 1.0016859469652887505173040814397197718981e-37L);
> +  TEST_f_f (tgamma, -0x2a.000000000002p0L, -1.0016859469652353766978684241048308120274e-37L);
> +  TEST_f_f (tgamma, -0x31.fffffffffffep0L, 4.6273774273632946947805289899230181990085e-51L);
> +  TEST_f_f (tgamma, -0x32.000000000002p0L, -4.6273774273630367887073532197576655720178e-51L);
> +  TEST_f_f (tgamma, -0x63.fffffffffffcp0L, 7.5400833348840965463348754984345825364294e-145L);
> +  TEST_f_f (tgamma, -0x64.000000000004p0L, -7.5400833348831085791638490135462230991587e-145L);
> +  TEST_f_f (tgamma, -0x95.fffffffffff8p0L, 6.1582369322723207086020016423767264008839e-250L);
> +  TEST_f_f (tgamma, -0x96.000000000008p0L, -6.1582369322705655439003240743176243138734e-250L);
> +  TEST_f_f (tgamma, -0xb4.fffffffffff8p0L, -9.6760879059917574597728750098636253931457e-319L, UNDERFLOW_EXCEPTION_DOUBLE);
> +  TEST_f_f (tgamma, -0xb5.000000000008p0L, 9.6760879059888966544677044221698800670218e-319L, UNDERFLOW_EXCEPTION_DOUBLE);
> +  TEST_f_f (tgamma, -0xb5.fffffffffff8p0L, 5.3165318164789884455066481673086605454904e-321L, UNDERFLOW_EXCEPTION_DOUBLE);
> +  TEST_f_f (tgamma, -0xb6.000000000008p0L, -5.3165318164774149139661976747137185876909e-321L, UNDERFLOW_EXCEPTION_DOUBLE);
> +  TEST_f_f (tgamma, -0xb6.fffffffffff8p0L, -2.9052086428846935908287469917922960610289e-323L, UNDERFLOW_EXCEPTION_DOUBLE);
> +  TEST_f_f (tgamma, -0xb7.000000000008p0L, 2.9052086428838328351126988236541632950925e-323L, UNDERFLOW_EXCEPTION_DOUBLE);
> +# if !defined TEST_LDOUBLE || LDBL_MAX_EXP <= 1024
> +  TEST_f_f (tgamma, -0xb7.fffffffffff8p0L, plus_zero, UNDERFLOW_EXCEPTION_DOUBLE);
> +  TEST_f_f (tgamma, -0xb8.000000000008p0L, minus_zero, UNDERFLOW_EXCEPTION_DOUBLE);
> +# else
> +  TEST_f_f (tgamma, -0xb7.fffffffffff8p0L, 1.5789177406982032823826953250736039527543e-325L);
> +  TEST_f_f (tgamma, -0xb8.000000000008p0L, -1.5789177406977349925854817486109369828857e-325L);
> +# endif
> +#endif
> +#if defined TEST_LDOUBLE && LDBL_MANT_DIG >= 106
> +  TEST_f_f (tgamma, 0x8.00000000000000000000000004p0L, 5.0400000000000000000000000000020034765587e+03L);
> +  TEST_f_f (tgamma, 0x7.fffffffffffffffffffffffffep0L, 5.0399999999999999999999999999989982617206e+03L);
> +  TEST_f_f (tgamma, 0x7.00000000000000000000000002p0L, 7.2000000000000000000000000000013296297113e+02L);
> +  TEST_f_f (tgamma, 0x6.fffffffffffffffffffffffffep0L, 7.1999999999999999999999999999986703702887e+02L);
> +  TEST_f_f (tgamma, 0x6.00000000000000000000000002p0L, 1.2000000000000000000000000000002018834292e+02L);
> +  TEST_f_f (tgamma, 0x5.fffffffffffffffffffffffffep0L, 1.1999999999999999999999999999997981165708e+02L);
> +  TEST_f_f (tgamma, 0x5.00000000000000000000000002p0L, 2.4000000000000000000000000000003564352042e+01L);
> +  TEST_f_f (tgamma, 0x4.fffffffffffffffffffffffffep0L, 2.3999999999999999999999999999996435647958e+01L);
> +  TEST_f_f (tgamma, 0x4.00000000000000000000000002p0L, 6.0000000000000000000000000000007431765907e+00L);
> +  TEST_f_f (tgamma, 0x3.ffffffffffffffffffffffffffp0L, 5.9999999999999999999999999999996284117046e+00L);
> +  TEST_f_f (tgamma, 0x3.00000000000000000000000001p0L, 2.0000000000000000000000000000000909935607e+00L);
> +  TEST_f_f (tgamma, 0x2.ffffffffffffffffffffffffffp0L, 1.9999999999999999999999999999999090064393e+00L);
> +  TEST_f_f (tgamma, 0x2.00000000000000000000000001p0L, 1.0000000000000000000000000000000208448771e+00L);
> +  TEST_f_f (tgamma, 0x1.ffffffffffffffffffffffffff8p0L, 9.9999999999999999999999999999998957756146e-01L);
> +  TEST_f_f (tgamma, 0x1.000000000000000000000000008p0L, 9.9999999999999999999999999999998577053525e-01L);
> +  TEST_f_f (tgamma, 0x0.ffffffffffffffffffffffffffcp0L, 1.0000000000000000000000000000000071147324e+00L);
> +  TEST_f_f (tgamma, -0x0.ffffffffffffffffffffffffffcp0L, -8.1129638414606681695789005144064422784335e+31L);
> +  TEST_f_f (tgamma, -0x1.000000000000000000000000008p0L, 4.0564819207303340847894502572031577215665e+31L);
> +  TEST_f_f (tgamma, -0x1.ffffffffffffffffffffffffff8p0L, 2.0282409603651670423947251286016461392168e+31L);
> +  TEST_f_f (tgamma, -0x2.00000000000000000000000001p0L, -1.0141204801825835211973625643007538607832e+31L);
> +  TEST_f_f (tgamma, -0x2.ffffffffffffffffffffffffffp0L, -3.3804016006086117373245418810028760196114e+30L);
> +  TEST_f_f (tgamma, -0x3.00000000000000000000000001p0L, 3.3804016006086117373245418810024573137219e+30L);
> +  TEST_f_f (tgamma, -0x3.ffffffffffffffffffffffffffp0L, 8.4510040015215293433113547025072942156952e+29L);
> +  TEST_f_f (tgamma, -0x4.00000000000000000000000002p0L, -4.2255020007607646716556773512527057843048e+29L);
> +  TEST_f_f (tgamma, -0x4.fffffffffffffffffffffffffep0L, -8.4510040015215293433113547025080884313904e+28L);
> +  TEST_f_f (tgamma, -0x5.00000000000000000000000002p0L, 8.4510040015215293433113547025052449019430e+28L);
> +  TEST_f_f (tgamma, -0x5.fffffffffffffffffffffffffep0L, 1.4085006669202548905518924504180378867132e+28L);
> +  TEST_f_f (tgamma, -0x6.00000000000000000000000002p0L, -1.4085006669202548905518924504175176688423e+28L);
> +  TEST_f_f (tgamma, -0x6.fffffffffffffffffffffffffep0L, -2.0121438098860784150741320720257967542615e+27L);
> +  TEST_f_f (tgamma, -0x7.00000000000000000000000002p0L, 2.0121438098860784150741320720249968965322e+27L);
> +  TEST_f_f (tgamma, -0x7.fffffffffffffffffffffffffep0L, 2.5151797623575980188426650900322769448110e+26L);
> +  TEST_f_f (tgamma, -0x8.00000000000000000000000004p0L, -1.2575898811787990094213325450153421028080e+26L);
> +  TEST_f_f (tgamma, -0x9.fffffffffffffffffffffffffcp0L, 1.3973220901986655660237028277960625420495e+24L);
> +  TEST_f_f (tgamma, -0xa.00000000000000000000000004p0L, -1.3973220901986655660237028277947663821128e+24L);
> +  TEST_f_f (tgamma, -0x13.fffffffffffffffffffffffff8p0L, 1.0420893204640670202556853709074896123293e+12L);
> +  TEST_f_f (tgamma, -0x14.00000000000000000000000008p0L, -1.0420893204640670202556853709050065497299e+12L);
> +  TEST_f_f (tgamma, -0x1d.fffffffffffffffffffffffff8p0L, 9.5580541610429641982963434151488827190079e-03L);
> +  TEST_f_f (tgamma, -0x1e.00000000000000000000000008p0L, -9.5580541610429641982963434151231128067179e-03L);
> +  TEST_f_f (tgamma, -0x27.fffffffffffffffffffffffffp0L, 1.5536546823612837336494536911280147806523e-18L);
> +  TEST_f_f (tgamma, -0x28.0000000000000000000000001p0L, -1.5536546823612837336494536911189419578973e-18L);
> +  TEST_f_f (tgamma, -0x28.fffffffffffffffffffffffffp0L, -3.7894016642958139845108626612879138384405e-20L);
> +  TEST_f_f (tgamma, -0x29.0000000000000000000000001p0L, 3.7894016642958139845108626612656391824122e-20L);
> +  TEST_f_f (tgamma, -0x29.fffffffffffffffffffffffffp0L, 9.0223849149900332964544349078285357440663e-22L);
> +  TEST_f_f (tgamma, -0x2a.0000000000000000000000001p0L, -9.0223849149900332964544349077751619246306e-22L);
> +  TEST_f_f (tgamma, -0x31.fffffffffffffffffffffffffp0L, 4.1679710515150795310771069868348482819424e-35L);
> +  TEST_f_f (tgamma, -0x32.0000000000000000000000001p0L, -4.1679710515150795310771069868090576746248e-35L);
> +  TEST_f_f (tgamma, -0x63.ffffffffffffffffffffffffep0L, 6.7915032994648558610510614163560656864280e-129L);
> +  TEST_f_f (tgamma, -0x64.0000000000000000000000002p0L, -6.7915032994648558610510614162572689693253e-129L);
> +  TEST_f_f (tgamma, -0x95.ffffffffffffffffffffffffcp0L, 5.5468467106873795353190582463444660395353e-234L);
> +  TEST_f_f (tgamma, -0x96.0000000000000000000000004p0L, -5.5468467106873795353190582461689495693675e-234L);
> +  TEST_f_f (tgamma, -0xb4.ffffffffffffffffffffffffcp0L, -8.7154451775644399729745472278589884205029e-303L, UNDERFLOW_EXCEPTION_LDOUBLE_IBM);
> +  TEST_f_f (tgamma, -0xb5.0000000000000000000000004p0L, 8.7154451775644399729745472275729078899858e-303L, UNDERFLOW_EXCEPTION_LDOUBLE_IBM);
> +  TEST_f_f (tgamma, -0xb5.ffffffffffffffffffffffffcp0L, 4.7887061415189230620739270482742524864974e-305L, UNDERFLOW_EXCEPTION_LDOUBLE_IBM);
> +  TEST_f_f (tgamma, -0xb6.0000000000000000000000004p0L, -4.7887061415189230620739270481168993324524e-305L, UNDERFLOW_EXCEPTION_LDOUBLE_IBM);
> +  TEST_f_f (tgamma, -0xb6.ffffffffffffffffffffffffcp0L, -2.6167793123054224382917634143575601407067e-307L, UNDERFLOW_EXCEPTION_LDOUBLE_IBM);
> +  TEST_f_f (tgamma, -0xb7.0000000000000000000000004p0L, 2.6167793123054224382917634142714845691019e-307L, UNDERFLOW_EXCEPTION_LDOUBLE_IBM);
> +  TEST_f_f (tgamma, -0xb7.ffffffffffffffffffffffffcp0L, 1.4221626697312078468976975078030462044826e-309L, UNDERFLOW_EXCEPTION_DOUBLE);
> +  TEST_f_f (tgamma, -0xb8.0000000000000000000000004p0L, -1.4221626697312078468976975077562172247612e-309L, UNDERFLOW_EXCEPTION_DOUBLE);
> +  TEST_f_f (tgamma, -0xbb.ffffffffffffffffffffffffcp0L, 1.1756150745511026776007338998283058108448e-318L, UNDERFLOW_EXCEPTION_DOUBLE);
> +  TEST_f_f (tgamma, -0xbc.0000000000000000000000004p0L, -1.1756150745511026776007338997894360215107e-318L, UNDERFLOW_EXCEPTION_DOUBLE);
> +  TEST_f_f (tgamma, -0xbc.ffffffffffffffffffffffffcp0L, -6.2201855796354639026493857133773832713743e-321L, UNDERFLOW_EXCEPTION_DOUBLE);
> +  TEST_f_f (tgamma, -0xbd.0000000000000000000000004p0L, 6.2201855796354639026493857131715153125175e-321L, UNDERFLOW_EXCEPTION_DOUBLE);
> +  TEST_f_f (tgamma, -0xbd.ffffffffffffffffffffffffcp0L, 3.2737818840186652119207293228302560914295e-323L, UNDERFLOW_EXCEPTION_DOUBLE);
> +  TEST_f_f (tgamma, -0xbe.0000000000000000000000004p0L, -3.2737818840186652119207293227217957948294e-323L, UNDERFLOW_EXCEPTION_DOUBLE);
> +# if LDBL_MAX_EXP <= 1024
> +  TEST_f_f (tgamma, -0xbe.ffffffffffffffffffffffffcp0L, minus_zero, UNDERFLOW_EXCEPTION_DOUBLE);
> +  TEST_f_f (tgamma, -0xbf.0000000000000000000000004p0L, plus_zero, UNDERFLOW_EXCEPTION_DOUBLE);
> +# else
> +  TEST_f_f (tgamma, -0xbe.ffffffffffffffffffffffffcp0L, -1.7140219288055838805867692789687487834686e-325L);
> +  TEST_f_f (tgamma, -0xbf.0000000000000000000000004p0L, 1.7140219288055838805867692789119066543632e-325L);
> +# endif
> +#endif
> +#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
> +  TEST_f_f (tgamma, 0x8.000000000000001p0L, 5.0400000000000000088113830892247051102283e+03L);
> +  TEST_f_f (tgamma, 0x7.fffffffffffffff8p0L, 5.0399999999999999955943084553876474508520e+03L);
> +  TEST_f_f (tgamma, 0x7.0000000000000008p0L, 7.2000000000000000058477733127664675369681e+02L);
> +  TEST_f_f (tgamma, 0x6.fffffffffffffff8p0L, 7.1999999999999999941522266872335324679893e+02L);
> +  TEST_f_f (tgamma, 0x6.0000000000000008p0L, 1.2000000000000000008878927116622375680433e+02L);
> +  TEST_f_f (tgamma, 0x5.fffffffffffffff8p0L, 1.1999999999999999991121072883377624326546e+02L);
> +  TEST_f_f (tgamma, 0x5.0000000000000008p0L, 2.4000000000000000015676186062072582846211e+01L);
> +  TEST_f_f (tgamma, 0x4.fffffffffffffff8p0L, 2.3999999999999999984323813937927417165027e+01L);
> +  TEST_f_f (tgamma, 0x4.0000000000000008p0L, 6.0000000000000000032685252120268430507939e+00L);
> +  TEST_f_f (tgamma, 0x3.fffffffffffffffcp0L, 5.9999999999999999983657373939865784753909e+00L);
> +  TEST_f_f (tgamma, 0x3.0000000000000004p0L, 2.0000000000000000004001939123397399171482e+00L);
> +  TEST_f_f (tgamma, 0x2.fffffffffffffffcp0L, 1.9999999999999999995998060876602600829690e+00L);
> +  TEST_f_f (tgamma, 0x2.0000000000000004p0L, 1.0000000000000000000916767389213195151634e+00L);
> +  TEST_f_f (tgamma, 0x1.fffffffffffffffep0L, 9.9999999999999999995416163053934024243282e-01L);
> +  TEST_f_f (tgamma, 0x1.0000000000000002p0L, 9.9999999999999999993741815221210931418291e-01L);
> +  TEST_f_f (tgamma, 0x0.ffffffffffffffffp0L, 1.0000000000000000000312909238939453429173e+00L);
> +  TEST_f_f (tgamma, -0x0.ffffffffffffffffp0L, -1.8446744073709551616422784335098467139470e+19L);
> +  TEST_f_f (tgamma, -0x1.0000000000000002p0L, 9.2233720368547758075772156649015328607596e+18L);
> +  TEST_f_f (tgamma, -0x1.fffffffffffffffep0L, 4.6116860184273879044613921675492335697983e+18L);
> +  TEST_f_f (tgamma, -0x2.0000000000000004p0L, -2.3058430092136939515386078324507664305064e+18L);
> +  TEST_f_f (tgamma, -0x2.fffffffffffffffcp0L, -7.6861433640456465087601961140530007887063e+17L);
> +  TEST_f_f (tgamma, -0x3.0000000000000004p0L, 7.6861433640456465045731372192803325462836e+17L);
> +  TEST_f_f (tgamma, -0x3.fffffffffffffffcp0L, 1.9215358410114116272942156951799168638773e+17L);
> +  TEST_f_f (tgamma, -0x4.0000000000000008p0L, -9.6076792050570581270578430482008313684602e+16L);
> +  TEST_f_f (tgamma, -0x4.fffffffffffffff8p0L, -1.9215358410114116280884313903598337283601e+16L);
> +  TEST_f_f (tgamma, -0x5.0000000000000008p0L, 1.9215358410114116252449019429734996071487e+16L);
> +  TEST_f_f (tgamma, -0x5.fffffffffffffff8p0L, 3.2025597350190193803788671320812043622696e+15L);
> +  TEST_f_f (tgamma, -0x6.0000000000000008p0L, -3.2025597350190193751766884234743511972877e+15L);
> +  TEST_f_f (tgamma, -0x6.fffffffffffffff8p0L, -4.5750853357414562579675426149912896787735e+14L);
> +  TEST_f_f (tgamma, -0x7.0000000000000008p0L, 4.5750853357414562499689653215166468353753e+14L);
> +  TEST_f_f (tgamma, -0x7.fffffffffffffff8p0L, 5.7188566696768203227694481100089533685959e+13L);
> +  TEST_f_f (tgamma, -0x8.000000000000001p0L, -2.8594283348384101534210280804672371201060e+13L);
> +  TEST_f_f (tgamma, -0x9.fffffffffffffffp0L, 3.1771425942649001828476427167843945971988e+11L);
> +  TEST_f_f (tgamma, -0xa.000000000000001p0L, -3.1771425942649001698860433502350057763905e+11L);
> +  TEST_f_f (tgamma, -0x13.ffffffffffffffep0L, 2.3694367893405502075347562184931828448654e-01L);
> +  TEST_f_f (tgamma, -0x14.000000000000002p0L, -2.3694367893405501827041302245929389013031e-01L);
> +  TEST_f_f (tgamma, -0x1d.ffffffffffffffep0L, 2.1732499046818166459536268654187775086902e-15L);
> +  TEST_f_f (tgamma, -0x1e.000000000000002p0L, -2.1732499046818166201837145753965837196590e-15L);
> +  TEST_f_f (tgamma, -0x27.ffffffffffffffcp0L, 3.5326017549807232935581894777156474496719e-31L);
> +  TEST_f_f (tgamma, -0x28.000000000000004p0L, -3.5326017549807232028299619269080956965758e-31L);
> +  TEST_f_f (tgamma, -0x28.ffffffffffffffcp0L, -8.6161018414163982777002940498289948893044e-33L);
> +  TEST_f_f (tgamma, -0x29.000000000000004p0L, 8.6161018414163980549537337663264762179535e-33L);
> +  TEST_f_f (tgamma, -0x29.ffffffffffffffcp0L, 2.0514528193848567329552463626090806737389e-34L);
> +  TEST_f_f (tgamma, -0x2a.000000000000004p0L, -2.0514528193848566795814269269517457847791e-34L);
> +  TEST_f_f (tgamma, -0x31.ffffffffffffffcp0L, 9.4768689712397635680446279661359728835046e-48L);
> +  TEST_f_f (tgamma, -0x32.000000000000004p0L, -9.4768689712397633101385547903658075308777e-48L);
> +  TEST_f_f (tgamma, -0x63.ffffffffffffff8p0L, 1.5442090669841618542494279375256856430049e-141L);
> +  TEST_f_f (tgamma, -0x64.000000000000008p0L, -1.5442090669841617554527108348771968070612e-141L);
> +  TEST_f_f (tgamma, -0x95.ffffffffffffffp0L, 1.2612069237291916400144732227892704713839e-246L);
> +  TEST_f_f (tgamma, -0x96.00000000000001p0L, -1.2612069237291914644980030550324645611752e-246L);
> +  TEST_f_f (tgamma, -0xb4.ffffffffffffffp0L, -1.9816628031468191243385005680879281767694e-315L);
> +  TEST_f_f (tgamma, -0xb5.00000000000001p0L, 1.9816628031468188382579700510291588022368e-315L);
> +  TEST_f_f (tgamma, -0xb5.ffffffffffffffp0L, 1.0888257160147357826865964233809723297472e-317L);
> +  TEST_f_f (tgamma, -0xb6.00000000000001p0L, -1.0888257160147356253334423783317128355514e-317L);
> +  TEST_f_f (tgamma, -0xb6.ffffffffffffffp0L, -5.9498673006269714905418984659220067091260e-320L);
> +  TEST_f_f (tgamma, -0xb7.00000000000001p0L, 5.9498673006269706297861824177538685763601e-320L);
> +  TEST_f_f (tgamma, -0xb7.ffffffffffffffp0L, 3.2336235329494410277123118903958061569834e-322L);
> +  TEST_f_f (tgamma, -0xb8.00000000000001p0L, -3.2336235329494405594225146768193434900135e-322L);
> +  TEST_f_f (tgamma, -0xbb.ffffffffffffffp0L, 2.6730392040715350119087465463119939092815e-331L);
> +  TEST_f_f (tgamma, -0xbc.00000000000001p0L, -2.6730392040715346232108532050343031951651e-331L);
> +  TEST_f_f (tgamma, -0xbc.ffffffffffffffp0L, -1.4143064571807063556111222197839950086445e-333L);
> +  TEST_f_f (tgamma, -0xbd.00000000000001p0L, 1.4143064571807061497431633629389135273431e-333L);
> +  TEST_f_f (tgamma, -0xbd.ffffffffffffffp0L, 7.4437181956879281879706555863416819210399e-336L);
> +  TEST_f_f (tgamma, -0xbe.00000000000001p0L, -7.4437181956879271033676895858841525581153e-336L);
> +  TEST_f_f (tgamma, -0xbe.ffffffffffffffp0L, -3.8972346574282346536709453101948570578636e-338L);
> +  TEST_f_f (tgamma, -0xbf.00000000000001p0L, 3.8972346574282340852496542564155275274974e-338L);
> +  TEST_f_f (tgamma, -0xf9.ffffffffffffffp0L, 2.2289142548411573883553287678043297937797e-476L);
> +  TEST_f_f (tgamma, -0xfa.00000000000001p0L, -2.2289142548411570466476165308364665814265e-476L);
> +  TEST_f_f (tgamma, -0x1f3.fffffffffffffep0L, 2.9528489142763141594943668922610562530068e-1118L);
> +  TEST_f_f (tgamma, -0x1f4.00000000000002p0L, -2.9528489142763131406565394149878256133744e-1118L);
> +  TEST_f_f (tgamma, -0x2ed.fffffffffffffcp0L, 6.9801511765871818502006905472380418430269e-1817L);
> +  TEST_f_f (tgamma, -0x2ee.00000000000004p0L, -6.9801511765871767194421856376592926002995e-1817L);
> +  TEST_f_f (tgamma, -0x3e7.fffffffffffffcp0L, 4.4768809295877296071892611539415773519036e-2552L);
> +  TEST_f_f (tgamma, -0x3e8.00000000000004p0L, -4.4768809295877261735541135972060089530309e-2552L);
> +  TEST_f_f (tgamma, -0x4e1.fffffffffffff8p0L, 5.4651488569236507565341414077911790252381e-3315L);
> +  TEST_f_f (tgamma, -0x4e2.00000000000008p0L, -5.4651488569236421026544487194247355967789e-3315L);
> +  TEST_f_f (tgamma, -0x5db.fffffffffffff8p0L, 1.8718211510339187689122114747834510481993e-4099L);
> +  TEST_f_f (tgamma, -0x5dc.00000000000008p0L, -1.8718211510339157291960718369454861898499e-4099L);
> +  TEST_f_f (tgamma, -0x6d5.fffffffffffff8p0L, 4.2925786447266492555651378780094011518063e-4902L);
> +  TEST_f_f (tgamma, -0x6d6.00000000000008p0L, -4.2925786447266421378134368786479937285900e-4902L);
> +  TEST_f_f (tgamma, -0x6e2.fffffffffffff8p0L, -2.8229173528168668283609231628661510187536e-4944L, UNDERFLOW_EXCEPTION);
> +  TEST_f_f (tgamma, -0x6e3.00000000000008p0L, 2.8229173528168621428945467532322397890424e-4944L, UNDERFLOW_EXCEPTION);
> +  TEST_f_f (tgamma, -0x6e3.fffffffffffff8p0L, 1.6002932839097884515524894602387278842659e-4947L, UNDERFLOW_EXCEPTION);
> +  TEST_f_f (tgamma, -0x6e4.00000000000008p0L, -1.6002932839097857951909742110188634438017e-4947L, UNDERFLOW_EXCEPTION);
> +  TEST_f_f (tgamma, -0x6e4.fffffffffffff8p0L, -9.0668174725767051085164382743343408095146e-4951L, UNDERFLOW_EXCEPTION);
> +  TEST_f_f (tgamma, -0x6e5.00000000000008p0L, 9.0668174725766900571689083050287838090189e-4951L, UNDERFLOW_EXCEPTION);
> +# if LDBL_MANT_DIG <= 64
> +  TEST_f_f (tgamma, -0x6e5.fffffffffffff8p0L, plus_zero, UNDERFLOW_EXCEPTION);
> +  TEST_f_f (tgamma, -0x6e6.00000000000008p0L, minus_zero, UNDERFLOW_EXCEPTION);
> +# else
> +  TEST_f_f (tgamma, -0x6e5.fffffffffffff8p0L, 5.1340982290921319983501912122292601780326e-4954L, UNDERFLOW_EXCEPTION);
> +  TEST_f_f (tgamma, -0x6e6.00000000000008p0L, -5.1340982290921234748578193083614111497314e-4954L, UNDERFLOW_EXCEPTION);
> +# endif
> +# if LDBL_MANT_DIG >= 113
> +  TEST_f_f (tgamma, 0x8.0000000000000000000000000008p0L, 5.0400000000000000000000000000000156521606e+03L);
> +  TEST_f_f (tgamma, 0x7.fffffffffffffffffffffffffffcp0L, 5.0399999999999999999999999999999921739197e+03L);
> +  TEST_f_f (tgamma, 0x7.0000000000000000000000000004p0L, 7.2000000000000000000000000000000103877321e+02L);
> +  TEST_f_f (tgamma, 0x6.fffffffffffffffffffffffffffcp0L, 7.1999999999999999999999999999999896122679e+02L);
> +  TEST_f_f (tgamma, 0x6.0000000000000000000000000004p0L, 1.2000000000000000000000000000000015772143e+02L);
> +  TEST_f_f (tgamma, 0x5.fffffffffffffffffffffffffffcp0L, 1.1999999999999999999999999999999984227857e+02L);
> +  TEST_f_f (tgamma, 0x5.0000000000000000000000000004p0L, 2.4000000000000000000000000000000027846500e+01L);
> +  TEST_f_f (tgamma, 0x4.fffffffffffffffffffffffffffcp0L, 2.3999999999999999999999999999999972153500e+01L);
> +  TEST_f_f (tgamma, 0x4.0000000000000000000000000004p0L, 6.0000000000000000000000000000000058060671e+00L);
> +  TEST_f_f (tgamma, 0x3.fffffffffffffffffffffffffffep0L, 5.9999999999999999999999999999999970969664e+00L);
> +  TEST_f_f (tgamma, 0x3.0000000000000000000000000002p0L, 2.0000000000000000000000000000000007108872e+00L);
> +  TEST_f_f (tgamma, 0x2.fffffffffffffffffffffffffffep0L, 1.9999999999999999999999999999999992891128e+00L);
> +  TEST_f_f (tgamma, 0x2.0000000000000000000000000002p0L, 1.0000000000000000000000000000000001628506e+00L);
> +  TEST_f_f (tgamma, 0x1.ffffffffffffffffffffffffffffp0L, 9.9999999999999999999999999999999991857470e-01L);
> +  TEST_f_f (tgamma, 0x1.0000000000000000000000000001p0L, 9.9999999999999999999999999999999988883231e-01L);
> +  TEST_f_f (tgamma, 0x0.ffffffffffffffffffffffffffff8p0L, 1.0000000000000000000000000000000000555838e+00L);
> +  TEST_f_f (tgamma, -0x0.ffffffffffffffffffffffffffff8p0L, -1.0384593717069655257060992658440192422784e+34L);
> +  TEST_f_f (tgamma, -0x1.0000000000000000000000000001p0L, 5.1922968585348276285304963292200955772157e+33L);
> +  TEST_f_f (tgamma, -0x1.ffffffffffffffffffffffffffffp0L, 2.5961484292674138142652481646100484613922e+33L);
> +  TEST_f_f (tgamma, -0x2.0000000000000000000000000002p0L, -1.2980742146337069071326240823050235386078e+33L);
> +  TEST_f_f (tgamma, -0x2.fffffffffffffffffffffffffffep0L, -4.3269140487790230237754136076834154268628e+32L);
> +  TEST_f_f (tgamma, -0x3.0000000000000000000000000002p0L, 4.3269140487790230237754136076834112398039e+32L);
> +  TEST_f_f (tgamma, -0x3.fffffffffffffffffffffffffffep0L, 1.0817285121947557559438534019208539608824e+32L);
> +  TEST_f_f (tgamma, -0x4.0000000000000000000000000004p0L, -5.4086425609737787797192670096042603911764e+31L);
> +  TEST_f_f (tgamma, -0x4.fffffffffffffffffffffffffffcp0L, -1.0817285121947557559438534019208547550981e+31L);
> +  TEST_f_f (tgamma, -0x5.0000000000000000000000000004p0L, 1.0817285121947557559438534019208519115686e+31L);
> +  TEST_f_f (tgamma, -0x5.fffffffffffffffffffffffffffcp0L, 1.8028808536579262599064223365347581566449e+30L);
> +  TEST_f_f (tgamma, -0x6.0000000000000000000000000004p0L, -1.8028808536579262599064223365347529544662e+30L);
> +  TEST_f_f (tgamma, -0x6.fffffffffffffffffffffffffffcp0L, -2.5755440766541803712948890521925119357966e+29L);
> +  TEST_f_f (tgamma, -0x7.0000000000000000000000000004p0L, 2.5755440766541803712948890521925039372193e+29L);
> +  TEST_f_f (tgamma, -0x7.fffffffffffffffffffffffffffcp0L, 3.2194300958177254641186113152406402297656e+28L);
> +  TEST_f_f (tgamma, -0x8.0000000000000000000000000008p0L, -1.6097150479088627320593056576203121511868e+28L);
> +  TEST_f_f (tgamma, -0x9.fffffffffffffffffffffffffff8p0L, 1.7885722754542919245103396195781369922635e+26L);
> +  TEST_f_f (tgamma, -0xa.0000000000000000000000000008p0L, -1.7885722754542919245103396195781240306642e+26L);
> +  TEST_f_f (tgamma, -0x13.fffffffffffffffffffffffffffp0L, 1.3338743301940057859272772747600099590309e+14L);
> +  TEST_f_f (tgamma, -0x14.000000000000000000000000001p0L, -1.3338743301940057859272772747599851284049e+14L);
> +  TEST_f_f (tgamma, -0x1d.fffffffffffffffffffffffffffp0L, 1.2234309326134994173819319571374205986026e+00L);
> +  TEST_f_f (tgamma, -0x1e.000000000000000000000000001p0L, -1.2234309326134994173819319571373948286903e+00L);
> +  TEST_f_f (tgamma, -0x27.ffffffffffffffffffffffffffep0L, 1.9886779934224431790713007246380976767855e-16L);
> +  TEST_f_f (tgamma, -0x28.000000000000000000000000002p0L, -1.9886779934224431790713007246380069485580e-16L);
> +  TEST_f_f (tgamma, -0x28.ffffffffffffffffffffffffffep0L, -4.8504341302986419001739042064343853066259e-18L);
> +  TEST_f_f (tgamma, -0x29.000000000000000000000000002p0L, 4.8504341302986419001739042064341625600656e-18L);
> +  TEST_f_f (tgamma, -0x29.ffffffffffffffffffffffffffep0L, 1.1548652691187242619461676681986633377063e-19L);
> +  TEST_f_f (tgamma, -0x2a.000000000000000000000000002p0L, -1.1548652691187242619461676681986099638869e-19L);
> +  TEST_f_f (tgamma, -0x31.ffffffffffffffffffffffffffep0L, 5.3350029459393017997786969431322287652396e-33L);
> +  TEST_f_f (tgamma, -0x32.000000000000000000000000002p0L, -5.3350029459393017997786969431319708591664e-33L);
> +  TEST_f_f (tgamma, -0x63.ffffffffffffffffffffffffffcp0L, 8.6931242233150155021453586128730281632676e-127L);
> +  TEST_f_f (tgamma, -0x64.000000000000000000000000004p0L, -8.6931242233150155021453586128720401960966e-127L);
> +  TEST_f_f (tgamma, -0x95.ffffffffffffffffffffffffff8p0L, 7.0999637896798458052083945552094635720486e-232L);
> +  TEST_f_f (tgamma, -0x96.000000000000000000000000008p0L, -7.0999637896798458052083945552077084073470e-232L);
> +  TEST_f_f (tgamma, -0xb4.ffffffffffffffffffffffffff8p0L, -1.1155769827282483165407420451477844041365e-300L);
> +  TEST_f_f (tgamma, -0xb5.000000000000000000000000008p0L, 1.1155769827282483165407420451474983236060e-300L);
> +  TEST_f_f (tgamma, -0xb5.ffffffffffffffffffffffffff8p0L, 6.1295438611442215194546266216911239298981e-303L);
> +  TEST_f_f (tgamma, -0xb6.000000000000000000000000008p0L, -6.1295438611442215194546266216895503983577e-303L);
> +  TEST_f_f (tgamma, -0xb6.ffffffffffffffffffffffffff8p0L, -3.3494775197509407210134571703230189921356e-305L);
> +  TEST_f_f (tgamma, -0xb7.000000000000000000000000008p0L, 3.3494775197509407210134571703221582364195e-305L);
> +  TEST_f_f (tgamma, -0xb7.ffffffffffffffffffffffffff8p0L, 1.8203682172559460440290528099581627396147e-307L);
> +  TEST_f_f (tgamma, -0xb8.000000000000000000000000008p0L, -1.8203682172559460440290528099576944498174e-307L);
> +  TEST_f_f (tgamma, -0xbb.ffffffffffffffffffffffffff8p0L, 1.5047872954254114273289393917555491216542e-316L);
> +  TEST_f_f (tgamma, -0xbc.000000000000000000000000008p0L, -1.5047872954254114273289393917551604237609e-316L);
> +  TEST_f_f (tgamma, -0xbc.ffffffffffffffffffffffffff8p0L, -7.9618375419333937953912137129923244334851e-319L);
> +  TEST_f_f (tgamma, -0xbd.000000000000000000000000008p0L, 7.9618375419333937953912137129902657538965e-319L);
> +  TEST_f_f (tgamma, -0xbd.ffffffffffffffffffffffffff8p0L, 4.1904408115438914712585335331538555086887e-321L);
> +  TEST_f_f (tgamma, -0xbe.000000000000000000000000008p0L, -4.1904408115438914712585335331527709057227e-321L);
> +  TEST_f_f (tgamma, -0xbe.ffffffffffffffffffffffffff8p0L, -2.1939480688711473671510646770439036908579e-323L);
> +  TEST_f_f (tgamma, -0xbf.000000000000000000000000008p0L, 2.1939480688711473671510646770433352695669e-323L);
> +  TEST_f_f (tgamma, -0xf9.ffffffffffffffffffffffffff8p0L, 1.2547671759429278005937024349858508569625e-461L);
> +  TEST_f_f (tgamma, -0xfa.000000000000000000000000008p0L, -1.2547671759429278005937024349855091492502e-461L);
> +  TEST_f_f (tgamma, -0x1f3.ffffffffffffffffffffffffffp0L, 1.6623061587520224800948170558276526986560e-1103L);
> +  TEST_f_f (tgamma, -0x1f4.00000000000000000000000001p0L, -1.6623061587520224800948170558266338608285e-1103L);
> +  TEST_f_f (tgamma, -0x2ed.fffffffffffffffffffffffffep0L, 3.9294757797334687313030998663853566763041e-1802L);
> +  TEST_f_f (tgamma, -0x2ee.00000000000000000000000002p0L, -3.9294757797334687313030998663802259177992e-1802L);
> +  TEST_f_f (tgamma, -0x3e7.fffffffffffffffffffffffffep0L, 2.5202599107841713834679953735597552855978e-2537L);
> +  TEST_f_f (tgamma, -0x3e8.00000000000000000000000002p0L, -2.5202599107841713834679953735563216504503e-2537L);
> +  TEST_f_f (tgamma, -0x4e1.fffffffffffffffffffffffffcp0L, 3.0766052944457032773369415396747635242705e-3300L);
> +  TEST_f_f (tgamma, -0x4e2.00000000000000000000000004p0L, -3.0766052944457032773369415396661096445778e-3300L);
> +  TEST_f_f (tgamma, -0x5db.fffffffffffffffffffffffffcp0L, 1.0537416297875703295453200836588944487917e-4084L);
> +  TEST_f_f (tgamma, -0x5dc.00000000000000000000000004p0L, -1.0537416297875703295453200836558547326520e-4084L);
> +  TEST_f_f (tgamma, -0x6d5.fffffffffffffffffffffffffcp0L, 2.4165069481061837867659591369998698607387e-4887L);
> +  TEST_f_f (tgamma, -0x6d6.00000000000000000000000004p0L, -2.4165069481061837867659591369927521090377e-4887L);
> +  TEST_f_f (tgamma, -0x6e2.fffffffffffffffffffffffffcp0L, -1.5891611922804672355414813803008070609431e-4929L);
> +  TEST_f_f (tgamma, -0x6e3.00000000000000000000000004p0L, 1.5891611922804672355414813802961215945667e-4929L);
> +  TEST_f_f (tgamma, -0x6e3.fffffffffffffffffffffffffcp0L, 9.0088502963745308137272187091882496519917e-4933L, UNDERFLOW_EXCEPTION);
> +  TEST_f_f (tgamma, -0x6e4.00000000000000000000000004p0L, -9.0088502963745308137272187091616860368392e-4933L, UNDERFLOW_EXCEPTION);
> +  TEST_f_f (tgamma, -0x6e4.fffffffffffffffffffffffffcp0L, -5.1041644738665896961627301468488672286746e-4936L, UNDERFLOW_EXCEPTION);
> +  TEST_f_f (tgamma, -0x6e5.00000000000000000000000004p0L, 5.1041644738665896961627301468338158811446e-4936L, UNDERFLOW_EXCEPTION);
> +  TEST_f_f (tgamma, -0x6e5.fffffffffffffffffffffffffcp0L, 2.8902403589278537350864836618623260468143e-4939L, UNDERFLOW_EXCEPTION);
> +  TEST_f_f (tgamma, -0x6e6.00000000000000000000000004p0L, -2.8902403589278537350864836618538025544424e-4939L, UNDERFLOW_EXCEPTION);
> +  TEST_f_f (tgamma, -0x6eb.fffffffffffffffffffffffffcp0L, 9.4152481960798287090819640211017219042636e-4959L, UNDERFLOW_EXCEPTION);
> +  TEST_f_f (tgamma, -0x6ec.00000000000000000000000004p0L, -9.4152481960798287090819640210739431797168e-4959L, UNDERFLOW_EXCEPTION);
> +  TEST_f_f (tgamma, -0x6ec.fffffffffffffffffffffffffcp0L, -5.3103486723518492437010513373388172315540e-4962L, UNDERFLOW_EXCEPTION);
> +  TEST_f_f (tgamma, -0x6ed.00000000000000000000000004p0L, 5.3103486723518492437010513373231484108489e-4962L, UNDERFLOW_EXCEPTION);
> +  TEST_f_f (tgamma, -0x6ed.fffffffffffffffffffffffffcp0L, 2.9934321715624854812294539669328172615019e-4965L, UNDERFLOW_EXCEPTION);
> +  TEST_f_f (tgamma, -0x6ee.00000000000000000000000004p0L, -2.9934321715624854812294539669239841152754e-4965L, UNDERFLOW_EXCEPTION);
> +  TEST_f_f (tgamma, -0x6ee.fffffffffffffffffffffffffcp0L, minus_zero, UNDERFLOW_EXCEPTION);
> +  TEST_f_f (tgamma, -0x6ef.00000000000000000000000004p0L, plus_zero, UNDERFLOW_EXCEPTION);
> +# endif
> +#endif
> +  TEST_f_f (tgamma, -0x1.0a32a2p+5L, 1.8125267978155035272941154746083439329912e-37L);
> +#ifndef TEST_FLOAT
> +  TEST_f_f (tgamma, -0x1.5800000080001p+7L, -3.1439271448823567326093363350637118195240e-304L, UNDERFLOW_EXCEPTION_LDOUBLE_IBM);
> +#endif
> +  TEST_f_f (tgamma, 18.5L, 1.4986120533153361177371791123515513270334e+15L);
> +  TEST_f_f (tgamma, 19.5L, 2.7724322986333718178137813578503699550119e+16L);
> +  TEST_f_f (tgamma, 23.5L, 5.3613035875444147334274983856108155717836e+21L);
> +  TEST_f_f (tgamma, 29.5L, 1.6348125198274266444378807806868221866931e+30L);
> +  TEST_f_f (tgamma, 30.5L, 4.8226969334909086010917483030261254507447e+31L);
> +  TEST_f_f (tgamma, 31.5L, 1.4709225647147271233329832324229682624771e+33L);
> +  TEST_f_f (tgamma, 32.5L, 4.6334060788513904384988971821323500268029e+34L);
> +  TEST_f_f (tgamma, 33.5L, 1.5058569756267018925121415841930137587110e+36L);
> +  TEST_f_f (tgamma, 34.5L, 5.0446208683494513399156743070465960916817e+37L);
> +  TEST_f_f (tgamma, 0x2.30a43cp+4L, 3.4027979115654976101247558405326779640190e+38L);
> +#ifdef TEST_FLOAT
> +  TEST_f_f (tgamma, 0x2.30a44p+4L, plus_infty, OVERFLOW_EXCEPTION);
> +#else
> +  TEST_f_f (tgamma, 0x2.30a44p+4L, 3.4028438913396451054667218138127983367181e+38L);
> +#endif
> +#ifndef TEST_FLOAT
> +  TEST_f_f (tgamma, 0xa.b9fd72b0fb238p+4L, 1.7976931348622298700886249281842651388250e+308L);
> +# if !defined TEST_LDOUBLE || LDBL_MAX_EXP <= 1024
> +  TEST_f_f (tgamma, 0xa.b9fd72b0fb24p+4L, plus_infty, OVERFLOW_EXCEPTION);
> +# else
> +  TEST_f_f (tgamma, 0xa.b9fd72b0fb24p+4L, 1.7976931348624926129589466917917870782484e+308L);
> +# endif
> +#endif
> +#if defined TEST_LDOUBLE && LDBL_MANT_DIG >= 106
> +  TEST_f_f (tgamma, 0xa.b9fd72b0fb23a9ddbf0d3804f4p+4L, 1.7976931348623158079372897140307719334857e+308L);
> +# if LDBL_MAX_EXP <= 1024
> +  TEST_f_f (tgamma, 0xa.b9fd72b0fb23a9ddbf0d3804f8p+4L, plus_infty, OVERFLOW_EXCEPTION);
> +# else
> +  TEST_f_f (tgamma, 0xa.b9fd72b0fb23a9ddbf0d3804f8p+4L, 1.7976931348623158079372897140599422519044e+308L);
> +# endif
> +#endif
> +#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
> +  TEST_f_f (tgamma, 0x6.db8c603359a97108p+8L, 1.1897314953572317517071551278058233700012e+4932L);
> +  TEST_f_f (tgamma, 0x6.db8c603359a9711p+8L, plus_infty, OVERFLOW_EXCEPTION);
> +# if LDBL_MANT_DIG >= 113
> +  TEST_f_f (tgamma, 0x6.db8c603359a971081bc4a2e9dfdp+8L, 1.1897314953572317650857593266265995494998e+4932L);
> +  TEST_f_f (tgamma, 0x6.db8c603359a971081bc4a2e9dfd4p+8L, plus_infty, OVERFLOW_EXCEPTION);
> +# endif
> +#endif
> +
>    END (tgamma);
>  }
>  
> diff --git a/sysdeps/generic/math_private.h b/sysdeps/generic/math_private.h
> index 7661788..9d6ecad 100644
> --- a/sysdeps/generic/math_private.h
> +++ b/sysdeps/generic/math_private.h
> @@ -371,6 +371,18 @@ extern float __x2y2m1f (float x, float y);
>  extern double __x2y2m1 (double x, double y);
>  extern long double __x2y2m1l (long double x, long double y);
>  
> +/* Compute the product of X + X_EPS, X + X_EPS + 1, ..., X + X_EPS + N
> +   - 1, in the form R * (1 + *EPS) where the return value R is an
> +   approximation to the product and *EPS is set to indicate the
> +   approximate error in the return value.  X is such that all the
> +   values X + 1, ..., X + N - 1 are exactly representable, and X_EPS /
> +   X is small enough that factors quadratic in it can be
> +   neglected.  */
> +extern float __gamma_productf (float x, float x_eps, int n, float *eps);
> +extern double __gamma_product (double x, double x_eps, int n, double *eps);
> +extern long double __gamma_productl (long double x, long double x_eps,
> +				     int n, long double *eps);
> +
>  #ifndef math_opt_barrier
>  # define math_opt_barrier(x) \
>  ({ __typeof (x) __x = (x); __asm ("" : "+m" (__x)); __x; })
> diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps
> index 0815592..8761d3a 100644
> --- a/sysdeps/i386/fpu/libm-test-ulps
> +++ b/sysdeps/i386/fpu/libm-test-ulps
> @@ -6165,6 +6165,379 @@ idouble: 2
>  ifloat: 1
>  ildouble: 1
>  ldouble: 1
> +Test "tgamma (-0x0.ffffffffffffffffp0) == -1.8446744073709551616422784335098467139470e+19":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x1.000002p0) == 8.3886075772158332060084424806449513922858e+06":
> +double: 1
> +idouble: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x1.0a32a2p+5) == 1.8125267978155035272941154746083439329912e-37":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x1.5800000080001p+7) == -3.1439271448823567326093363350637118195240e-304":
> +double: 1
> +idouble: 1
> +ildouble: 2
> +ldouble: 2
> +Test "tgamma (-0x1.fffffffffffffp0) == 2.2517998136852484613921675492337776673289e+15":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x13.ffffep0) == 2.1550026214525536756224040483579183652119e-13":
> +float: 1
> +ifloat: 1
> +Test "tgamma (-0x13.ffffffffffffffep0) == 2.3694367893405502075347562184931828448654e-01":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x13.ffffffffffffp0) == 1.1569515572952029402736625857313236848570e-04":
> +double: 1
> +idouble: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x14.000000000001p0) == -1.1569515572951781096476686854873801225397e-04":
> +ildouble: 2
> +ldouble: 2
> +Test "tgamma (-0x1d.ffffep0) == 1.9765721589464867957912772592816027583176e-27":
> +double: 1
> +idouble: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x1d.ffffffffffffffep0) == 2.1732499046818166459536268654187775086902e-15":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x1e.000000000000002p0) == -2.1732499046818166201837145753965837196590e-15":
> +ildouble: 2
> +ldouble: 2
> +Test "tgamma (-0x1e.00002p0) == -1.9765463890341964384070157599286498212650e-27":
> +float: 1
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x1f4.00000000000002p0) == -2.9528489142763131406565394149878256133744e-1118":
> +ildouble: 3
> +ldouble: 3
> +Test "tgamma (-0x1p-24) == -1.6777216577215723853867349114260580375249e+07":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x2.0000000000000004p0) == -2.3058430092136939515386078324507664305064e+18":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x2.0000000000002p0) == -1.1258999068426235386078324507668462444260e+15":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x2.000004p0) == -2.0971515386080557574407223895988378776747e+06":
> +double: 1
> +idouble: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x2.fffffcp0) == -6.9905087601970247876992248591045142913324e+05":
> +double: 1
> +idouble: 1
> +Test "tgamma (-0x27.ffffcp0) == 3.2129279441390812141195076945616975790225e-43":
> +double: 1
> +idouble: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x27.fffffffffffep0) == 1.7249032006742266376460389310340465554361e-34":
> +double: 1
> +idouble: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x28.000000000002p0) == -1.7249032006741359094184881234822934593822e-34":
> +double: 1
> +idouble: 1
> +Test "tgamma (-0x28.00004p0) == -3.2128372159115252365699015758097981155793e-43":
> +double: 1
> +idouble: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x28.ffffcp0) == -7.8364103489619817539676737414096652170685e-45":
> +double: 1
> +idouble: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x28.fffffffffffep0) == -4.2070809772542120404320040128839297118648e-36":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x28.ffffffffffffffcp0) == -8.6161018414163982777002940498289948893044e-33":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x29.000000000000004p0) == 8.6161018414163980549537337663264762179535e-33":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x29.000000000002p0) == 4.2070809772539892938717205103652583609422e-36":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x29.00004p0) == 7.8361876024016854597745353972619195760515e-45":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x29.ffffcp0) == 1.8658121573125798145204120066590953505132e-46":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x2a.00004p0) == -1.8657587834931410688246126853566488626385e-46":
> +double: 1
> +idouble: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x2ed.fffffffffffffcp0) == 6.9801511765871818502006905472380418430269e-1817":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x3.0000000000002p0) == 3.7529996894754112398038859470009084971438e+14":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x3.000004p0) == 6.9905045731381300146131914617735687322025e+05":
> +double: 1
> +idouble: 1
> +Test "tgamma (-0x3.fffffcp0) == 1.7476272942159602684441970627092458855771e+05":
> +float: 1
> +ifloat: 1
> +Test "tgamma (-0x3.ffffffffffffep0) == 9.3824992236885396088236184658402406857503e+13":
> +double: 1
> +idouble: 1
> +Test "tgamma (-0x3.fffffffffffffffcp0) == 1.9215358410114116272942156951799168638773e+17":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x31.fffffffffffep0) == 4.6273774273632946947805289899230181990085e-51":
> +double: 1
> +idouble: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x32.000000000000004p0) == -9.4768689712397633101385547903658075308777e-48":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x32.000000000002p0) == -4.6273774273630367887073532197576655720178e-51":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x3e7.fffffffffffffcp0) == 4.4768809295877296071892611539415773519036e-2552":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x3e8.00000000000004p0) == -4.4768809295877261735541135972060089530309e-2552":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x4.0000000000000008p0) == -9.6076792050570581270578430482008313684602e+16":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x4.000008p0) == -8.7381270578483499672965708923121931082305e+04":
> +float: 2
> +ifloat: 2
> +Test "tgamma (-0x4.fffff8p0) == -1.7476280884325863043793087474680780379554e+04":
> +float: 1
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x4.ffffffffffffcp0) == -9.3824992236885475509805702650262155809819e+12":
> +double: 1
> +idouble: 1
> +Test "tgamma (-0x4e2.00000000000008p0) == -5.4651488569236421026544487194247355967789e-3315":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x5.0000000000000008p0) == 1.9215358410114116252449019429734996071487e+16":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x5.0000000000004p0) == 9.3824992236885191156860964016850034672946e+12":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x5.000008p0) == 1.7476252449031389167286893378510439443844e+04":
> +float: 1
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x5.fffff8p0) == 2.9127137122026653716311560165769071985443e+03":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +Test "tgamma (-0x5.ffffffffffffcp0) == 1.5637498706147581566449098589862357835505e+12":
> +double: 1
> +idouble: 1
> +Test "tgamma (-0x5db.fffffffffffff8p0) == 1.8718211510339187689122114747834510481993e-4099":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x6.0000000000000008p0) == -3.2025597350190193751766884234743511972877e+15":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x6.000008p0) == -2.9127085100239567622341538102130981196910e+03":
> +double: 1
> +idouble: 1
> +Test "tgamma (-0x6.fffff8p0) == -4.1610198723079349791939054365613377035519e+02":
> +double: 2
> +float: 1
> +idouble: 2
> +ifloat: 1
> +ildouble: 2
> +ldouble: 2
> +Test "tgamma (-0x6.ffffffffffffcp0) == -2.2339283865925119357965832452642909859289e+11":
> +double: 3
> +idouble: 3
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x63.fffffffffffcp0) == 7.5400833348840965463348754984345825364294e-145":
> +double: 1
> +idouble: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x63.ffffffffffffff8p0) == 1.5442090669841618542494279375256856430049e-141":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x64.000000000000008p0) == -1.5442090669841617554527108348771968070612e-141":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x64.000000000004p0) == -7.5400833348831085791638490135462230991587e-145":
> +ildouble: 2
> +ldouble: 2
> +Test "tgamma (-0x6d6.00000000000008p0) == -4.2925786447266421378134368786479937285900e-4902":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x7.0000000000000008p0) == 4.5750853357414562499689653215166468353753e+14":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x7.0000000000004p0) == 2.2339283865925039372192897706214475877342e+11":
> +double: 4
> +idouble: 4
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x7.000008p0) == 4.1610118737306415004517215226199741948733e+02":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x7.fffff8p0) == 5.2012751504050764429534086402871289946986e+01":
> +double: 1
> +idouble: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x7.ffffffffffffcp0) == 2.7924104832406402297655703264222230055898e+10":
> +double: 2
> +idouble: 2
> +Test "tgamma (-0x7.fffffffffffffff8p0) == 5.7188566696768203227694481100089533685959e+13":
> +ildouble: 4
> +ldouble: 4
> +Test "tgamma (-0x8.000000000000001p0) == -2.8594283348384101534210280804672371201060e+13":
> +ildouble: 2
> +ldouble: 2
> +Test "tgamma (-0x8.00001p0) == -2.6006296115134418896533598545925084576702e+01":
> +float: 1
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x9.ffffffffffff8p0) == 1.5513391573559147700413058496716749249803e+08":
> +double: 1
> +idouble: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x95.fffffffffff8p0) == 6.1582369322723207086020016423767264008839e-250":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x95.ffffffffffffffp0) == 1.2612069237291916400144732227892704713839e-246":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x96.000000000008p0) == -6.1582369322705655439003240743176243138734e-250":
> +double: 1
> +idouble: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0xa.000000000000001p0) == -3.1771425942649001698860433502350057763905e+11":
> +ildouble: 2
> +ldouble: 2
> +Test "tgamma (-0xa.0000000000008p0) == -1.5513391573559018084419393002828541166901e+08":
> +double: 1
> +idouble: 1
> +Test "tgamma (-0xa.00001p0) == -2.8895878754728051776830454190076999107021e-01":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0xb4.ffffffffffffffp0) == -1.9816628031468191243385005680879281767694e-315":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0xb5.00000000000001p0) == 1.9816628031468188382579700510291588022368e-315":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0xb5.000000000008p0) == 9.6760879059888966544677044221698800670218e-319":
> +ildouble: 2
> +ldouble: 2
> +Test "tgamma (-0xb5.ffffffffffffffp0) == 1.0888257160147357826865964233809723297472e-317":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0xb6.00000000000001p0) == -1.0888257160147356253334423783317128355514e-317":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0xb6.fffffffffff8p0) == -2.9052086428846935908287469917922960610289e-323":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0xb6.ffffffffffffffp0) == -5.9498673006269714905418984659220067091260e-320":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0xb7.fffffffffff8p0) == 1.5789177406982032823826953250736039527543e-325":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0xb7.ffffffffffffffp0) == 3.2336235329494410277123118903958061569834e-322":
> +ildouble: 2
> +ldouble: 2
> +Test "tgamma (-0xb8.000000000008p0) == -1.5789177406977349925854817486109369828857e-325":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0xbb.ffffffffffffffp0) == 2.6730392040715350119087465463119939092815e-331":
> +ildouble: 2
> +ldouble: 2
> +Test "tgamma (-0xbc.00000000000001p0) == -2.6730392040715346232108532050343031951651e-331":
> +ildouble: 3
> +ldouble: 3
> +Test "tgamma (-0xbd.00000000000001p0) == 1.4143064571807061497431633629389135273431e-333":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0xbe.00000000000001p0) == -7.4437181956879271033676895858841525581153e-336":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0xbe.ffffffffffffffp0) == -3.8972346574282346536709453101948570578636e-338":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0xbf.00000000000001p0) == 3.8972346574282340852496542564155275274974e-338":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0xf9.ffffffffffffffp0) == 2.2289142548411573883553287678043297937797e-476":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0xfa.00000000000001p0) == -2.2289142548411570466476165308364665814265e-476":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-1.5) == 2.3632718012073547030642233111215269103967e+00":
> +float: 1
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-2.5) == -9.4530872048294188122568932444861076415869e-01":
> +double: 1
> +idouble: 1
> +Test "tgamma (-4.5) == -6.0019601300504246427027893615784810422774e-02":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-5.5) == 1.0912654781909862986732344293779056440504e-02":
> +float: 1
> +ifloat: 1
> +Test "tgamma (-6.5) == -1.6788699664476712287280529682737009908468e-03":
> +float: 1
> +ifloat: 1
> +Test "tgamma (-7.5) == 2.2384932885968949716374039576982679877958e-04":
> +double: 2
> +idouble: 2
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-8.5) == -2.6335215159963470254557693619979623385833e-05":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-9.5) == 2.7721279115751021320587045915768024616666e-06":
> +ildouble: 2
> +ldouble: 2
>  Test "tgamma (0.5) == sqrt (pi)":
>  float: 1
>  ifloat: 1
> @@ -6173,7 +6546,218 @@ double: 1
>  float: 1
>  idouble: 1
>  ifloat: 1
> +Test "tgamma (0x1.fffffep0) == 9.9999994960018563231526611134590489120697e-01":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x1.fffffffffffffffep0) == 9.9999999999999999995416163053934024243282e-01":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x1p-24) == 1.6777215422784394050795179874582764575261e+07":
> +float: 1
> +ifloat: 1
> +Test "tgamma (0x1p-53) == 9.0071992547409914227843350984672492007618e+15":
> +double: 1
> +idouble: 1
> +Test "tgamma (0x1p-64) == 1.8446744073709551615422784335098467139447e+19":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x2.0000000000002p0) == 1.0000000000000001877539613108624482361963e+00":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x2.000004p0) == 1.0000001007996638509889062631687945799175e+00":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x2.30a43cp+4) == 3.4027979115654976101247558405326779640190e+38":
> +double: 1
> +idouble: 1
> +Test "tgamma (0x2.fffffcp0) == 1.9999995599822108706107786027549565954046e+00":
> +float: 2
> +ifloat: 2
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x2.ffffffffffffep0) == 1.9999999999999991804028675282128956223990e+00":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x3.0000000000002p0) == 2.0000000000000008195971324717875960213536e+00":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x3.fffffcp0) == 5.9999982031095793171233994481968816873643e+00":
> +float: 1
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x3.fffffffffffffffcp0) == 5.9999999999999999983657373939865784753909e+00":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x4.0000000000000008p0) == 6.0000000000000000032685252120268430507939e+00":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x4.000008p0) == 6.0000035937827461765660468073471093546129e+00":
> +float: 1
> +ifloat: 1
> +Test "tgamma (0x4.fffff8p0) == 2.3999982763857938712639837029547357501709e+01":
> +float: 1
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x4.fffffffffffffff8p0) == 2.3999999999999999984323813937927417165027e+01":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x5.0000000000000008p0) == 2.4000000000000000015676186062072582846211e+01":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x5.000008p0) == 2.4000017236155647574166073485628713443799e+01":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x5.fffff8p0) == 1.1999990237520611552119807476573441975106e+02":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x5.ffffffffffffcp0) == 1.1999999999999981815957265157389249327533e+02":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x6.0000000000000008p0) == 1.2000000000000000008878927116622375680433e+02":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x6.000008p0) == 1.2000009762487825358530770343720418162783e+02":
> +float: 1
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x6.fffff8p0) == 7.1999935703082425988147448928288557689866e+02":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x6.ffffffffffffcp0) == 7.1999999999999880237602554542848858572672e+02":
> +double: 3
> +idouble: 3
> +Test "tgamma (0x7.0000000000000008p0) == 7.2000000000000000058477733127664675369681e+02":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x7.0000000000004p0) == 7.2000000000000119762397445457359071259652e+02":
> +double: 4
> +idouble: 4
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x7.000008p0) == 7.2000064296977505705636258629805621178587e+02":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x7.fffff8p0) == 5.0399951558933225045148935487583089307135e+03":
> +double: 1
> +idouble: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x7.ffffffffffffcp0) == 5.0399999999999909771437166339103165198442e+03":
> +double: 2
> +idouble: 2
> +Test "tgamma (0x7.fffffffffffffff8p0) == 5.0399999999999999955943084553876474508520e+03":
> +ildouble: 3
> +ldouble: 3
> +Test "tgamma (0x8.000000000000001p0) == 5.0400000000000000088113830892247051102283e+03":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x8.0000000000008p0) == 5.0400000000000180457125667322294144477136e+03":
> +double: 1
> +idouble: 1
> +Test "tgamma (0x8.00001p0) == 5.0400096882277802019946778420223050233915e+03":
> +double: 1
> +idouble: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0xa.b9fd72b0fb238p+4) == 1.7976931348622298700886249281842651388250e+308":
> +double: 1
> +idouble: 1
> +ildouble: 2
> +ldouble: 2
> +Test "tgamma (10) == 362880":
> +double: 1
> +idouble: 1
> +Test "tgamma (18.5) == 1.4986120533153361177371791123515513270334e+15":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (19.5) == 2.7724322986333718178137813578503699550119e+16":
> +double: 1
> +idouble: 1
> +Test "tgamma (2.5) == 1.3293403881791370204736256125058588870982e+00":
> +float: 1
> +ifloat: 1
> +Test "tgamma (23.5) == 5.3613035875444147334274983856108155717836e+21":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (29.5) == 1.6348125198274266444378807806868221866931e+30":
> +double: 1
> +idouble: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (3) == 2":
> +float: 1
> +ifloat: 1
> +Test "tgamma (3.5) == 3.3233509704478425511840640312646472177454e+00":
> +float: 1
> +ifloat: 1
> +Test "tgamma (30.5) == 4.8226969334909086010917483030261254507447e+31":
> +float: 1
> +ifloat: 1
> +Test "tgamma (31.5) == 1.4709225647147271233329832324229682624771e+33":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +Test "tgamma (32.5) == 4.6334060788513904384988971821323500268029e+34":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (34.5) == 5.0446208683494513399156743070465960916817e+37":
> +ildouble: 1
> +ldouble: 1
>  Test "tgamma (4) == 6":
> +float: 1
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (4.5) == 1.1631728396567448929144224109426265262109e+01":
> +double: 1
> +idouble: 1
> +Test "tgamma (5.5) == 5.2342777784553520181149008492418193679490e+01":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (6.5) == 2.8788527781504436099631954670830006523720e+02":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (7.5) == 1.8712543057977883464760770536039504240418e+03":
> +double: 2
> +float: 1
> +idouble: 2
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (8) == 5040":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (8.5) == 1.4034407293483412598570577902029628180313e+04":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +ildouble: 2
> +ldouble: 2
> +Test "tgamma (9) == 40320":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (9.5) == 1.1929246199460900708784991216725183953266e+05":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
>  ildouble: 1
>  ldouble: 1
>  
> @@ -7103,12 +7687,12 @@ ildouble: 2
>  ldouble: 2
>  
>  Function: "tgamma":
> -double: 2
> -float: 1
> -idouble: 2
> -ifloat: 1
> -ildouble: 1
> -ldouble: 1
> +double: 4
> +float: 2
> +idouble: 4
> +ifloat: 2
> +ildouble: 4
> +ldouble: 4
>  
>  Function: "y0":
>  double: 2
> diff --git a/sysdeps/ieee754/dbl-64/e_gamma_r.c b/sysdeps/ieee754/dbl-64/e_gamma_r.c
> index 9873551..5b17f7b 100644
> --- a/sysdeps/ieee754/dbl-64/e_gamma_r.c
> +++ b/sysdeps/ieee754/dbl-64/e_gamma_r.c
> @@ -19,14 +19,104 @@
>  
>  #include <math.h>
>  #include <math_private.h>
> +#include <float.h>
>  
> +/* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's
> +   approximation to gamma function.  */
> +
> +static const double gamma_coeff[] =
> +  {
> +    0x1.5555555555555p-4,
> +    -0xb.60b60b60b60b8p-12,
> +    0x3.4034034034034p-12,
> +    -0x2.7027027027028p-12,
> +    0x3.72a3c5631fe46p-12,
> +    -0x7.daac36664f1f4p-12,
> +  };
> +
> +#define NCOEFF (sizeof (gamma_coeff) / sizeof (gamma_coeff[0]))
> +
> +/* Return gamma (X), for positive X less than 184, in the form R *
> +   2^(*EXP2_ADJ), where R is the return value and *EXP2_ADJ is set to
> +   avoid overflow or underflow in intermediate calculations.  */
> +
> +static double
> +gamma_positive (double x, int *exp2_adj)
> +{
> +  int local_signgam;
> +  if (x < 0.5)
> +    {
> +      *exp2_adj = 0;
> +      return __ieee754_exp (__ieee754_lgamma_r (x + 1, &local_signgam)) / x;
> +    }
> +  else if (x <= 1.5)
> +    {
> +      *exp2_adj = 0;
> +      return __ieee754_exp (__ieee754_lgamma_r (x, &local_signgam));
> +    }
> +  else if (x < 6.5)
> +    {
> +      /* Adjust into the range for using exp (lgamma).  */
> +      *exp2_adj = 0;
> +      double n = __ceil (x - 1.5);
> +      double x_adj = x - n;
> +      double eps;
> +      double prod = __gamma_product (x_adj, 0, n, &eps);
> +      return (__ieee754_exp (__ieee754_lgamma_r (x_adj, &local_signgam))
> +	      * prod * (1.0 + eps));
> +    }
> +  else
> +    {
> +      double eps = 0;
> +      double x_eps = 0;
> +      double x_adj = x;
> +      double prod = 1;
> +      if (x < 12.0)
> +	{
> +	  /* Adjust into the range for applying Stirling's
> +	     approximation.  */
> +	  double n = __ceil (12.0 - x);
> +#if FLT_EVAL_METHOD != 0
> +	  volatile
> +#endif
> +	  double x_tmp = x + n;
> +	  x_adj = x_tmp;
> +	  x_eps = (x - (x_adj - n));
> +	  prod = __gamma_product (x_adj - n, x_eps, n, &eps);
> +	}
> +      /* The result is now gamma (X_ADJ + X_EPS) / (PROD * (1 + EPS)).
> +	 Compute gamma (X_ADJ + X_EPS) using Stirling's approximation,
> +	 starting by computing pow (X_ADJ, X_ADJ) with a power of 2
> +	 factored out.  */
> +      double exp_adj = -eps;
> +      double x_adj_int = __round (x_adj);
> +      double x_adj_frac = x_adj - x_adj_int;
> +      int x_adj_log2;
> +      double x_adj_mant = __frexp (x_adj, &x_adj_log2);
> +      if (x_adj_mant < M_SQRT1_2)
> +	{
> +	  x_adj_log2--;
> +	  x_adj_mant *= 2.0;
> +	}
> +      *exp2_adj = x_adj_log2 * (int) x_adj_int;
> +      double ret = (__ieee754_pow (x_adj_mant, x_adj)
> +		    * __ieee754_exp2 (x_adj_log2 * x_adj_frac)
> +		    * __ieee754_exp (-x_adj)
> +		    * __ieee754_sqrt (2 * M_PI / x_adj)
> +		    / prod);
> +      exp_adj += x_eps * __ieee754_log (x);
> +      double bsum = gamma_coeff[NCOEFF - 1];
> +      double x_adj2 = x_adj * x_adj;
> +      for (size_t i = 1; i <= NCOEFF - 1; i++)
> +	bsum = bsum / x_adj2 + gamma_coeff[NCOEFF - 1 - i];
> +      exp_adj += bsum / x_adj;
> +      return ret + ret * __expm1 (exp_adj);
> +    }
> +}
>  
>  double
>  __ieee754_gamma_r (double x, int *signgamp)
>  {
> -  /* We don't have a real gamma implementation now.  We'll use lgamma
> -     and the exp function.  But due to the required boundary
> -     conditions we must check some values separately.  */
>    int32_t hx;
>    u_int32_t lx;
>  
> @@ -51,8 +141,48 @@ __ieee754_gamma_r (double x, int *signgamp)
>        *signgamp = 0;
>        return x - x;
>      }
> +  if (__builtin_expect ((hx & 0x7ff00000) == 0x7ff00000, 0))
> +    {
> +      /* Positive infinity (return positive infinity) or NaN (return
> +	 NaN).  */
> +      *signgamp = 0;
> +      return x + x;
> +    }
>  
> -  /* XXX FIXME.  */
> -  return __ieee754_exp (__ieee754_lgamma_r (x, signgamp));
> +  if (x >= 172.0)
> +    {
> +      /* Overflow.  */
> +      *signgamp = 0;
> +      return DBL_MAX * DBL_MAX;
> +    }
> +  else if (x > 0.0)
> +    {
> +      *signgamp = 0;
> +      int exp2_adj;
> +      double ret = gamma_positive (x, &exp2_adj);
> +      return __scalbn (ret, exp2_adj);
> +    }
> +  else if (x >= -DBL_EPSILON / 4.0)
> +    {
> +      *signgamp = 0;
> +      return 1.0 / x;
> +    }
> +  else
> +    {
> +      double tx = __trunc (x);
> +      *signgamp = (tx == 2.0 * __trunc (tx / 2.0)) ? -1 : 1;
> +      if (x <= -184.0)
> +	/* Underflow.  */
> +	return DBL_MIN * DBL_MIN;
> +      double frac = tx - x;
> +      if (frac > 0.5)
> +	frac = 1.0 - frac;
> +      double sinpix = (frac <= 0.25
> +		       ? __sin (M_PI * frac)
> +		       : __cos (M_PI * (0.5 - frac)));
> +      int exp2_adj;
> +      double ret = M_PI / (-x * sinpix * gamma_positive (-x, &exp2_adj));
> +      return __scalbn (ret, -exp2_adj);
> +    }
>  }
>  strong_alias (__ieee754_gamma_r, __gamma_r_finite)
> diff --git a/sysdeps/ieee754/dbl-64/gamma_product.c b/sysdeps/ieee754/dbl-64/gamma_product.c
> new file mode 100644
> index 0000000..2a3fc1a
> --- /dev/null
> +++ b/sysdeps/ieee754/dbl-64/gamma_product.c
> @@ -0,0 +1,75 @@
> +/* Compute a product of X, X+1, ..., with an error estimate.
> +   Copyright (C) 2013 Free Software Foundation, Inc.
> +   This file is part of the GNU C Library.
> +
> +   The GNU C Library is free software; you can redistribute it and/or
> +   modify it under the terms of the GNU Lesser General Public
> +   License as published by the Free Software Foundation; either
> +   version 2.1 of the License, or (at your option) any later version.
> +
> +   The GNU C Library is distributed in the hope that it will be useful,
> +   but WITHOUT ANY WARRANTY; without even the implied warranty of
> +   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
> +   Lesser General Public License for more details.
> +
> +   You should have received a copy of the GNU Lesser General Public
> +   License along with the GNU C Library; if not, see
> +   <http://www.gnu.org/licenses/>.  */
> +
> +#include <math.h>
> +#include <math_private.h>
> +#include <float.h>
> +
> +/* Calculate X * Y exactly and store the result in *HI + *LO.  It is
> +   given that the values are small enough that no overflow occurs and
> +   large enough (or zero) that no underflow occurs.  */
> +
> +static void
> +mul_split (double *hi, double *lo, double x, double y)
> +{
> +#ifdef __FP_FAST_FMA
> +  /* Fast built-in fused multiply-add.  */
> +  *hi = x * y;
> +  *lo = __builtin_fma (x, y, -*hi);
> +#elif defined FP_FAST_FMA
> +  /* Fast library fused multiply-add, compiler before GCC 4.6.  */
> +  *hi = x * y;
> +  *lo = __fma (x, y, -*hi);
> +#else
> +  /* Apply Dekker's algorithm.  */
> +  *hi = x * y;
> +# define C ((1 << (DBL_MANT_DIG + 1) / 2) + 1)
> +  double x1 = x * C;
> +  double y1 = y * C;
> +# undef C
> +  x1 = (x - x1) + x1;
> +  y1 = (y - y1) + y1;
> +  double x2 = x - x1;
> +  double y2 = y - y1;
> +  *lo = (((x1 * y1 - *hi) + x1 * y2) + x2 * y1) + x2 * y2;
> +#endif
> +}
> +
> +/* Compute the product of X + X_EPS, X + X_EPS + 1, ..., X + X_EPS + N
> +   - 1, in the form R * (1 + *EPS) where the return value R is an
> +   approximation to the product and *EPS is set to indicate the
> +   approximate error in the return value.  X is such that all the
> +   values X + 1, ..., X + N - 1 are exactly representable, and X_EPS /
> +   X is small enough that factors quadratic in it can be
> +   neglected.  */
> +
> +double
> +__gamma_product (double x, double x_eps, int n, double *eps)
> +{
> +  SET_RESTORE_ROUND (FE_TONEAREST);
> +  double ret = x;
> +  *eps = x_eps / x;
> +  for (int i = 1; i < n; i++)
> +    {
> +      *eps += x_eps / (x + i);
> +      double lo;
> +      mul_split (&ret, &lo, ret, x + i);
> +      *eps += lo / ret;
> +    }
> +  return ret;
> +}
> diff --git a/sysdeps/ieee754/dbl-64/gamma_productf.c b/sysdeps/ieee754/dbl-64/gamma_productf.c
> new file mode 100644
> index 0000000..46072f1
> --- /dev/null
> +++ b/sysdeps/ieee754/dbl-64/gamma_productf.c
> @@ -0,0 +1,46 @@
> +/* Compute a product of X, X+1, ..., with an error estimate.
> +   Copyright (C) 2013 Free Software Foundation, Inc.
> +   This file is part of the GNU C Library.
> +
> +   The GNU C Library is free software; you can redistribute it and/or
> +   modify it under the terms of the GNU Lesser General Public
> +   License as published by the Free Software Foundation; either
> +   version 2.1 of the License, or (at your option) any later version.
> +
> +   The GNU C Library is distributed in the hope that it will be useful,
> +   but WITHOUT ANY WARRANTY; without even the implied warranty of
> +   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
> +   Lesser General Public License for more details.
> +
> +   You should have received a copy of the GNU Lesser General Public
> +   License along with the GNU C Library; if not, see
> +   <http://www.gnu.org/licenses/>.  */
> +
> +#include <math.h>
> +#include <math_private.h>
> +#include <float.h>
> +
> +/* Compute the product of X + X_EPS, X + X_EPS + 1, ..., X + X_EPS + N
> +   - 1, in the form R * (1 + *EPS) where the return value R is an
> +   approximation to the product and *EPS is set to indicate the
> +   approximate error in the return value.  X is such that all the
> +   values X + 1, ..., X + N - 1 are exactly representable, and X_EPS /
> +   X is small enough that factors quadratic in it can be
> +   neglected.  */
> +
> +float
> +__gamma_productf (float x, float x_eps, int n, float *eps)
> +{
> +  double x_full = (double) x + (double) x_eps;
> +  double ret = x_full;
> +  for (int i = 1; i < n; i++)
> +    ret *= x_full + i;
> +
> +#if FLT_EVAL_METHOD != 0
> +  volatile
> +#endif
> +  float fret = ret;
> +  *eps = (ret - fret) / fret;
> +
> +  return fret;
> +}
> diff --git a/sysdeps/ieee754/flt-32/e_gammaf_r.c b/sysdeps/ieee754/flt-32/e_gammaf_r.c
> index a312957..f58f4c8 100644
> --- a/sysdeps/ieee754/flt-32/e_gammaf_r.c
> +++ b/sysdeps/ieee754/flt-32/e_gammaf_r.c
> @@ -19,14 +19,97 @@
>  
>  #include <math.h>
>  #include <math_private.h>
> +#include <float.h>
>  
> +/* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's
> +   approximation to gamma function.  */
> +
> +static const float gamma_coeff[] =
> +  {
> +    0x1.555556p-4f,
> +    -0xb.60b61p-12f,
> +    0x3.403404p-12f,
> +  };
> +
> +#define NCOEFF (sizeof (gamma_coeff) / sizeof (gamma_coeff[0]))
> +
> +/* Return gamma (X), for positive X less than 42, in the form R *
> +   2^(*EXP2_ADJ), where R is the return value and *EXP2_ADJ is set to
> +   avoid overflow or underflow in intermediate calculations.  */
> +
> +static float
> +gammaf_positive (float x, int *exp2_adj)
> +{
> +  int local_signgam;
> +  if (x < 0.5f)
> +    {
> +      *exp2_adj = 0;
> +      return __ieee754_expf (__ieee754_lgammaf_r (x + 1, &local_signgam)) / x;
> +    }
> +  else if (x <= 1.5f)
> +    {
> +      *exp2_adj = 0;
> +      return __ieee754_expf (__ieee754_lgammaf_r (x, &local_signgam));
> +    }
> +  else if (x < 2.5f)
> +    {
> +      *exp2_adj = 0;
> +      float x_adj = x - 1;
> +      return (__ieee754_expf (__ieee754_lgammaf_r (x_adj, &local_signgam))
> +	      * x_adj);
> +    }
> +  else
> +    {
> +      float eps = 0;
> +      float x_eps = 0;
> +      float x_adj = x;
> +      float prod = 1;
> +      if (x < 4.0f)
> +	{
> +	  /* Adjust into the range for applying Stirling's
> +	     approximation.  */
> +	  float n = __ceilf (4.0f - x);
> +#if FLT_EVAL_METHOD != 0
> +	  volatile
> +#endif
> +	  float x_tmp = x + n;
> +	  x_adj = x_tmp;
> +	  x_eps = (x - (x_adj - n));
> +	  prod = __gamma_productf (x_adj - n, x_eps, n, &eps);
> +	}
> +      /* The result is now gamma (X_ADJ + X_EPS) / (PROD * (1 + EPS)).
> +	 Compute gamma (X_ADJ + X_EPS) using Stirling's approximation,
> +	 starting by computing pow (X_ADJ, X_ADJ) with a power of 2
> +	 factored out.  */
> +      float exp_adj = -eps;
> +      float x_adj_int = __roundf (x_adj);
> +      float x_adj_frac = x_adj - x_adj_int;
> +      int x_adj_log2;
> +      float x_adj_mant = __frexpf (x_adj, &x_adj_log2);
> +      if (x_adj_mant < (float) M_SQRT1_2)
> +	{
> +	  x_adj_log2--;
> +	  x_adj_mant *= 2.0f;
> +	}
> +      *exp2_adj = x_adj_log2 * (int) x_adj_int;
> +      float ret = (__ieee754_powf (x_adj_mant, x_adj)
> +		   * __ieee754_exp2f (x_adj_log2 * x_adj_frac)
> +		   * __ieee754_expf (-x_adj)
> +		   * __ieee754_sqrtf (2 * (float) M_PI / x_adj)
> +		   / prod);
> +      exp_adj += x_eps * __ieee754_logf (x);
> +      float bsum = gamma_coeff[NCOEFF - 1];
> +      float x_adj2 = x_adj * x_adj;
> +      for (size_t i = 1; i <= NCOEFF - 1; i++)
> +	bsum = bsum / x_adj2 + gamma_coeff[NCOEFF - 1 - i];
> +      exp_adj += bsum / x_adj;
> +      return ret + ret * __expm1f (exp_adj);
> +    }
> +}
>  
>  float
>  __ieee754_gammaf_r (float x, int *signgamp)
>  {
> -  /* We don't have a real gamma implementation now.  We'll use lgamma
> -     and the exp function.  But due to the required boundary
> -     conditions we must check some values separately.  */
>    int32_t hx;
>  
>    GET_FLOAT_WORD (hx, x);
> @@ -50,8 +133,49 @@ __ieee754_gammaf_r (float x, int *signgamp)
>        *signgamp = 0;
>        return x - x;
>      }
> +  if (__builtin_expect ((hx & 0x7f800000) == 0x7f800000, 0))
> +    {
> +      /* Positive infinity (return positive infinity) or NaN (return
> +	 NaN).  */
> +      *signgamp = 0;
> +      return x + x;
> +    }
>  
> -  /* XXX FIXME.  */
> -  return __ieee754_expf (__ieee754_lgammaf_r (x, signgamp));
> +  if (x >= 36.0f)
> +    {
> +      /* Overflow.  */
> +      *signgamp = 0;
> +      return FLT_MAX * FLT_MAX;
> +    }
> +  else if (x > 0.0f)
> +    {
> +      *signgamp = 0;
> +      int exp2_adj;
> +      float ret = gammaf_positive (x, &exp2_adj);
> +      return __scalbnf (ret, exp2_adj);
> +    }
> +  else if (x >= -FLT_EPSILON / 4.0f)
> +    {
> +      *signgamp = 0;
> +      return 1.0f / x;
> +    }
> +  else
> +    {
> +      float tx = __truncf (x);
> +      *signgamp = (tx == 2.0f * __truncf (tx / 2.0f)) ? -1 : 1;
> +      if (x <= -42.0f)
> +	/* Underflow.  */
> +	return FLT_MIN * FLT_MIN;
> +      float frac = tx - x;
> +      if (frac > 0.5f)
> +	frac = 1.0f - frac;
> +      float sinpix = (frac <= 0.25f
> +		      ? __sinf ((float) M_PI * frac)
> +		      : __cosf ((float) M_PI * (0.5f - frac)));
> +      int exp2_adj;
> +      float ret = (float) M_PI / (-x * sinpix
> +				  * gammaf_positive (-x, &exp2_adj));
> +      return __scalbnf (ret, -exp2_adj);
> +    }
>  }
>  strong_alias (__ieee754_gammaf_r, __gammaf_r_finite)
> diff --git a/sysdeps/ieee754/k_standard.c b/sysdeps/ieee754/k_standard.c
> index cd31230..150921f 100644
> --- a/sysdeps/ieee754/k_standard.c
> +++ b/sysdeps/ieee754/k_standard.c
> @@ -837,7 +837,7 @@ __kernel_standard(double x, double y, int type)
>  		exc.type = OVERFLOW;
>  		exc.name = type < 100 ? "tgamma" : (type < 200
>  						   ? "tgammaf" : "tgammal");
> -		exc.retval = HUGE_VAL;
> +		exc.retval = __copysign (HUGE_VAL, x);
>  		if (_LIB_VERSION == _POSIX_)
>  		  __set_errno (ERANGE);
>  		else if (!matherr(&exc)) {
> diff --git a/sysdeps/ieee754/ldbl-128/e_gammal_r.c b/sysdeps/ieee754/ldbl-128/e_gammal_r.c
> index b6da31c..e8d49e9 100644
> --- a/sysdeps/ieee754/ldbl-128/e_gammal_r.c
> +++ b/sysdeps/ieee754/ldbl-128/e_gammal_r.c
> @@ -20,14 +20,108 @@
>  
>  #include <math.h>
>  #include <math_private.h>
> +#include <float.h>
>  
> +/* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's
> +   approximation to gamma function.  */
> +
> +static const long double gamma_coeff[] =
> +  {
> +    0x1.5555555555555555555555555555p-4L,
> +    -0xb.60b60b60b60b60b60b60b60b60b8p-12L,
> +    0x3.4034034034034034034034034034p-12L,
> +    -0x2.7027027027027027027027027028p-12L,
> +    0x3.72a3c5631fe46ae1d4e700dca8f2p-12L,
> +    -0x7.daac36664f1f207daac36664f1f4p-12L,
> +    0x1.a41a41a41a41a41a41a41a41a41ap-8L,
> +    -0x7.90a1b2c3d4e5f708192a3b4c5d7p-8L,
> +    0x2.dfd2c703c0cfff430edfd2c703cp-4L,
> +    -0x1.6476701181f39edbdb9ce625987dp+0L,
> +    0xd.672219167002d3a7a9c886459cp+0L,
> +    -0x9.cd9292e6660d55b3f712eb9e07c8p+4L,
> +    0x8.911a740da740da740da740da741p+8L,
> +    -0x8.d0cc570e255bf59ff6eec24b49p+12L,
> +  };


If all three of these tables are of type long double, why can't
we define the table once and use it everywhere?


> +
> +#define NCOEFF (sizeof (gamma_coeff) / sizeof (gamma_coeff[0]))
> +
> +/* Return gamma (X), for positive X less than 1775, in the form R *
> +   2^(*EXP2_ADJ), where R is the return value and *EXP2_ADJ is set to
> +   avoid overflow or underflow in intermediate calculations.  */
> +
> +static long double
> +gammal_positive (long double x, int *exp2_adj)
> +{
> +  int local_signgam;
> +  if (x < 0.5L)
> +    {
> +      *exp2_adj = 0;
> +      return __ieee754_expl (__ieee754_lgammal_r (x + 1, &local_signgam)) / x;
> +    }
> +  else if (x <= 1.5L)
> +    {
> +      *exp2_adj = 0;
> +      return __ieee754_expl (__ieee754_lgammal_r (x, &local_signgam));
> +    }
> +  else if (x < 12.5L)
> +    {
> +      /* Adjust into the range for using exp (lgamma).  */
> +      *exp2_adj = 0;
> +      long double n = __ceill (x - 1.5L);
> +      long double x_adj = x - n;
> +      long double eps;
> +      long double prod = __gamma_productl (x_adj, 0, n, &eps);
> +      return (__ieee754_expl (__ieee754_lgammal_r (x_adj, &local_signgam))
> +	      * prod * (1.0L + eps));
> +    }
> +  else
> +    {
> +      long double eps = 0;
> +      long double x_eps = 0;
> +      long double x_adj = x;
> +      long double prod = 1;
> +      if (x < 24.0L)
> +	{
> +	  /* Adjust into the range for applying Stirling's
> +	     approximation.  */
> +	  long double n = __ceill (24.0L - x);
> +	  x_adj = x + n;
> +	  x_eps = (x - (x_adj - n));
> +	  prod = __gamma_productl (x_adj - n, x_eps, n, &eps);
> +	}
> +      /* The result is now gamma (X_ADJ + X_EPS) / (PROD * (1 + EPS)).
> +	 Compute gamma (X_ADJ + X_EPS) using Stirling's approximation,
> +	 starting by computing pow (X_ADJ, X_ADJ) with a power of 2
> +	 factored out.  */
> +      long double exp_adj = -eps;
> +      long double x_adj_int = __roundl (x_adj);
> +      long double x_adj_frac = x_adj - x_adj_int;
> +      int x_adj_log2;
> +      long double x_adj_mant = __frexpl (x_adj, &x_adj_log2);
> +      if (x_adj_mant < M_SQRT1_2l)
> +	{
> +	  x_adj_log2--;
> +	  x_adj_mant *= 2.0L;
> +	}
> +      *exp2_adj = x_adj_log2 * (int) x_adj_int;
> +      long double ret = (__ieee754_powl (x_adj_mant, x_adj)
> +			 * __ieee754_exp2l (x_adj_log2 * x_adj_frac)
> +			 * __ieee754_expl (-x_adj)
> +			 * __ieee754_sqrtl (2 * M_PIl / x_adj)
> +			 / prod);
> +      exp_adj += x_eps * __ieee754_logl (x);
> +      long double bsum = gamma_coeff[NCOEFF - 1];
> +      long double x_adj2 = x_adj * x_adj;
> +      for (size_t i = 1; i <= NCOEFF - 1; i++)
> +	bsum = bsum / x_adj2 + gamma_coeff[NCOEFF - 1 - i];
> +      exp_adj += bsum / x_adj;
> +      return ret + ret * __expm1l (exp_adj);
> +    }
> +}
>  
>  long double
>  __ieee754_gammal_r (long double x, int *signgamp)
>  {
> -  /* We don't have a real gamma implementation now.  We'll use lgamma
> -     and the exp function.  But due to the required boundary
> -     conditions we must check some values separately.  */
>    int64_t hx;
>    u_int64_t lx;
>  
> @@ -51,8 +145,49 @@ __ieee754_gammal_r (long double x, int *signgamp)
>        *signgamp = 0;
>        return x - x;
>      }
> +  if ((hx & 0x7fff000000000000ULL) == 0x7fff000000000000ULL)
> +    {
> +      /* Positive infinity (return positive infinity) or NaN (return
> +	 NaN).  */
> +      *signgamp = 0;
> +      return x + x;
> +    }
>  
> -  /* XXX FIXME.  */
> -  return __ieee754_expl (__ieee754_lgammal_r (x, signgamp));
> +  if (x >= 1756.0L)
> +    {
> +      /* Overflow.  */
> +      *signgamp = 0;
> +      return LDBL_MAX * LDBL_MAX;
> +    }
> +  else if (x > 0.0L)
> +    {
> +      *signgamp = 0;
> +      int exp2_adj;
> +      long double ret = gammal_positive (x, &exp2_adj);
> +      return __scalbnl (ret, exp2_adj);
> +    }
> +  else if (x >= -LDBL_EPSILON / 4.0L)
> +    {
> +      *signgamp = 0;
> +      return 1.0f / x;
> +    }
> +  else
> +    {
> +      long double tx = __truncl (x);
> +      *signgamp = (tx == 2.0L * __truncl (tx / 2.0L)) ? -1 : 1;
> +      if (x <= -1775.0L)
> +	/* Underflow.  */
> +	return LDBL_MIN * LDBL_MIN;
> +      long double frac = tx - x;
> +      if (frac > 0.5L)
> +	frac = 1.0L - frac;
> +      long double sinpix = (frac <= 0.25L
> +			    ? __sinl (M_PIl * frac)
> +			    : __cosl (M_PIl * (0.5L - frac)));
> +      int exp2_adj;
> +      long double ret = M_PIl / (-x * sinpix
> +				 * gammal_positive (-x, &exp2_adj));
> +      return __scalbnl (ret, -exp2_adj);
> +    }
>  }
>  strong_alias (__ieee754_gammal_r, __gammal_r_finite)
> diff --git a/sysdeps/ieee754/ldbl-128/gamma_productl.c b/sysdeps/ieee754/ldbl-128/gamma_productl.c
> new file mode 100644
> index 0000000..157dbab
> --- /dev/null
> +++ b/sysdeps/ieee754/ldbl-128/gamma_productl.c
> @@ -0,0 +1,75 @@
> +/* Compute a product of X, X+1, ..., with an error estimate.
> +   Copyright (C) 2013 Free Software Foundation, Inc.
> +   This file is part of the GNU C Library.
> +
> +   The GNU C Library is free software; you can redistribute it and/or
> +   modify it under the terms of the GNU Lesser General Public
> +   License as published by the Free Software Foundation; either
> +   version 2.1 of the License, or (at your option) any later version.
> +
> +   The GNU C Library is distributed in the hope that it will be useful,
> +   but WITHOUT ANY WARRANTY; without even the implied warranty of
> +   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
> +   Lesser General Public License for more details.
> +
> +   You should have received a copy of the GNU Lesser General Public
> +   License along with the GNU C Library; if not, see
> +   <http://www.gnu.org/licenses/>.  */
> +
> +#include <math.h>
> +#include <math_private.h>
> +#include <float.h>
> +
> +/* Calculate X * Y exactly and store the result in *HI + *LO.  It is
> +   given that the values are small enough that no overflow occurs and
> +   large enough (or zero) that no underflow occurs.  */
> +
> +static inline void
> +mul_split (long double *hi, long double *lo, long double x, long double y)
> +{
> +#ifdef __FP_FAST_FMAL
> +  /* Fast built-in fused multiply-add.  */
> +  *hi = x * y;
> +  *lo = __builtin_fmal (x, y, -*hi);
> +#elif defined FP_FAST_FMAL
> +  /* Fast library fused multiply-add, compiler before GCC 4.6.  */
> +  *hi = x * y;
> +  *lo = __fmal (x, y, -*hi);
> +#else
> +  /* Apply Dekker's algorithm.  */
> +  *hi = x * y;
> +# define C ((1LL << (LDBL_MANT_DIG + 1) / 2) + 1)
> +  long double x1 = x * C;
> +  long double y1 = y * C;
> +# undef C
> +  x1 = (x - x1) + x1;
> +  y1 = (y - y1) + y1;
> +  long double x2 = x - x1;
> +  long double y2 = y - y1;
> +  *lo = (((x1 * y1 - *hi) + x1 * y2) + x2 * y1) + x2 * y2;
> +#endif
> +}
> +
> +/* Compute the product of X + X_EPS, X + X_EPS + 1, ..., X + X_EPS + N
> +   - 1, in the form R * (1 + *EPS) where the return value R is an
> +   approximation to the product and *EPS is set to indicate the
> +   approximate error in the return value.  X is such that all the
> +   values X + 1, ..., X + N - 1 are exactly representable, and X_EPS /
> +   X is small enough that factors quadratic in it can be
> +   neglected.  */
> +
> +long double
> +__gamma_productl (long double x, long double x_eps, int n, long double *eps)
> +{
> +  SET_RESTORE_ROUNDL (FE_TONEAREST);
> +  long double ret = x;
> +  *eps = x_eps / x;
> +  for (int i = 1; i < n; i++)
> +    {
> +      *eps += x_eps / (x + i);
> +      long double lo;
> +      mul_split (&ret, &lo, ret, x + i);
> +      *eps += lo / ret;
> +    }
> +  return ret;
> +}
> diff --git a/sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c b/sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c
> index 52ade9e..90d8e3f 100644
> --- a/sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c
> +++ b/sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c
> @@ -20,14 +20,107 @@
>  
>  #include <math.h>
>  #include <math_private.h>
> +#include <float.h>
>  
> +/* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's
> +   approximation to gamma function.  */
> +
> +static const long double gamma_coeff[] =
> +  {
> +    0x1.555555555555555555555555558p-4L,
> +    -0xb.60b60b60b60b60b60b60b60b6p-12L,
> +    0x3.4034034034034034034034034p-12L,
> +    -0x2.7027027027027027027027027p-12L,
> +    0x3.72a3c5631fe46ae1d4e700dca9p-12L,
> +    -0x7.daac36664f1f207daac36664f2p-12L,
> +    0x1.a41a41a41a41a41a41a41a41a4p-8L,
> +    -0x7.90a1b2c3d4e5f708192a3b4c5ep-8L,
> +    0x2.dfd2c703c0cfff430edfd2c704p-4L,
> +    -0x1.6476701181f39edbdb9ce625988p+0L,
> +    0xd.672219167002d3a7a9c886459cp+0L,
> +    -0x9.cd9292e6660d55b3f712eb9e08p+4L,
> +    0x8.911a740da740da740da740da74p+8L,
> +  };
> +
> +#define NCOEFF (sizeof (gamma_coeff) / sizeof (gamma_coeff[0]))
> +
> +/* Return gamma (X), for positive X less than 191, in the form R *
> +   2^(*EXP2_ADJ), where R is the return value and *EXP2_ADJ is set to
> +   avoid overflow or underflow in intermediate calculations.  */
> +
> +static long double
> +gammal_positive (long double x, int *exp2_adj)
> +{
> +  int local_signgam;
> +  if (x < 0.5L)
> +    {
> +      *exp2_adj = 0;
> +      return __ieee754_expl (__ieee754_lgammal_r (x + 1, &local_signgam)) / x;
> +    }
> +  else if (x <= 1.5L)
> +    {
> +      *exp2_adj = 0;
> +      return __ieee754_expl (__ieee754_lgammal_r (x, &local_signgam));
> +    }
> +  else if (x < 11.5L)
> +    {
> +      /* Adjust into the range for using exp (lgamma).  */
> +      *exp2_adj = 0;
> +      long double n = __ceill (x - 1.5L);
> +      long double x_adj = x - n;
> +      long double eps;
> +      long double prod = __gamma_productl (x_adj, 0, n, &eps);
> +      return (__ieee754_expl (__ieee754_lgammal_r (x_adj, &local_signgam))
> +	      * prod * (1.0L + eps));
> +    }
> +  else
> +    {
> +      long double eps = 0;
> +      long double x_eps = 0;
> +      long double x_adj = x;
> +      long double prod = 1;
> +      if (x < 23.0L)
> +	{
> +	  /* Adjust into the range for applying Stirling's
> +	     approximation.  */
> +	  long double n = __ceill (23.0L - x);
> +	  x_adj = x + n;
> +	  x_eps = (x - (x_adj - n));
> +	  prod = __gamma_productl (x_adj - n, x_eps, n, &eps);
> +	}
> +      /* The result is now gamma (X_ADJ + X_EPS) / (PROD * (1 + EPS)).
> +	 Compute gamma (X_ADJ + X_EPS) using Stirling's approximation,
> +	 starting by computing pow (X_ADJ, X_ADJ) with a power of 2
> +	 factored out.  */
> +      long double exp_adj = -eps;
> +      long double x_adj_int = __roundl (x_adj);
> +      long double x_adj_frac = x_adj - x_adj_int;
> +      int x_adj_log2;
> +      long double x_adj_mant = __frexpl (x_adj, &x_adj_log2);
> +      if (x_adj_mant < M_SQRT1_2l)
> +	{
> +	  x_adj_log2--;
> +	  x_adj_mant *= 2.0L;
> +	}
> +      *exp2_adj = x_adj_log2 * (int) x_adj_int;
> +      long double ret = (__ieee754_powl (x_adj_mant, x_adj)
> +			 * __ieee754_exp2l (x_adj_log2 * x_adj_frac)
> +			 * __ieee754_expl (-x_adj)
> +			 * __ieee754_sqrtl (2 * M_PIl / x_adj)
> +			 / prod);
> +      exp_adj += x_eps * __ieee754_logl (x);
> +      long double bsum = gamma_coeff[NCOEFF - 1];
> +      long double x_adj2 = x_adj * x_adj;
> +      for (size_t i = 1; i <= NCOEFF - 1; i++)
> +	bsum = bsum / x_adj2 + gamma_coeff[NCOEFF - 1 - i];
> +      exp_adj += bsum / x_adj;
> +      return ret + ret * __expm1l (exp_adj);
> +    }
> +}
>  
>  long double
>  __ieee754_gammal_r (long double x, int *signgamp)
>  {
> -  /* We don't have a real gamma implementation now.  We'll use lgamma
> -     and the exp function.  But due to the required boundary
> -     conditions we must check some values separately.  */
>    int64_t hx;
>    u_int64_t lx;
>  
> @@ -51,8 +144,49 @@ __ieee754_gammal_r (long double x, int *signgamp)
>        *signgamp = 0;
>        return x - x;
>      }
> +  if ((hx & 0x7ff0000000000000ULL) == 0x7ff0000000000000ULL)
> +    {
> +      /* Positive infinity (return positive infinity) or NaN (return
> +	 NaN).  */
> +      *signgamp = 0;
> +      return x + x;
> +    }
>  
> -  /* XXX FIXME.  */
> -  return __ieee754_expl (__ieee754_lgammal_r (x, signgamp));
> +  if (x >= 172.0L)
> +    {
> +      /* Overflow.  */
> +      *signgamp = 0;
> +      return LDBL_MAX * LDBL_MAX;
> +    }
> +  else if (x > 0.0L)
> +    {
> +      *signgamp = 0;
> +      int exp2_adj;
> +      long double ret = gammal_positive (x, &exp2_adj);
> +      return __scalbnl (ret, exp2_adj);
> +    }
> +  else if (x >= -0x1p-110L)
> +    {
> +      *signgamp = 0;
> +      return 1.0f / x;
> +    }
> +  else
> +    {
> +      long double tx = __truncl (x);
> +      *signgamp = (tx == 2.0L * __truncl (tx / 2.0L)) ? -1 : 1;
> +      if (x <= -191.0L)
> +	/* Underflow.  */
> +	return LDBL_MIN * LDBL_MIN;
> +      long double frac = tx - x;
> +      if (frac > 0.5L)
> +	frac = 1.0L - frac;
> +      long double sinpix = (frac <= 0.25L
> +			    ? __sinl (M_PIl * frac)
> +			    : __cosl (M_PIl * (0.5L - frac)));
> +      int exp2_adj;
> +      long double ret = M_PIl / (-x * sinpix
> +				 * gammal_positive (-x, &exp2_adj));
> +      return __scalbnl (ret, -exp2_adj);
> +    }
>  }
>  strong_alias (__ieee754_gammal_r, __gammal_r_finite)
> diff --git a/sysdeps/ieee754/ldbl-128ibm/gamma_productl.c b/sysdeps/ieee754/ldbl-128ibm/gamma_productl.c
> new file mode 100644
> index 0000000..7c6186d
> --- /dev/null
> +++ b/sysdeps/ieee754/ldbl-128ibm/gamma_productl.c
> @@ -0,0 +1,42 @@
> +/* Compute a product of X, X+1, ..., with an error estimate.
> +   Copyright (C) 2013 Free Software Foundation, Inc.
> +   This file is part of the GNU C Library.
> +
> +   The GNU C Library is free software; you can redistribute it and/or
> +   modify it under the terms of the GNU Lesser General Public
> +   License as published by the Free Software Foundation; either
> +   version 2.1 of the License, or (at your option) any later version.
> +
> +   The GNU C Library is distributed in the hope that it will be useful,
> +   but WITHOUT ANY WARRANTY; without even the implied warranty of
> +   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
> +   Lesser General Public License for more details.
> +
> +   You should have received a copy of the GNU Lesser General Public
> +   License along with the GNU C Library; if not, see
> +   <http://www.gnu.org/licenses/>.  */
> +
> +#include <math.h>
> +#include <math_private.h>
> +
> +/* Compute the product of X + X_EPS, X + X_EPS + 1, ..., X + X_EPS + N
> +   - 1, in the form R * (1 + *EPS) where the return value R is an
> +   approximation to the product and *EPS is set to indicate the
> +   approximate error in the return value.  X is such that all the
> +   values X + 1, ..., X + N - 1 are exactly representable, and X_EPS /
> +   X is small enough that factors quadratic in it can be
> +   neglected.  */
> +
> +long double
> +__gamma_productl (long double x, long double x_eps, int n, long double *eps)
> +{
> +  long double ret = x;
> +  *eps = x_eps / x;
> +  for (int i = 1; i < n; i++)
> +    {
> +      *eps += x_eps / (x + i);
> +      ret *= x + i;
> +      /* FIXME: no error estimates for the multiplication.  */
> +    }
> +  return ret;
> +}
> diff --git a/sysdeps/ieee754/ldbl-96/e_gammal_r.c b/sysdeps/ieee754/ldbl-96/e_gammal_r.c
> index 0974351..7cb3e85 100644
> --- a/sysdeps/ieee754/ldbl-96/e_gammal_r.c
> +++ b/sysdeps/ieee754/ldbl-96/e_gammal_r.c
> @@ -19,14 +19,102 @@
>  
>  #include <math.h>
>  #include <math_private.h>
> +#include <float.h>
>  
> +/* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's
> +   approximation to gamma function.  */
> +
> +static const long double gamma_coeff[] =
> +  {
> +    0x1.5555555555555556p-4L,
> +    -0xb.60b60b60b60b60bp-12L,
> +    0x3.4034034034034034p-12L,
> +    -0x2.7027027027027028p-12L,
> +    0x3.72a3c5631fe46aep-12L,
> +    -0x7.daac36664f1f208p-12L,
> +    0x1.a41a41a41a41a41ap-8L,
> +    -0x7.90a1b2c3d4e5f708p-8L,
> +  };
> +
> +#define NCOEFF (sizeof (gamma_coeff) / sizeof (gamma_coeff[0]))
> +
> +/* Return gamma (X), for positive X less than 1766, in the form R *
> +   2^(*EXP2_ADJ), where R is the return value and *EXP2_ADJ is set to
> +   avoid overflow or underflow in intermediate calculations.  */
> +
> +static long double
> +gammal_positive (long double x, int *exp2_adj)
> +{
> +  int local_signgam;
> +  if (x < 0.5L)
> +    {
> +      *exp2_adj = 0;
> +      return __ieee754_expl (__ieee754_lgammal_r (x + 1, &local_signgam)) / x;
> +    }
> +  else if (x <= 1.5L)
> +    {
> +      *exp2_adj = 0;
> +      return __ieee754_expl (__ieee754_lgammal_r (x, &local_signgam));
> +    }
> +  else if (x < 7.5L)
> +    {
> +      /* Adjust into the range for using exp (lgamma).  */
> +      *exp2_adj = 0;
> +      long double n = __ceill (x - 1.5L);
> +      long double x_adj = x - n;
> +      long double eps;
> +      long double prod = __gamma_productl (x_adj, 0, n, &eps);
> +      return (__ieee754_expl (__ieee754_lgammal_r (x_adj, &local_signgam))
> +	      * prod * (1.0L + eps));
> +    }
> +  else
> +    {
> +      long double eps = 0;
> +      long double x_eps = 0;
> +      long double x_adj = x;
> +      long double prod = 1;
> +      if (x < 13.0L)
> +	{
> +	  /* Adjust into the range for applying Stirling's
> +	     approximation.  */
> +	  long double n = __ceill (13.0L - x);
> +	  x_adj = x + n;
> +	  x_eps = (x - (x_adj - n));
> +	  prod = __gamma_productl (x_adj - n, x_eps, n, &eps);
> +	}
> +      /* The result is now gamma (X_ADJ + X_EPS) / (PROD * (1 + EPS)).
> +	 Compute gamma (X_ADJ + X_EPS) using Stirling's approximation,
> +	 starting by computing pow (X_ADJ, X_ADJ) with a power of 2
> +	 factored out.  */
> +      long double exp_adj = -eps;
> +      long double x_adj_int = __roundl (x_adj);
> +      long double x_adj_frac = x_adj - x_adj_int;
> +      int x_adj_log2;
> +      long double x_adj_mant = __frexpl (x_adj, &x_adj_log2);
> +      if (x_adj_mant < M_SQRT1_2l)
> +	{
> +	  x_adj_log2--;
> +	  x_adj_mant *= 2.0L;
> +	}
> +      *exp2_adj = x_adj_log2 * (int) x_adj_int;
> +      long double ret = (__ieee754_powl (x_adj_mant, x_adj)
> +			 * __ieee754_exp2l (x_adj_log2 * x_adj_frac)
> +			 * __ieee754_expl (-x_adj)
> +			 * __ieee754_sqrtl (2 * M_PIl / x_adj)
> +			 / prod);
> +      exp_adj += x_eps * __ieee754_logl (x);
> +      long double bsum = gamma_coeff[NCOEFF - 1];
> +      long double x_adj2 = x_adj * x_adj;
> +      for (size_t i = 1; i <= NCOEFF - 1; i++)
> +	bsum = bsum / x_adj2 + gamma_coeff[NCOEFF - 1 - i];
> +      exp_adj += bsum / x_adj;
> +      return ret + ret * __expm1l (exp_adj);
> +    }
> +}
>  
>  long double
>  __ieee754_gammal_r (long double x, int *signgamp)
>  {
> -  /* We don't have a real gamma implementation now.  We'll use lgamma
> -     and the exp function.  But due to the required boundary
> -     conditions we must check some values separately.  */
>    u_int32_t es, hx, lx;
>  
>    GET_LDOUBLE_WORDS (es, hx, lx, x);
> @@ -43,22 +131,55 @@ __ieee754_gammal_r (long double x, int *signgamp)
>        *signgamp = 0;
>        return x - x;
>      }
> -  if (__builtin_expect ((es & 0x7fff) == 0x7fff, 0)
> -      && ((hx & 0x7fffffff) | lx) != 0)
> +  if (__builtin_expect ((es & 0x7fff) == 0x7fff, 0))
>      {
> -      /* NaN, return it.  */
> +      /* Positive infinity (return positive infinity) or NaN (return
> +	 NaN).  */
>        *signgamp = 0;
> -      return x;
> +      return x + x;
>      }
> -  if (__builtin_expect ((es & 0x8000) != 0, 0)
> -      && x < 0xffffffff && __rintl (x) == x)
> +  if (__builtin_expect ((es & 0x8000) != 0, 0) && __rintl (x) == x)
>      {
>        /* Return value for integer x < 0 is NaN with invalid exception.  */
>        *signgamp = 0;
>        return (x - x) / (x - x);
>      }
>  
> -  /* XXX FIXME.  */
> -  return __ieee754_expl (__ieee754_lgammal_r (x, signgamp));
> +  if (x >= 1756.0L)
> +    {
> +      /* Overflow.  */
> +      *signgamp = 0;
> +      return LDBL_MAX * LDBL_MAX;
> +    }
> +  else if (x > 0.0L)
> +    {
> +      *signgamp = 0;
> +      int exp2_adj;
> +      long double ret = gammal_positive (x, &exp2_adj);
> +      return __scalbnl (ret, exp2_adj);
> +    }
> +  else if (x >= -LDBL_EPSILON / 4.0L)
> +    {
> +      *signgamp = 0;
> +      return 1.0f / x;
> +    }
> +  else
> +    {
> +      long double tx = __truncl (x);
> +      *signgamp = (tx == 2.0L * __truncl (tx / 2.0L)) ? -1 : 1;
> +      if (x <= -1766.0L)
> +	/* Underflow.  */
> +	return LDBL_MIN * LDBL_MIN;
> +      long double frac = tx - x;
> +      if (frac > 0.5L)
> +	frac = 1.0L - frac;
> +      long double sinpix = (frac <= 0.25L
> +			    ? __sinl (M_PIl * frac)
> +			    : __cosl (M_PIl * (0.5L - frac)));
> +      int exp2_adj;
> +      long double ret = M_PIl / (-x * sinpix
> +				 * gammal_positive (-x, &exp2_adj));
> +      return __scalbnl (ret, -exp2_adj);
> +    }
>  }
>  strong_alias (__ieee754_gammal_r, __gammal_r_finite)
> diff --git a/sysdeps/ieee754/ldbl-96/gamma_product.c b/sysdeps/ieee754/ldbl-96/gamma_product.c
> new file mode 100644
> index 0000000..d464e70
> --- /dev/null
> +++ b/sysdeps/ieee754/ldbl-96/gamma_product.c
> @@ -0,0 +1,46 @@
> +/* Compute a product of X, X+1, ..., with an error estimate.
> +   Copyright (C) 2013 Free Software Foundation, Inc.
> +   This file is part of the GNU C Library.
> +
> +   The GNU C Library is free software; you can redistribute it and/or
> +   modify it under the terms of the GNU Lesser General Public
> +   License as published by the Free Software Foundation; either
> +   version 2.1 of the License, or (at your option) any later version.
> +
> +   The GNU C Library is distributed in the hope that it will be useful,
> +   but WITHOUT ANY WARRANTY; without even the implied warranty of
> +   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
> +   Lesser General Public License for more details.
> +
> +   You should have received a copy of the GNU Lesser General Public
> +   License along with the GNU C Library; if not, see
> +   <http://www.gnu.org/licenses/>.  */
> +
> +#include <math.h>
> +#include <math_private.h>
> +#include <float.h>
> +
> +/* Compute the product of X + X_EPS, X + X_EPS + 1, ..., X + X_EPS + N
> +   - 1, in the form R * (1 + *EPS) where the return value R is an
> +   approximation to the product and *EPS is set to indicate the
> +   approximate error in the return value.  X is such that all the
> +   values X + 1, ..., X + N - 1 are exactly representable, and X_EPS /
> +   X is small enough that factors quadratic in it can be
> +   neglected.  */
> +
> +double
> +__gamma_product (double x, double x_eps, int n, double *eps)
> +{
> +  long double x_full = (long double) x + (long double) x_eps;
> +  long double ret = x_full;
> +  for (int i = 1; i < n; i++)
> +    ret *= x_full + i;
> +
> +#if FLT_EVAL_METHOD != 0
> +  volatile
> +#endif
> +  double fret = ret;
> +  *eps = (ret - fret) / fret;
> +
> +  return fret;
> +}
> diff --git a/sysdeps/ieee754/ldbl-96/gamma_productl.c b/sysdeps/ieee754/ldbl-96/gamma_productl.c
> new file mode 100644
> index 0000000..157dbab
> --- /dev/null
> +++ b/sysdeps/ieee754/ldbl-96/gamma_productl.c
> @@ -0,0 +1,75 @@
> +/* Compute a product of X, X+1, ..., with an error estimate.
> +   Copyright (C) 2013 Free Software Foundation, Inc.
> +   This file is part of the GNU C Library.
> +
> +   The GNU C Library is free software; you can redistribute it and/or
> +   modify it under the terms of the GNU Lesser General Public
> +   License as published by the Free Software Foundation; either
> +   version 2.1 of the License, or (at your option) any later version.
> +
> +   The GNU C Library is distributed in the hope that it will be useful,
> +   but WITHOUT ANY WARRANTY; without even the implied warranty of
> +   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
> +   Lesser General Public License for more details.
> +
> +   You should have received a copy of the GNU Lesser General Public
> +   License along with the GNU C Library; if not, see
> +   <http://www.gnu.org/licenses/>.  */
> +
> +#include <math.h>
> +#include <math_private.h>
> +#include <float.h>
> +
> +/* Calculate X * Y exactly and store the result in *HI + *LO.  It is
> +   given that the values are small enough that no overflow occurs and
> +   large enough (or zero) that no underflow occurs.  */
> +
> +static inline void
> +mul_split (long double *hi, long double *lo, long double x, long double y)
> +{
> +#ifdef __FP_FAST_FMAL
> +  /* Fast built-in fused multiply-add.  */
> +  *hi = x * y;
> +  *lo = __builtin_fmal (x, y, -*hi);
> +#elif defined FP_FAST_FMAL
> +  /* Fast library fused multiply-add, compiler before GCC 4.6.  */
> +  *hi = x * y;
> +  *lo = __fmal (x, y, -*hi);
> +#else
> +  /* Apply Dekker's algorithm.  */
> +  *hi = x * y;
> +# define C ((1LL << (LDBL_MANT_DIG + 1) / 2) + 1)
> +  long double x1 = x * C;
> +  long double y1 = y * C;
> +# undef C
> +  x1 = (x - x1) + x1;
> +  y1 = (y - y1) + y1;
> +  long double x2 = x - x1;
> +  long double y2 = y - y1;
> +  *lo = (((x1 * y1 - *hi) + x1 * y2) + x2 * y1) + x2 * y2;
> +#endif
> +}
> +
> +/* Compute the product of X + X_EPS, X + X_EPS + 1, ..., X + X_EPS + N
> +   - 1, in the form R * (1 + *EPS) where the return value R is an
> +   approximation to the product and *EPS is set to indicate the
> +   approximate error in the return value.  X is such that all the
> +   values X + 1, ..., X + N - 1 are exactly representable, and X_EPS /
> +   X is small enough that factors quadratic in it can be
> +   neglected.  */
> +
> +long double
> +__gamma_productl (long double x, long double x_eps, int n, long double *eps)
> +{
> +  SET_RESTORE_ROUNDL (FE_TONEAREST);
> +  long double ret = x;
> +  *eps = x_eps / x;
> +  for (int i = 1; i < n; i++)
> +    {
> +      *eps += x_eps / (x + i);
> +      long double lo;
> +      mul_split (&ret, &lo, ret, x + i);
> +      *eps += lo / ret;
> +    }
> +  return ret;
> +}
> diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
> index d84a898..3827b9d 100644
> --- a/sysdeps/x86_64/fpu/libm-test-ulps
> +++ b/sysdeps/x86_64/fpu/libm-test-ulps
> @@ -7142,6 +7142,417 @@ idouble: 1
>  ifloat: 1
>  ildouble: 1
>  ldouble: 1
> +Test "tgamma (-0x0.fffffffffffff8p0) == -9.0071992547409924227843350984672961392521e+15":
> +double: 1
> +idouble: 1
> +Test "tgamma (-0x0.ffffffffffffffffp0) == -1.8446744073709551616422784335098467139470e+19":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x0.ffffffp0) == -1.6777216422784419250710305882992376932423e+07":
> +float: 1
> +ifloat: 1
> +Test "tgamma (-0x1.000002p0) == 8.3886075772158332060084424806449513922858e+06":
> +double: 2
> +idouble: 2
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x1.0a32a2p+5) == 1.8125267978155035272941154746083439329912e-37":
> +float: 2
> +ifloat: 2
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x1.5800000080001p+7) == -3.1439271448823567326093363350637118195240e-304":
> +ildouble: 2
> +ldouble: 2
> +Test "tgamma (-0x1.fffffffffffffp0) == 2.2517998136852484613921675492337776673289e+15":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x13.ffffep0) == 2.1550026214525536756224040483579183652119e-13":
> +float: 2
> +ifloat: 2
> +Test "tgamma (-0x13.ffffffffffffffep0) == 2.3694367893405502075347562184931828448654e-01":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x13.ffffffffffffp0) == 1.1569515572952029402736625857313236848570e-04":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x14.000000000001p0) == -1.1569515572951781096476686854873801225397e-04":
> +double: 1
> +idouble: 1
> +ildouble: 2
> +ldouble: 2
> +Test "tgamma (-0x14.00002p0) == -2.1549777908265594916405421768142757507179e-13":
> +float: 1
> +ifloat: 1
> +Test "tgamma (-0x1d.ffffep0) == 1.9765721589464867957912772592816027583176e-27":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x1d.ffffffffffffffep0) == 2.1732499046818166459536268654187775086902e-15":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x1e.000000000000002p0) == -2.1732499046818166201837145753965837196590e-15":
> +ildouble: 2
> +ldouble: 2
> +Test "tgamma (-0x1e.000000000001p0) == -1.0611571800204053929094168642022073530425e-18":
> +double: 3
> +idouble: 3
> +Test "tgamma (-0x1e.00002p0) == -1.9765463890341964384070157599286498212650e-27":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x1f4.00000000000002p0) == -2.9528489142763131406565394149878256133744e-1118":
> +ildouble: 3
> +ldouble: 3
> +Test "tgamma (-0x1p-24) == -1.6777216577215723853867349114260580375249e+07":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x2.0000000000000004p0) == -2.3058430092136939515386078324507664305064e+18":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x2.0000000000002p0) == -1.1258999068426235386078324507668462444260e+15":
> +double: 1
> +idouble: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x2.000004p0) == -2.0971515386080557574407223895988378776747e+06":
> +double: 2
> +float: 1
> +idouble: 2
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x2.fffffcp0) == -6.9905087601970247876992248591045142913324e+05":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +Test "tgamma (-0x27.ffffcp0) == 3.2129279441390812141195076945616975790225e-43":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x27.fffffffffffep0) == 1.7249032006742266376460389310340465554361e-34":
> +double: 1
> +idouble: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x28.000000000002p0) == -1.7249032006741359094184881234822934593822e-34":
> +double: 1
> +idouble: 1
> +Test "tgamma (-0x28.00004p0) == -3.2128372159115252365699015758097981155793e-43":
> +double: 2
> +idouble: 2
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x28.ffffcp0) == -7.8364103489619817539676737414096652170685e-45":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x28.fffffffffffep0) == -4.2070809772542120404320040128839297118648e-36":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x28.ffffffffffffffcp0) == -8.6161018414163982777002940498289948893044e-33":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x29.000000000000004p0) == 8.6161018414163980549537337663264762179535e-33":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x29.000000000002p0) == 4.2070809772539892938717205103652583609422e-36":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x29.00004p0) == 7.8361876024016854597745353972619195760515e-45":
> +double: 1
> +idouble: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x29.ffffcp0) == 1.8658121573125798145204120066590953505132e-46":
> +double: 1
> +idouble: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x2a.00004p0) == -1.8657587834931410688246126853566488626385e-46":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x2ed.fffffffffffffcp0) == 6.9801511765871818502006905472380418430269e-1817":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x3.0000000000002p0) == 3.7529996894754112398038859470009084971438e+14":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x3.000004p0) == 6.9905045731381300146131914617735687322025e+05":
> +double: 2
> +float: 1
> +idouble: 2
> +ifloat: 1
> +Test "tgamma (-0x3.fffffcp0) == 1.7476272942159602684441970627092458855771e+05":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +Test "tgamma (-0x3.ffffffffffffep0) == 9.3824992236885396088236184658402406857503e+13":
> +double: 2
> +idouble: 2
> +Test "tgamma (-0x3.fffffffffffffffcp0) == 1.9215358410114116272942156951799168638773e+17":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x31.fffffffffffep0) == 4.6273774273632946947805289899230181990085e-51":
> +double: 3
> +idouble: 3
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x32.000000000000004p0) == -9.4768689712397633101385547903658075308777e-48":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x32.000000000002p0) == -4.6273774273630367887073532197576655720178e-51":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x3e7.fffffffffffffcp0) == 4.4768809295877296071892611539415773519036e-2552":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x3e8.00000000000004p0) == -4.4768809295877261735541135972060089530309e-2552":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x4.0000000000000008p0) == -9.6076792050570581270578430482008313684602e+16":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x4.000008p0) == -8.7381270578483499672965708923121931082305e+04":
> +float: 1
> +ifloat: 1
> +Test "tgamma (-0x4.fffff8p0) == -1.7476280884325863043793087474680780379554e+04":
> +double: 1
> +idouble: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x4.ffffffffffffcp0) == -9.3824992236885475509805702650262155809819e+12":
> +double: 1
> +idouble: 1
> +Test "tgamma (-0x4e2.00000000000008p0) == -5.4651488569236421026544487194247355967789e-3315":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x5.0000000000000008p0) == 1.9215358410114116252449019429734996071487e+16":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x5.0000000000004p0) == 9.3824992236885191156860964016850034672946e+12":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x5.000008p0) == 1.7476252449031389167286893378510439443844e+04":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x5.ffffffffffffcp0) == 1.5637498706147581566449098589862357835505e+12":
> +double: 1
> +idouble: 1
> +Test "tgamma (-0x5db.fffffffffffff8p0) == 1.8718211510339187689122114747834510481993e-4099":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x6.0000000000000008p0) == -3.2025597350190193751766884234743511972877e+15":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x6.000008p0) == -2.9127085100239567622341538102130981196910e+03":
> +float: 2
> +ifloat: 2
> +Test "tgamma (-0x6.fffff8p0) == -4.1610198723079349791939054365613377035519e+02":
> +double: 2
> +float: 1
> +idouble: 2
> +ifloat: 1
> +ildouble: 2
> +ldouble: 2
> +Test "tgamma (-0x6.ffffffffffffcp0) == -2.2339283865925119357965832452642909859289e+11":
> +double: 2
> +idouble: 2
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x63.fffffffffffcp0) == 7.5400833348840965463348754984345825364294e-145":
> +double: 2
> +idouble: 2
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x63.ffffffffffffff8p0) == 1.5442090669841618542494279375256856430049e-141":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x64.000000000000008p0) == -1.5442090669841617554527108348771968070612e-141":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x64.000000000004p0) == -7.5400833348831085791638490135462230991587e-145":
> +double: 1
> +idouble: 1
> +ildouble: 2
> +ldouble: 2
> +Test "tgamma (-0x6d6.00000000000008p0) == -4.2925786447266421378134368786479937285900e-4902":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x7.0000000000000008p0) == 4.5750853357414562499689653215166468353753e+14":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x7.0000000000004p0) == 2.2339283865925039372192897706214475877342e+11":
> +double: 3
> +idouble: 3
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x7.000008p0) == 4.1610118737306415004517215226199741948733e+02":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x7.fffff8p0) == 5.2012751504050764429534086402871289946986e+01":
> +double: 3
> +float: 1
> +idouble: 3
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x7.ffffffffffffcp0) == 2.7924104832406402297655703264222230055898e+10":
> +double: 3
> +idouble: 3
> +Test "tgamma (-0x7.fffffffffffffff8p0) == 5.7188566696768203227694481100089533685959e+13":
> +ildouble: 4
> +ldouble: 4
> +Test "tgamma (-0x8.000000000000001p0) == -2.8594283348384101534210280804672371201060e+13":
> +ildouble: 2
> +ldouble: 2
> +Test "tgamma (-0x8.00001p0) == -2.6006296115134418896533598545925084576702e+01":
> +double: 2
> +idouble: 2
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x9.ffffffffffff8p0) == 1.5513391573559147700413058496716749249803e+08":
> +double: 1
> +idouble: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x9.fffffp0) == 2.8896008370721717567612135720915723136310e-01":
> +float: 1
> +ifloat: 1
> +Test "tgamma (-0x95.fffffffffff8p0) == 6.1582369322723207086020016423767264008839e-250":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x95.ffffffffffffffp0) == 1.2612069237291916400144732227892704713839e-246":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0x96.000000000008p0) == -6.1582369322705655439003240743176243138734e-250":
> +double: 1
> +idouble: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0xa.000000000000001p0) == -3.1771425942649001698860433502350057763905e+11":
> +ildouble: 2
> +ldouble: 2
> +Test "tgamma (-0xa.00001p0) == -2.8895878754728051776830454190076999107021e-01":
> +double: 1
> +idouble: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0xb4.ffffffffffffffp0) == -1.9816628031468191243385005680879281767694e-315":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0xb5.00000000000001p0) == 1.9816628031468188382579700510291588022368e-315":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0xb5.000000000008p0) == 9.6760879059888966544677044221698800670218e-319":
> +ildouble: 2
> +ldouble: 2
> +Test "tgamma (-0xb5.ffffffffffffffp0) == 1.0888257160147357826865964233809723297472e-317":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0xb6.00000000000001p0) == -1.0888257160147356253334423783317128355514e-317":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0xb6.fffffffffff8p0) == -2.9052086428846935908287469917922960610289e-323":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0xb6.ffffffffffffffp0) == -5.9498673006269714905418984659220067091260e-320":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0xb7.fffffffffff8p0) == 1.5789177406982032823826953250736039527543e-325":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0xb7.ffffffffffffffp0) == 3.2336235329494410277123118903958061569834e-322":
> +ildouble: 2
> +ldouble: 2
> +Test "tgamma (-0xb8.000000000008p0) == -1.5789177406977349925854817486109369828857e-325":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0xbb.ffffffffffffffp0) == 2.6730392040715350119087465463119939092815e-331":
> +ildouble: 2
> +ldouble: 2
> +Test "tgamma (-0xbc.00000000000001p0) == -2.6730392040715346232108532050343031951651e-331":
> +ildouble: 3
> +ldouble: 3
> +Test "tgamma (-0xbd.00000000000001p0) == 1.4143064571807061497431633629389135273431e-333":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0xbe.00000000000001p0) == -7.4437181956879271033676895858841525581153e-336":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0xbe.ffffffffffffffp0) == -3.8972346574282346536709453101948570578636e-338":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0xbf.00000000000001p0) == 3.8972346574282340852496542564155275274974e-338":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0xf9.ffffffffffffffp0) == 2.2289142548411573883553287678043297937797e-476":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-0xfa.00000000000001p0) == -2.2289142548411570466476165308364665814265e-476":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-1.5) == 2.3632718012073547030642233111215269103967e+00":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-2.5) == -9.4530872048294188122568932444861076415869e-01":
> +double: 1
> +float: 2
> +idouble: 1
> +ifloat: 2
> +Test "tgamma (-3.5) == 2.7008820585226910892162552127103164690248e-01":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +Test "tgamma (-4.5) == -6.0019601300504246427027893615784810422774e-02":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-5.5) == 1.0912654781909862986732344293779056440504e-02":
> +double: 1
> +idouble: 1
> +Test "tgamma (-6.5) == -1.6788699664476712287280529682737009908468e-03":
> +float: 1
> +ifloat: 1
> +Test "tgamma (-7.5) == 2.2384932885968949716374039576982679877958e-04":
> +double: 2
> +float: 1
> +idouble: 2
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-8.5) == -2.6335215159963470254557693619979623385833e-05":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (-9.5) == 2.7721279115751021320587045915768024616666e-06":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +ildouble: 2
> +ldouble: 2
>  Test "tgamma (0.5) == sqrt (pi)":
>  float: 1
>  ifloat: 1
> @@ -7150,7 +7561,260 @@ double: 1
>  float: 1
>  idouble: 1
>  ifloat: 1
> +Test "tgamma (0x1.fffffep0) == 9.9999994960018563231526611134590489120697e-01":
> +float: 1
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x1.fffffffffffffffep0) == 9.9999999999999999995416163053934024243282e-01":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x1.fffffffffffffp0) == 9.9999999999999990612301934456883679778984e-01":
> +double: 1
> +idouble: 1
> +Test "tgamma (0x1p-24) == 1.6777215422784394050795179874582764575261e+07":
> +float: 1
> +ifloat: 1
> +Test "tgamma (0x1p-53) == 9.0071992547409914227843350984672492007618e+15":
> +double: 1
> +idouble: 1
> +Test "tgamma (0x1p-64) == 1.8446744073709551615422784335098467139447e+19":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x2.0000000000002p0) == 1.0000000000000001877539613108624482361963e+00":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x2.000004p0) == 1.0000001007996638509889062631687945799175e+00":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x2.30a43cp+4) == 3.4027979115654976101247558405326779640190e+38":
> +double: 1
> +float: 2
> +idouble: 1
> +ifloat: 2
> +Test "tgamma (0x2.fffffcp0) == 1.9999995599822108706107786027549565954046e+00":
> +float: 3
> +ifloat: 3
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x2.ffffffffffffep0) == 1.9999999999999991804028675282128956223990e+00":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x3.0000000000002p0) == 2.0000000000000008195971324717875960213536e+00":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x3.fffffcp0) == 5.9999982031095793171233994481968816873643e+00":
> +float: 1
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x3.ffffffffffffep0) == 5.9999999999999966530301828845138185025345e+00":
> +double: 1
> +idouble: 1
> +Test "tgamma (0x3.fffffffffffffffcp0) == 5.9999999999999999983657373939865784753909e+00":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x4.0000000000000008p0) == 6.0000000000000000032685252120268430507939e+00":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x4.0000000000004p0) == 6.0000000000000066939396342309789716341613e+00":
> +double: 1
> +idouble: 1
> +Test "tgamma (0x4.fffff8p0) == 2.3999982763857938712639837029547357501709e+01":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x4.ffffffffffffcp0) == 2.3999999999999967895170944875373910918544e+01":
> +double: 1
> +idouble: 1
> +Test "tgamma (0x4.fffffffffffffff8p0) == 2.3999999999999999984323813937927417165027e+01":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x5.0000000000000008p0) == 2.4000000000000000015676186062072582846211e+01":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x5.0000000000004p0) == 2.4000000000000032104829055124673225982803e+01":
> +double: 1
> +idouble: 1
> +Test "tgamma (0x5.000008p0) == 2.4000017236155647574166073485628713443799e+01":
> +float: 2
> +ifloat: 2
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x5.fffff8p0) == 1.1999990237520611552119807476573441975106e+02":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x5.ffffffffffffcp0) == 1.1999999999999981815957265157389249327533e+02":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x6.0000000000000008p0) == 1.2000000000000000008878927116622375680433e+02":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x6.0000000000004p0) == 1.2000000000000018184042734842640022086408e+02":
> +double: 1
> +idouble: 1
> +Test "tgamma (0x6.000008p0) == 1.2000009762487825358530770343720418162783e+02":
> +float: 2
> +ifloat: 2
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x6.fffff8p0) == 7.1999935703082425988147448928288557689866e+02":
> +double: 1
> +idouble: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x6.ffffffffffffcp0) == 7.1999999999999880237602554542848858572672e+02":
> +double: 3
> +idouble: 3
> +Test "tgamma (0x7.0000000000000008p0) == 7.2000000000000000058477733127664675369681e+02":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x7.0000000000004p0) == 7.2000000000000119762397445457359071259652e+02":
> +double: 4
> +idouble: 4
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x7.000008p0) == 7.2000064296977505705636258629805621178587e+02":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x7.fffff8p0) == 5.0399951558933225045148935487583089307135e+03":
> +double: 2
> +float: 1
> +idouble: 2
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x7.ffffffffffffcp0) == 5.0399999999999909771437166339103165198442e+03":
> +double: 2
> +idouble: 2
> +Test "tgamma (0x7.fffffffffffffff8p0) == 5.0399999999999999955943084553876474508520e+03":
> +ildouble: 3
> +ldouble: 3
> +Test "tgamma (0x8.000000000000001p0) == 5.0400000000000000088113830892247051102283e+03":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0x8.00001p0) == 5.0400096882277802019946778420223050233915e+03":
> +double: 2
> +idouble: 2
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (0xa.b9fd72b0fb238p+4) == 1.7976931348622298700886249281842651388250e+308":
> +double: 1
> +idouble: 1
> +ildouble: 2
> +ldouble: 2
> +Test "tgamma (10) == 362880":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +Test "tgamma (18.5) == 1.4986120533153361177371791123515513270334e+15":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (19.5) == 2.7724322986333718178137813578503699550119e+16":
> +double: 2
> +idouble: 2
> +Test "tgamma (2.5) == 1.3293403881791370204736256125058588870982e+00":
> +float: 2
> +ifloat: 2
> +Test "tgamma (23.5) == 5.3613035875444147334274983856108155717836e+21":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (29.5) == 1.6348125198274266444378807806868221866931e+30":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (3) == 2":
> +float: 1
> +ifloat: 1
> +Test "tgamma (3.5) == 3.3233509704478425511840640312646472177454e+00":
> +float: 2
> +ifloat: 2
> +Test "tgamma (30.5) == 4.8226969334909086010917483030261254507447e+31":
> +float: 1
> +ifloat: 1
> +Test "tgamma (32.5) == 4.6334060788513904384988971821323500268029e+34":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (33.5) == 1.5058569756267018925121415841930137587110e+36":
> +float: 1
> +ifloat: 1
> +Test "tgamma (34.5) == 5.0446208683494513399156743070465960916817e+37":
> +double: 1
> +float: 2
> +idouble: 1
> +ifloat: 2
> +ildouble: 1
> +ldouble: 1
>  Test "tgamma (4) == 6":
> +float: 1
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (4.5) == 1.1631728396567448929144224109426265262109e+01":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +Test "tgamma (5.5) == 5.2342777784553520181149008492418193679490e+01":
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (6) == 120":
> +float: 1
> +ifloat: 1
> +Test "tgamma (6.5) == 2.8788527781504436099631954670830006523720e+02":
> +float: 1
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (7) == 720":
> +double: 1
> +idouble: 1
> +Test "tgamma (7.5) == 1.8712543057977883464760770536039504240418e+03":
> +double: 2
> +float: 1
> +idouble: 2
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (8) == 5040":
> +double: 1
> +idouble: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (8.5) == 1.4034407293483412598570577902029628180313e+04":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +ildouble: 2
> +ldouble: 2
> +Test "tgamma (9) == 40320":
> +double: 1
> +idouble: 1
> +ildouble: 1
> +ldouble: 1
> +Test "tgamma (9.5) == 1.1929246199460900708784991216725183953266e+05":
> +double: 1
> +idouble: 1
>  ildouble: 1
>  ldouble: 1
>  
> @@ -8052,12 +8716,12 @@ ildouble: 2
>  ldouble: 2
>  
>  Function: "tgamma":
> -double: 1
> -float: 1
> -idouble: 1
> -ifloat: 1
> -ildouble: 1
> -ldouble: 1
> +double: 4
> +float: 3
> +idouble: 4
> +ifloat: 3
> +ildouble: 4
> +ldouble: 4
>  
>  Function: "y0":
>  double: 2
> 

Everything else looks fine.

Cheers,
Carlos.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]