diff options
| author | Szabolcs Nagy <nsz@port70.net> | 2013-09-04 15:52:23 +0000 |
|---|---|---|
| committer | Szabolcs Nagy <nsz@port70.net> | 2013-09-05 11:30:08 +0000 |
| commit | 34660d73bd0db29469d2758e1b48d2360edf3a2f (patch) | |
| tree | d803e3693b46e4bfe258c4d284fc5c3f7878daae /src/math/fma.c | |
| parent | 535104ab6a2d6f22098f79e7107963e3fc3448a3 (diff) | |
| download | musl-34660d73bd0db29469d2758e1b48d2360edf3a2f.tar.gz musl-34660d73bd0db29469d2758e1b48d2360edf3a2f.tar.xz musl-34660d73bd0db29469d2758e1b48d2360edf3a2f.zip | |
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.c | 38 |
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) |