diff --git a/sysdeps/ieee754/flt-32/e_gammaf_r.c b/sysdeps/ieee754/flt-32/e_gammaf_r.c index 90ed3b4890..6b1f95d50f 100644 --- a/sysdeps/ieee754/flt-32/e_gammaf_r.c +++ b/sysdeps/ieee754/flt-32/e_gammaf_r.c @@ -37,9 +37,7 @@ SOFTWARE. #include #include #include - -typedef union {float f; uint32_t u;} b32u32_u; -typedef union {double f; uint64_t u;} b64u64_u; +#include "math_config.h" float __ieee754_gammaf_r (float x, int *signgamp) @@ -54,97 +52,125 @@ __ieee754_gammaf_r (float x, int *signgamp) /* List of exceptional cases. Each entry contains the 32-bit encoding u of x, a binary32 approximation f of gamma(x), and a correction term df. */ - static const struct {uint32_t u; float f, df;} tb[] = { - {0x27de86a9u, 0x1.268266p+47f, 0x1p22f}, // x = 0x1.bd0d52p-48 - {0x27e05475u, 0x1.242422p+47f, 0x1p22f}, // x = 0x1.c0a8eap-48 - {0xb63befb3u, -0x1.5cb6e4p+18f, 0x1p-7f}, // x = -0x1.77df66p-19 - {0x3c7bb570u, 0x1.021d9p+6f, 0x1p-19f}, // x = 0x1.f76aep-7 - {0x41e886d1u, 0x1.33136ap+98f, 0x1p73f}, // x = 0x1.d10da2p+4 - {0xc067d177u, 0x1.f6850cp-3f, 0x1p-28f}, // x = -0x1.cfa2eep+1 - {0xbd99da31u, -0x1.befe66p+3, -0x1p-22f}, // x = -0x1.33b462p-4 - {0xbf54c45au, -0x1.a6b4ecp+2, +0x1p-23f}, // x = -0x1.a988b4p-1 - {0x41ee77feu, 0x1.d3631cp+101, -0x1p-76f}, // x = 0x1.dceffcp+4 - {0x3f843a64u, 0x1.f6c638p-1, 0x1p-26f}, // x = 0x1.0874c8p+0 + static const struct + { + uint32_t u; + float f, df; + } tb[] = { + { 0x27de86a9u, 0x1.268266p+47f, 0x1p22f }, /* x = 0x1.bd0d52p-48 */ + { 0x27e05475u, 0x1.242422p+47f, 0x1p22f }, /* x = 0x1.c0a8eap-48 */ + { 0xb63befb3u, -0x1.5cb6e4p+18f, 0x1p-7f }, /* x = -0x1.77df66p-19 */ + { 0x3c7bb570u, 0x1.021d9p+6f, 0x1p-19f }, /* x = 0x1.f76aep-7 */ + { 0x41e886d1u, 0x1.33136ap+98f, 0x1p73f }, /* x = 0x1.d10da2p+4 */ + { 0xc067d177u, 0x1.f6850cp-3f, 0x1p-28f }, /* x = -0x1.cfa2eep+1 */ + { 0xbd99da31u, -0x1.befe66p+3, -0x1p-22f }, /* x = -0x1.33b462p-4 */ + { 0xbf54c45au, -0x1.a6b4ecp+2, 0x1p-23f }, /* x = -0x1.a988b4p-1 */ + { 0x41ee77feu, 0x1.d3631cp+101, -0x1p-76f }, /* x = 0x1.dceffcp+4 */ + { 0x3f843a64u, 0x1.f6c638p-1, 0x1p-26f }, /* x = 0x1.0874c8p+0 */ }; - b32u32_u t = {.f = x}; - uint32_t ax = t.u<<1; - if(__builtin_expect(ax>=(0xffu<<24), 0)){ /* x=NaN or +/-Inf */ - if(ax==(0xffu<<24)){ /* x=+/-Inf */ - if(t.u>>31){ /* x=-Inf */ - return x / x; /* will raise the "Invalid operation" exception */ - } - return x; /* x=+Inf */ + uint32_t t = asuint (x); + uint32_t ax = t << 1; + if (__glibc_unlikely (ax >= (0xffu << 24))) + { /* x=NaN or +/-Inf */ + if (ax == (0xffu << 24)) + { /* x=+/-Inf */ + if (t >> 31) /* x=-Inf */ + return __math_invalidf (x); + return x; /* x=+Inf */ + } + return x + x; /* x=NaN, where x+x ensures the "Invalid operation" + exception is set if x is sNaN */ } - return x + x; /* x=NaN, where x+x ensures the "Invalid operation" - exception is set if x is sNaN */ - } double z = x; - if(__builtin_expect(ax<0x6d000000u, 0)){ /* |x| < 0x1p-18 */ - volatile double d = (0x1.fa658c23b1578p-1 - 0x1.d0a118f324b63p-1*z)*z - 0x1.2788cfc6fb619p-1; - double f = 1.0/z + d; - float r = f; - b64u64_u rt = {.f = f}; - if(((rt.u+2)&0xfffffff) < 4){ - for(unsigned i=0;i= 0x1.18522p+5f)) + { + /* Overflow case. The original CORE-MATH code returns + 0x1p127f * 0x1p127f, but apparently some compilers replace this + by +Inf. */ + return math_narrow_eval (x * 0x1p127f); } - return r; - } - float fx = __builtin_floorf(x); - if(__builtin_expect(x >= 0x1.18522p+5f, 0)){ - /* Overflow case. The original CORE-MATH code returns 0x1p127f * 0x1p127f, - but apparently some compilers replace this by +Inf. */ - return math_narrow_eval (x * 0x1p127f); - } /* compute k only after the overflow check, otherwise the case to integer might overflow */ int k = fx; - if(__builtin_expect(fx==x, 0)){ /* x is integer */ - if(x == 0.0f){ - return 1.0f/x; + if (__glibc_unlikely (fx == x)) + { /* x is integer */ + if (x == 0.0f) + return 1.0f / x; + if (x < 0.0f) + return __math_invalidf (0.0f); + double t0 = 1, x0 = 1; + for (int i = 1; i < k; i++, x0 += 1.0) + t0 *= x0; + return t0; } - if(x < 0.0f){ - return 0.0f / 0.0f; /* should raise the "Invalid operation" exception */ + if (__glibc_unlikely (x < -42.0f)) + { /* negative non-integer */ + /* For x < -42, x non-integer, |gamma(x)| < 2^-151. */ + static const float sgn[2] = { 0x1p-127f, -0x1p-127f }; + /* Underflows always happens */ + return math_narrow_eval (0x1p-127f * sgn[k & 1]); } - double t0 = 1, x0 = 1; - for(int i=1; i