diff options
-rw-r--r-- | src/internal/libm.h | 28 | ||||
-rw-r--r-- | src/internal/longdbl.h | 26 | ||||
-rw-r--r-- | src/math/__fpclassifyl.c | 28 | ||||
-rw-r--r-- | src/math/__signbitl.c | 4 | ||||
-rw-r--r-- | src/math/copysignl.c | 6 | ||||
-rw-r--r-- | src/math/fabsl.c | 4 | ||||
-rw-r--r-- | src/math/ilogbl.c | 8 | ||||
-rw-r--r-- | src/math/nextafterl.c | 81 |
8 files changed, 89 insertions, 96 deletions
diff --git a/src/internal/libm.h b/src/internal/libm.h index 64cc8388..52ac61ec 100644 --- a/src/internal/libm.h +++ b/src/internal/libm.h @@ -17,6 +17,7 @@ #include <float.h> #include <math.h> #include <complex.h> +#include <endian.h> #include "longdbl.h" @@ -32,6 +33,33 @@ union dshape { uint64_t bits; }; +#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 +#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 && __BYTE_ORDER == __LITTLE_ENDIAN +union ldshape { + long double f; + struct { + uint64_t m; + uint16_t se; + } i; +}; +#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 && __BYTE_ORDER == __LITTLE_ENDIAN +union ldshape { + long double f; + struct { + uint64_t lo; + uint32_t mid; + uint16_t top; + uint16_t se; + } i; + struct { + uint64_t lo; + uint64_t hi; + } i2; +}; +#else +#error Unsupported long double representation +#endif + #define FORCE_EVAL(x) do { \ if (sizeof(x) == sizeof(float)) { \ volatile float __x; \ diff --git a/src/internal/longdbl.h b/src/internal/longdbl.h index 25ec8021..e93fb4ff 100644 --- a/src/internal/longdbl.h +++ b/src/internal/longdbl.h @@ -4,32 +4,6 @@ #include <float.h> #include <stdint.h> -#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 -#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 -union ldshape { - long double value; - struct { - uint64_t m; - uint16_t exp:15; - uint16_t sign:1; - uint16_t pad; - } bits; -}; -#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 -union ldshape { - long double value; - struct { - uint64_t mlo; - uint64_t mhi:48; - uint16_t exp:15; - uint16_t sign:1; - } bits; -}; -#else -#error Unsupported long double representation -#endif - - // FIXME: hacks to make freebsd+openbsd long double code happy // union and macros for freebsd diff --git a/src/math/__fpclassifyl.c b/src/math/__fpclassifyl.c index e4d231b5..6365c588 100644 --- a/src/math/__fpclassifyl.c +++ b/src/math/__fpclassifyl.c @@ -3,28 +3,28 @@ #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 +/* invalid representations (top bit of u.i.m is wrong) are not handled */ int __fpclassifyl(long double x) { - union ldshape u = { x }; - int e = u.bits.exp; - if (!e) { - if (u.bits.m >> 63) return FP_NAN; - else if (u.bits.m) return FP_SUBNORMAL; - else return FP_ZERO; - } + union ldshape u = {x}; + int e = u.i.se & 0x7fff; + if (!e) + return u.i.m ? FP_SUBNORMAL : FP_ZERO; if (e == 0x7fff) - return u.bits.m & (uint64_t)-1>>1 ? FP_NAN : FP_INFINITE; - return u.bits.m & (uint64_t)1<<63 ? FP_NORMAL : FP_NAN; + return u.i.m << 1 ? FP_NAN : FP_INFINITE; + return FP_NORMAL; } #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 int __fpclassifyl(long double x) { - union ldshape u = { x }; - int e = u.bits.exp; + union ldshape u = {x}; + int e = u.i.se & 0x7fff; if (!e) - return u.bits.mlo | u.bits.mhi ? FP_SUBNORMAL : FP_ZERO; - if (e == 0x7fff) - return u.bits.mlo | u.bits.mhi ? FP_NAN : FP_INFINITE; + return u.i2.lo | u.i2.hi ? FP_SUBNORMAL : FP_ZERO; + if (e == 0x7fff) { + u.i.se = 0; + return u.i2.lo | u.i2.hi ? FP_NAN : FP_INFINITE; + } return FP_NORMAL; } #endif diff --git a/src/math/__signbitl.c b/src/math/__signbitl.c index 81adb6ce..c52c87bb 100644 --- a/src/math/__signbitl.c +++ b/src/math/__signbitl.c @@ -1,11 +1,9 @@ #include "libm.h" -// FIXME: should be a macro #if (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 int __signbitl(long double x) { union ldshape u = {x}; - - return u.bits.sign; + return u.i.se >> 15; } #endif diff --git a/src/math/copysignl.c b/src/math/copysignl.c index 72a21488..9dd933cf 100644 --- a/src/math/copysignl.c +++ b/src/math/copysignl.c @@ -9,8 +9,8 @@ long double copysignl(long double x, long double y) long double copysignl(long double x, long double y) { union ldshape ux = {x}, uy = {y}; - - ux.bits.sign = uy.bits.sign; - return ux.value; + ux.i.se &= 0x7fff; + ux.i.se |= uy.i.se & 0x8000; + return ux.f; } #endif diff --git a/src/math/fabsl.c b/src/math/fabsl.c index 711d908a..c4f36ec2 100644 --- a/src/math/fabsl.c +++ b/src/math/fabsl.c @@ -9,7 +9,7 @@ long double fabsl(long double x) { union ldshape u = {x}; - u.bits.sign = 0; - return u.value; + u.i.se &= 0x7fff; + return u.f; } #endif diff --git a/src/math/ilogbl.c b/src/math/ilogbl.c index 1512934f..7df6eb6c 100644 --- a/src/math/ilogbl.c +++ b/src/math/ilogbl.c @@ -10,12 +10,12 @@ int ilogbl(long double x) int ilogbl(long double x) { union ldshape u = {x}; - uint64_t m = u.bits.m; - int e = u.bits.exp; + uint64_t m = u.i.m; + int e = u.i.se & 0x7fff; if (!e) { if (m == 0) { - FORCE_EVAL(0/0.0f); + FORCE_EVAL(x/x); return FP_ILOGB0; } /* subnormal x */ @@ -25,7 +25,7 @@ int ilogbl(long double x) if (e == 0x7fff) { FORCE_EVAL(0/0.0f); /* in ld80 msb is set in inf */ - return m & (uint64_t)-1>>1 ? FP_ILOGBNAN : INT_MAX; + return m << 1 ? FP_ILOGBNAN : INT_MAX; } return e - 0x3fff; } diff --git a/src/math/nextafterl.c b/src/math/nextafterl.c index edc3cc9c..37e858fb 100644 --- a/src/math/nextafterl.c +++ b/src/math/nextafterl.c @@ -6,7 +6,6 @@ long double nextafterl(long double x, long double y) return nextafter(x, y); } #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 -#define MSB ((uint64_t)1<<63) long double nextafterl(long double x, long double y) { union ldshape ux, uy; @@ -15,32 +14,32 @@ long double nextafterl(long double x, long double y) return x + y; if (x == y) return y; - ux.value = x; + ux.f = x; if (x == 0) { - uy.value = y; - ux.bits.m = 1; - ux.bits.sign = uy.bits.sign; - } else if (x < y ^ ux.bits.sign) { - ux.bits.m++; - if ((ux.bits.m & ~MSB) == 0) { - ux.bits.m = MSB; - ux.bits.exp++; + uy.f = y; + ux.i.m = 1; + ux.i.se = uy.i.se & 0x8000; + } else if ((x < y) == !(ux.i.se & 0x8000)) { + ux.i.m++; + if (ux.i.m << 1 == 0) { + ux.i.m = 1ULL << 63; + ux.i.se++; } } else { - if ((ux.bits.m & ~MSB) == 0) { - ux.bits.exp--; - if (ux.bits.exp) - ux.bits.m = 0; + if (ux.i.m << 1 == 0) { + ux.i.se--; + if (ux.i.se) + ux.i.m = 0; } - ux.bits.m--; + ux.i.m--; } - /* raise overflow if ux.value is infinite and x is finite */ - if (ux.bits.exp == 0x7fff) + /* raise overflow if ux is infinite and x is finite */ + if ((ux.i.se & 0x7fff) == 0x7fff) return x + x; - /* raise underflow if ux.value is subnormal or zero */ - if (ux.bits.exp == 0) - FORCE_EVAL(x*x + ux.value*ux.value); - return ux.value; + /* raise underflow if ux is subnormal or zero */ + if ((ux.i.se & 0x7fff) == 0) + FORCE_EVAL(x*x + ux.f*ux.f); + return ux.f; } #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 long double nextafterl(long double x, long double y) @@ -51,32 +50,26 @@ long double nextafterl(long double x, long double y) return x + y; if (x == y) return y; - ux.value = x; + ux.f = x; if (x == 0) { - uy.value = y; - ux.bits.mlo = 1; - ux.bits.sign = uy.bits.sign; - } else if (x < y ^ ux.bits.sign) { - ux.bits.mlo++; - if (ux.bits.mlo == 0) { - ux.bits.mhi++; - if (ux.bits.mhi == 0) - ux.bits.exp++; - } + uy.f = y; + ux.i.lo = 1; + ux.i.se = uy.i.se & 0x8000; + } else if ((x < y) == !(ux.i.se & 0x8000)) { + ux.i2.lo++; + if (ux.i2.lo == 0) + ux.i2.hi++; } else { - if (ux.bits.mlo == 0) { - if (ux.bits.mhi == 0) - ux.bits.exp--; - ux.bits.mhi--; - } - ux.bits.mlo--; + if (ux.i2.lo == 0) + ux.i2.hi--; + ux.i2.lo--; } - /* raise overflow if ux.value is infinite and x is finite */ - if (ux.bits.exp == 0x7fff) + /* raise overflow if ux is infinite and x is finite */ + if ((ux.i.se & 0x7fff) == 0x7fff) return x + x; - /* raise underflow if ux.value is subnormal or zero */ - if (ux.bits.exp == 0) - FORCE_EVAL(x*x + ux.value*ux.value); - return ux.value; + /* raise underflow if ux is subnormal or zero */ + if ((ux.i.se & 0x7fff) == 0) + FORCE_EVAL(x*x + ux.f*ux.f); + return ux.f; } #endif |