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)

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 <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 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;