From 345e9c7d0b36922e790e43bc4a75c40664e7981a Mon Sep 17 00:00:00 2001 From: Adhemerval Zanella Date: Fri, 25 Oct 2024 15:21:39 -0300 Subject: [PATCH] math: Add e_gammaf_r to glibc code and style Also remove the use of builtins in favor of standard names, compiler already inline them (if supported) with current compiler options. It also fixes and issue where __builtin_roundeven is not support on gcc older than version 10. Checked on x86_64-linux-gnu and i686-linux_gnu. Signed-off-by: Adhemerval Zanella Reviewed-by: DJ Delorie --- sysdeps/ieee754/flt-32/e_gammaf_r.c | 180 ++++++++++++++++------------ sysdeps/m68k/m680x0/fpu/math_errf.c | 1 - 2 files changed, 103 insertions(+), 78 deletions(-) delete mode 100644 sysdeps/m68k/m680x0/fpu/math_errf.c 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