From aa0c4a204e28cfc90a9ee955691a6cbe014c9fde Mon Sep 17 00:00:00 2001 From: Szabolcs Nagy Date: Wed, 4 Sep 2013 15:54:02 +0000 Subject: math: long double fix (use ldshape union) * use new ldshape union consistently * add ld128 support to frexpl * simplify sqrtl comment (ld64 is not just arm) --- src/math/acoshl.c | 10 ++++------ src/math/asinhl.c | 5 +---- src/math/atanhl.c | 5 +---- src/math/coshl.c | 5 +---- src/math/frexpl.c | 36 ++++++++++++++---------------------- src/math/sinhl.c | 5 +---- src/math/sqrtl.c | 4 +--- src/math/tanhl.c | 5 +---- 8 files changed, 24 insertions(+), 51 deletions(-) diff --git a/src/math/acoshl.c b/src/math/acoshl.c index 472c71cb..de31fb75 100644 --- a/src/math/acoshl.c +++ b/src/math/acoshl.c @@ -9,15 +9,13 @@ long double acoshl(long double x) /* acosh(x) = log(x + sqrt(x*x-1)) */ long double acoshl(long double x) { - union { - long double f; - struct{uint64_t m; int16_t se; uint16_t pad;} i; - } u = {.f = x}; + union ldshape u = {x}; + int e = u.i.se & 0x7fff; - if (u.i.se < 0x3fff + 1) + if (e < 0x3fff + 1) /* x < 2, invalid if x < 1 or nan */ return log1pl(x-1 + sqrtl((x-1)*(x-1)+2*(x-1))); - if (u.i.se < 0x3fff + 32) + if (e < 0x3fff + 32) /* x < 0x1p32 */ return logl(2*x - 1/(x+sqrtl(x*x-1))); return logl(x) + 0.693147180559945309417232121458176568L; diff --git a/src/math/asinhl.c b/src/math/asinhl.c index 3ea88745..e5f31751 100644 --- a/src/math/asinhl.c +++ b/src/math/asinhl.c @@ -9,10 +9,7 @@ long double asinhl(long double x) /* asinh(x) = sign(x)*log(|x|+sqrt(x*x+1)) ~= x - x^3/6 + o(x^5) */ long double asinhl(long double x) { - union { - long double f; - struct{uint64_t m; uint16_t se; uint16_t pad;} i; - } u = {.f = x}; + union ldshape u = {x}; unsigned e = u.i.se & 0x7fff; unsigned s = u.i.se >> 15; diff --git a/src/math/atanhl.c b/src/math/atanhl.c index b4c5e58b..72a1e345 100644 --- a/src/math/atanhl.c +++ b/src/math/atanhl.c @@ -9,10 +9,7 @@ long double atanhl(long double x) /* atanh(x) = log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2 ~= x + x^3/3 + o(x^5) */ long double atanhl(long double x) { - union { - long double f; - struct{uint64_t m; uint16_t se; uint16_t pad;} i; - } u = {.f = x}; + union ldshape u = {x}; unsigned e = u.i.se & 0x7fff; unsigned s = u.i.se >> 15; diff --git a/src/math/coshl.c b/src/math/coshl.c index d09070bb..080e5eb0 100644 --- a/src/math/coshl.c +++ b/src/math/coshl.c @@ -8,10 +8,7 @@ long double coshl(long double x) #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 long double coshl(long double x) { - union { - long double f; - struct{uint64_t m; uint16_t se; uint16_t pad;} i; - } u = {.f = x}; + union ldshape u = {x}; unsigned ex = u.i.se & 0x7fff; uint32_t w; long double t; diff --git a/src/math/frexpl.c b/src/math/frexpl.c index f9d90a6d..3c1b5537 100644 --- a/src/math/frexpl.c +++ b/src/math/frexpl.c @@ -1,20 +1,20 @@ -#include -#include -#include - -#if LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 - -/* This version is for 80-bit little endian long double */ +#include "libm.h" +#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 long double frexpl(long double x, int *e) { - union { long double ld; uint16_t hw[5]; } y = { x }; - int ee = y.hw[4]&0x7fff; + return frexp(x, e); +} +#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 +long double frexpl(long double x, int *e) +{ + union ldshape u = {x}; + int ee = u.i.se & 0x7fff; if (!ee) { if (x) { - x = frexpl(x*0x1p64, e); - *e -= 64; + x = frexpl(x*0x1p120, e); + *e -= 120; } else *e = 0; return x; } else if (ee == 0x7fff) { @@ -22,16 +22,8 @@ long double frexpl(long double x, int *e) } *e = ee - 0x3ffe; - y.hw[4] &= 0x8000; - y.hw[4] |= 0x3ffe; - return y.ld; + u.i.se &= 0x8000; + u.i.se |= 0x3ffe; + return u.f; } - -#else - -long double frexpl(long double x, int *e) -{ - return frexp(x, e); -} - #endif diff --git a/src/math/sinhl.c b/src/math/sinhl.c index 14e3371b..4864ddfa 100644 --- a/src/math/sinhl.c +++ b/src/math/sinhl.c @@ -8,10 +8,7 @@ long double sinhl(long double x) #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 long double sinhl(long double x) { - union { - long double f; - struct{uint64_t m; uint16_t se; uint16_t pad;} i; - } u = {.f = x}; + union ldshape u = {x}; unsigned ex = u.i.se & 0x7fff; long double h, t, absx; diff --git a/src/math/sqrtl.c b/src/math/sqrtl.c index 0645cca0..83a8f80c 100644 --- a/src/math/sqrtl.c +++ b/src/math/sqrtl.c @@ -2,8 +2,6 @@ long double sqrtl(long double x) { - /* FIXME: implement sqrtl in C. At least this works for now on - * ARM (which uses ld64), the only arch without sqrtl asm - * that's supported so far. */ + /* FIXME: implement in C, this is for LDBL_MANT_DIG == 64 only */ return sqrt(x); } diff --git a/src/math/tanhl.c b/src/math/tanhl.c index 66559e9f..f594b85f 100644 --- a/src/math/tanhl.c +++ b/src/math/tanhl.c @@ -8,10 +8,7 @@ long double tanhl(long double x) #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 long double tanhl(long double x) { - union { - long double f; - struct{uint64_t m; uint16_t se; uint16_t pad;} i; - } u = {.f = x}; + union ldshape u = {x}; unsigned ex = u.i.se & 0x7fff; unsigned sign = u.i.se & 0x8000; uint32_t w; -- cgit v1.2.1