summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorSzabolcs Nagy <nsz@port70.net>2013-09-04 15:54:02 +0000
committerSzabolcs Nagy <nsz@port70.net>2013-09-05 11:30:08 +0000
commitaa0c4a204e28cfc90a9ee955691a6cbe014c9fde (patch)
treee75c43e934d087e203dc999ba75168ee78e6ee7b
parent2eaed464e2080d8321d3903b71086a1ecfc4ee4a (diff)
downloadmusl-aa0c4a204e28cfc90a9ee955691a6cbe014c9fde.tar.gz
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)
-rw-r--r--src/math/acoshl.c10
-rw-r--r--src/math/asinhl.c5
-rw-r--r--src/math/atanhl.c5
-rw-r--r--src/math/coshl.c5
-rw-r--r--src/math/frexpl.c36
-rw-r--r--src/math/sinhl.c5
-rw-r--r--src/math/sqrtl.c4
-rw-r--r--src/math/tanhl.c5
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;