@@ -37,9 +37,7 @@ SOFTWARE.
#include <stddef.h>
#include <libm-alias-finite.h>
#include <math-narrow-eval.h>
-
-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,123 @@ __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<sizeof(tb)/sizeof(tb[0]);i++)
- if(t.u==tb[i].u) return tb[i].f + tb[i].df;
+ if (__glibc_unlikely (ax < 0x6d000000u))
+ { /* |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;
+ uint64_t rt = asuint64 (f);
+ if (((rt + 2) & 0xfffffff) < 4)
+ {
+ for (unsigned i = 0; i < sizeof (tb) / sizeof (tb[0]); i++)
+ if (t == tb[i].u)
+ return tb[i].f + tb[i].df;
+ }
+ return r;
+ }
+ float fx = floorf (x);
+ if (__glibc_unlikely (x >= 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<k; i++, x0 += 1.0) t0 *= x0;
- return t0;
- }
- if(__builtin_expect(x<-42.0f, 0)){ /* 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]);
- }
- /* The array c[] stores a degree-15 polynomial approximation for gamma(x). */
+ /* The array c[] stores a degree-15 polynomial approximation for gamma(x). */
static const double c[] =
- {0x1.c9a76be577123p+0, 0x1.8f2754ddcf90dp+0, 0x1.0d1191949419bp+0, 0x1.e1f42cf0ae4a1p-2,
- 0x1.82b358a3ab638p-3, 0x1.e1f2b30cd907bp-5, 0x1.240f6d4071bd8p-6, 0x1.1522c9f3cd012p-8,
- 0x1.1fd0051a0525bp-10, 0x1.9808a8b96c37ep-13, 0x1.b3f78e01152b5p-15, 0x1.49c85a7e1fd04p-18,
- 0x1.471ca49184475p-19, -0x1.368f0b7ed9e36p-23, 0x1.882222f9049efp-23, -0x1.a69ed2042842cp-25};
+ {
+ 0x1.c9a76be577123p+0, 0x1.8f2754ddcf90dp+0, 0x1.0d1191949419bp+0,
+ 0x1.e1f42cf0ae4a1p-2, 0x1.82b358a3ab638p-3, 0x1.e1f2b30cd907bp-5,
+ 0x1.240f6d4071bd8p-6, 0x1.1522c9f3cd012p-8, 0x1.1fd0051a0525bp-10,
+ 0x1.9808a8b96c37ep-13, 0x1.b3f78e01152b5p-15, 0x1.49c85a7e1fd04p-18,
+ 0x1.471ca49184475p-19, -0x1.368f0b7ed9e36p-23, 0x1.882222f9049efp-23,
+ -0x1.a69ed2042842cp-25
+ };
- double m = z - 0x1.7p+1, i = __builtin_roundeven(m), step = __builtin_copysign(1.0,i);
- double d = m - i, d2 = d*d, d4 = d2*d2, d8 = d4*d4;
- double f = (c[0] + d*c[1]) + d2*(c[2] + d*c[3]) + d4*((c[4] + d*c[5]) + d2*(c[6] + d*c[7]))
- + d8*((c[8] + d*c[9]) + d2*(c[10] + d*c[11]) + d4*((c[12] + d*c[13]) + d2*(c[14] + d*c[15])));
- int jm = __builtin_fabs(i);
+ double m = z - 0x1.7p+1;
+ double i = roundeven (m);
+ double step = copysign (1.0, i);
+ double d = m - i, d2 = d * d, d4 = d2 * d2, d8 = d4 * d4;
+ double f = (c[0] + d * c[1]) + d2 * (c[2] + d * c[3])
+ + d4 * ((c[4] + d * c[5]) + d2 * (c[6] + d * c[7]))
+ + d8 * ((c[8] + d * c[9]) + d2 * (c[10] + d * c[11])
+ + d4 * ((c[12] + d * c[13]) + d2 * (c[14] + d * c[15])));
+ int jm = fabs (i);
double w = 1;
- if(jm){
- z -= 0.5 + step*0.5;
- w = z;
- for(int j=jm-1; j; j--) {z -= step; w *= z;}
- }
- if(i<=-0.5) w = 1/w;
+ if (jm)
+ {
+ z -= 0.5 + step * 0.5;
+ w = z;
+ for (int j = jm - 1; j; j--)
+ {
+ z -= step;
+ w *= z;
+ }
+ }
+ if (i <= -0.5)
+ w = 1 / w;
f *= w;
- b64u64_u rt = {.f = f};
+ uint64_t rt = asuint64 (f);
float r = f;
/* Deal with exceptional cases. */
- if(__builtin_expect(((rt.u+2)&0xfffffff) < 8, 0)){
- for(unsigned j=0;j<sizeof(tb)/sizeof(tb[0]);j++) {
- if(t.u==tb[j].u) return tb[j].f + tb[j].df;
+ if (__glibc_unlikely (((rt + 2) & 0xfffffff) < 8))
+ {
+ for (unsigned j = 0; j < sizeof (tb) / sizeof (tb[0]); j++)
+ if (t == tb[j].u)
+ return tb[j].f + tb[j].df;
}
- }
return r;
}
libm_alias_finite (__ieee754_gammaf_r, __gammaf_r)
deleted file mode 100644
@@ -1 +0,0 @@
-/* Not needed. */
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 <adhemerval.zanella@linaro.org> --- sysdeps/ieee754/flt-32/e_gammaf_r.c | 178 ++++++++++++++++------------ sysdeps/m68k/m680x0/fpu/math_errf.c | 1 - 2 files changed, 101 insertions(+), 78 deletions(-) delete mode 100644 sysdeps/m68k/m680x0/fpu/math_errf.c