diff --git a/libquadmath/ChangeLog b/libquadmath/ChangeLog index efb95d179d2..f1b4ce9b816 100644 --- a/libquadmath/ChangeLog +++ b/libquadmath/ChangeLog @@ -1,3 +1,10 @@ +2012-11-01 Tobias Burnus + Joseph Myers + + * math/fmaq.c (fmaq): Merge from GLIBC. Handle cases + with small x * y using scaling, not as x * y + z. + * math/lgammaq.c (lgammaq): Fix signgam handling. + 2012-11-01 Tobias Burnus * Makefile.am (libquadmath_la_SOURCES): Add new math/* files. diff --git a/libquadmath/math/fmaq.c b/libquadmath/math/fmaq.c index 5616c1a2781..e5a9d37627d 100644 --- a/libquadmath/math/fmaq.c +++ b/libquadmath/math/fmaq.c @@ -73,6 +73,37 @@ fmaq (__float128 x, __float128 y, __float128 z) || u.ieee.exponent + v.ieee.exponent < IEEE854_FLOAT128_BIAS - FLT128_MANT_DIG - 2) return x * y + z; + /* If x * y is less than 1/4 of FLT128_DENORM_MIN, neither the + result nor whether there is underflow depends on its exact + value, only on its sign. */ + if (u.ieee.exponent + v.ieee.exponent + < IEEE854_FLT128_DOUBLE_BIAS - FLT128_MANT_DIG - 2) + { + int neg = u.ieee.negative ^ v.ieee.negative; + __float128 tiny = neg ? -0x1p-16494L : 0x1p-16494L; + if (w.ieee.exponent >= 3) + return tiny + z; + /* Scaling up, adding TINY and scaling down produces the + correct result, because in round-to-nearest mode adding + TINY has no effect and in other modes double rounding is + harmless. But it may not produce required underflow + exceptions. */ + v.value = z * 0x1p114L + tiny; + if (TININESS_AFTER_ROUNDING + ? v.ieee.exponent < 115 + : (w.ieee.exponent == 0 + || (w.ieee.exponent == 1 + && w.ieee.negative != neg + && w.ieee.mantissa3 == 0 + && w.ieee.mantissa2 == 0 + && w.ieee.mantissa1 == 0 + && w.ieee.mantissa0 == 0))) + { + volatile __float128 force_underflow = x * y; + (void) force_underflow; + } + return v.value * 0x1p-114L; + } if (u.ieee.exponent + v.ieee.exponent >= 0x7fff + IEEE854_FLOAT128_BIAS - FLT128_MANT_DIG) { @@ -228,17 +259,17 @@ fmaq (__float128 x, __float128 y, __float128 z) for proper rounding. */ if (v.ieee.exponent == 226) { - /* If the exponent would be in the normal range when - rounding to normal precision with unbounded exponent - range, the exact result is known and spurious underflows - must be avoided on systems detecting tininess after - rounding. */ - if (TININESS_AFTER_ROUNDING) - { - w.value = a1 + u.value; - if (w.ieee.exponent == 227) - return w.value * 0x1p-226L; - } + /* If the exponent would be in the normal range when + rounding to normal precision with unbounded exponent + range, the exact result is known and spurious underflows + must be avoided on systems detecting tininess after + rounding. */ + if (TININESS_AFTER_ROUNDING) + { + w.value = a1 + u.value; + if (w.ieee.exponent == 227) + return w.value * 0x1p-226L; + } /* v.ieee.mant_low & 2 is LSB bit of the result before rounding, v.ieee.mant_low & 1 is the round bit and j is our sticky bit. */ diff --git a/libquadmath/math/lgammaq.c b/libquadmath/math/lgammaq.c index 361f7037bc3..485939af1b4 100644 --- a/libquadmath/math/lgammaq.c +++ b/libquadmath/math/lgammaq.c @@ -759,7 +759,8 @@ lgammaq (__float128 x) { __float128 p, q, w, z, nx; int i, nn; - int sign; + + signgam = 1; if (! finiteq (x)) return x * x; @@ -767,11 +768,9 @@ lgammaq (__float128 x) if (x == 0.0Q) { if (signbitq (x)) - sign = -1; + signgam = -1; } - signgam = sign; - if (x < 0.0Q) { q = -x; @@ -780,9 +779,9 @@ lgammaq (__float128 x) return (one / (p - p)); i = p; if ((i & 1) == 0) - sign = -1; + signgam = -1; else - sign = 1; + signgam = 1; z = q - p; if (z > 0.5Q) { @@ -790,10 +789,8 @@ lgammaq (__float128 x) z = p - q; } z = q * sinq (PIQ * z); - signgam = sign; - if (z == 0.0Q) - return (sign * huge * huge); + return (signgam * huge * huge); w = lgammaq (q); z = logq (PIQ / z) - w; return (z); @@ -1025,7 +1022,7 @@ lgammaq (__float128 x) } if (x > MAXLGM) - return (sign * huge * huge); + return (signgam * huge * huge); q = ls2pi - x; q = (x - 0.5Q) * logq (x) + q;