summaryrefslogtreecommitdiff
path: root/src/math/fmal.c
diff options
context:
space:
mode:
authorSzabolcs Nagy <nsz@port70.net>2013-09-04 15:52:23 +0000
committerSzabolcs Nagy <nsz@port70.net>2013-09-05 11:30:08 +0000
commit34660d73bd0db29469d2758e1b48d2360edf3a2f (patch)
treed803e3693b46e4bfe258c4d284fc5c3f7878daae /src/math/fmal.c
parent535104ab6a2d6f22098f79e7107963e3fc3448a3 (diff)
downloadmusl-34660d73bd0db29469d2758e1b48d2360edf3a2f.tar.gz
math: fix remaining old long double code (erfl, fmal, lgammal, scalbnl)
in lgammal don't handle 1 and 2 specially, in fma use the new ldshape union instead of ld80 one.
Diffstat (limited to 'src/math/fmal.c')
-rw-r--r--src/math/fmal.c30
1 files changed, 16 insertions, 14 deletions
diff --git a/src/math/fmal.c b/src/math/fmal.c
index 87e30fcf..c68db255 100644
--- a/src/math/fmal.c
+++ b/src/math/fmal.c
@@ -34,6 +34,13 @@ long double fmal(long double x, long double y, long double z)
}
#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
#include <fenv.h>
+#if LDBL_MANT_DIG == 64
+#define LASTBIT(u) (u.i.m & 1)
+#define SPLIT (0x1p32L + 1)
+#elif LDBL_MANT_DIG == 113
+#define LASTBIT(u) (u.i.lo & 1)
+#define SPLIT (0x1p57L + 1)
+#endif
/*
* A struct dd represents a floating-point number with twice the precision
@@ -75,12 +82,12 @@ static inline struct dd dd_add(long double a, long double b)
static inline long double add_adjusted(long double a, long double b)
{
struct dd sum;
- union IEEEl2bits u;
+ union ldshape u;
sum = dd_add(a, b);
if (sum.lo != 0) {
- u.e = sum.hi;
- if ((u.bits.manl & 1) == 0)
+ u.f = sum.hi;
+ if (!LASTBIT(u))
sum.hi = nextafterl(sum.hi, INFINITY * sum.lo);
}
return (sum.hi);
@@ -95,7 +102,7 @@ static inline long double add_and_denormalize(long double a, long double b, int
{
struct dd sum;
int bits_lost;
- union IEEEl2bits u;
+ union ldshape u;
sum = dd_add(a, b);
@@ -110,9 +117,9 @@ static inline long double add_and_denormalize(long double a, long double b, int
* break the ties manually.
*/
if (sum.lo != 0) {
- u.e = sum.hi;
- bits_lost = -u.bits.exp - scale + 1;
- if (bits_lost != 1 ^ (int)(u.bits.manl & 1))
+ u.f = sum.hi;
+ bits_lost = -u.i.se - scale + 1;
+ if ((bits_lost != 1) ^ LASTBIT(u))
sum.hi = nextafterl(sum.hi, INFINITY * sum.lo);
}
return scalbnl(sum.hi, scale);
@@ -125,20 +132,15 @@ static inline long double add_and_denormalize(long double a, long double b, int
*/
static inline struct dd dd_mul(long double a, long double b)
{
-#if LDBL_MANT_DIG == 64
- static const long double split = 0x1p32L + 1.0;
-#elif LDBL_MANT_DIG == 113
- static const long double split = 0x1p57L + 1.0;
-#endif
struct dd ret;
long double ha, hb, la, lb, p, q;
- p = a * split;
+ p = a * SPLIT;
ha = a - p;
ha += p;
la = a - ha;
- p = b * split;
+ p = b * SPLIT;
hb = b - p;
hb += p;
lb = b - hb;