summaryrefslogtreecommitdiff
path: root/src/math/fma.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/fma.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/fma.c')
-rw-r--r--src/math/fma.c38
1 files changed, 14 insertions, 24 deletions
diff --git a/src/math/fma.c b/src/math/fma.c
index 850a4a6c..1b1c1321 100644
--- a/src/math/fma.c
+++ b/src/math/fma.c
@@ -2,16 +2,6 @@
#include "libm.h"
#if LDBL_MANT_DIG==64 && LDBL_MAX_EXP==16384
-union ld80 {
- long double x;
- struct {
- uint64_t m;
- uint16_t e : 15;
- uint16_t s : 1;
- uint16_t pad;
- } bits;
-};
-
/* exact add, assumes exponent_x >= exponent_y */
static void add(long double *hi, long double *lo, long double x, long double y)
{
@@ -45,25 +35,25 @@ return an adjusted hi so that rounding it to double (or less) precision is corre
*/
static long double adjust(long double hi, long double lo)
{
- union ld80 uhi, ulo;
+ union ldshape uhi, ulo;
if (lo == 0)
return hi;
- uhi.x = hi;
- if (uhi.bits.m & 0x3ff)
+ uhi.f = hi;
+ if (uhi.i.m & 0x3ff)
return hi;
- ulo.x = lo;
- if (uhi.bits.s == ulo.bits.s)
- uhi.bits.m++;
+ ulo.f = lo;
+ if ((uhi.i.se & 0x8000) == (ulo.i.se & 0x8000))
+ uhi.i.m++;
else {
- uhi.bits.m--;
/* handle underflow and take care of ld80 implicit msb */
- if (uhi.bits.m == (uint64_t)-1/2) {
- uhi.bits.m *= 2;
- uhi.bits.e--;
+ if (uhi.i.m << 1 == 0) {
+ uhi.i.m = 0;
+ uhi.i.se--;
}
+ uhi.i.m--;
}
- return uhi.x;
+ return uhi.f;
}
/* adjusted add so the result is correct when rounded to double (or less) precision */
@@ -82,9 +72,9 @@ static long double dmul(long double x, long double y)
static int getexp(long double x)
{
- union ld80 u;
- u.x = x;
- return u.bits.e;
+ union ldshape u;
+ u.f = x;
+ return u.i.se & 0x7fff;
}
double fma(double x, double y, double z)