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)

diff --git a/src/math/acoshl.c b/src/math/acoshl.c
index 472c71cbc81749cd22c1886fd02ebfd77b5e628f..de31fb752b76262fcfaac7af9fcd6790179e2fc8 100644 (file)
--- 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 3ea887455925f1098dafb4cce417ecf056dacabe..e5f3175121cc6145d8d6dc7590d1ddb189065894 100644 (file)
--- 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 b4c5e58b97f2d5070073f3b418ba74edd9a89c18..72a1e345b88e383824267aefb0e7280e915debc9 100644 (file)
--- 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 d09070bbf7a4962d7f1f66826ef58cbaafce6158..080e5eb0425c0f8bd1ddcf1b0539c2ac9ba081d1 100644 (file)
--- 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 f9d90a6d3017ba1c1dc9659cefd82c78368c321d..3c1b553748eec19a18804362acd508d3877d7fb6 100644 (file)
--- a/src/math/frexpl.c
+++ b/src/math/frexpl.c
@@ -1,20 +1,20 @@
-#include <math.h>
-#include <stdint.h>
-#include <float.h>
-
-#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 14e3371b35aea410eddb2e40818155a26da48071..4864ddfa6a77c215a2dc1158ca1bb75137730063 100644 (file)
--- 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 0645cca0cf95cdf23afee4cd03b87ca61fa88f62..83a8f80c9984dc7b02e2de6d73f87ed58b566b8b 100644 (file)
--- 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 66559e9fb005886564f220624d1d463518426fb0..f594b85f3d6ce3a6722ebcf49eb6584410a41ad0 100644 (file)
--- 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;