fmaq.c (fmaq): Merge from GLIBC.

2012-11-01  Tobias Burnus  <burnus@net-b.de>
            Joseph Myers  <joseph@codesourcery.com>

        * 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.


Co-Authored-By: Joseph Myers <joseph@codesourcery.com>

From-SVN: r193099
This commit is contained in:
Tobias Burnus 2012-11-02 17:59:30 +01:00 committed by Tobias Burnus
parent f2a1b46987
commit e4320d7da0
3 changed files with 56 additions and 21 deletions

View File

@ -1,3 +1,10 @@
2012-11-01 Tobias Burnus <burnus@net-b.de>
Joseph Myers <joseph@codesourcery.com>
* 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 <burnus@net-b.de>
* Makefile.am (libquadmath_la_SOURCES): Add new math/* files.

View File

@ -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)
{

View File

@ -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;