mirror of
https://sourceware.org/git/glibc.git
synced 2024-11-25 06:20:06 +00:00
Implement proper fmal for ldbl-128ibm (bug 13304).
ldbl-128ibm had an implementation of fmal that just did (x * y) + z in most cases, with no attempt at actually being a fused operation. This patch replaces it with a genuine fused operation. It is not necessarily correctly rounding, but should produce a result at least as accurate as the long double arithmetic operations in libgcc, which I think is all that can reasonably be expected for such a non-IEEE format where arithmetic is approximate rather than rounded according to any particular rule for determining the exact result. Like the libgcc arithmetic, it may produce spurious overflow and underflow results, and it falls back to the libgcc multiplication in the case of (finite, finite, zero). This concludes the fixes for bug 13304; any subsequently found fma issues should go in separate Bugzilla bugs. Various other pieces of bug 13304 were fixed in past releases over the past several years. Tested for powerpc. [BZ #13304] * sysdeps/ieee754/ldbl-128ibm/s_fmal.c: Include <fenv.h>, <float.h>, <math_private.h> and <stdlib.h>. (add_split): New function. (mul_split): Likewise. (ext_val): New typedef. (store_ext_val): New function. (mul_ext_val): New function. (compare): New function. (add_split_ext): New function. (__fmal): After checking for Inf, NaN and zero, compute result as an exact sum of scaled double values in round-to-nearest before adding those up and adjusting for other rounding modes. * math/auto-libm-test-in: Remove xfail-rounding:ldbl-128ibm from tests of fma. * math/auto-libm-test-out: Regenerated.
This commit is contained in:
Joseph Myers
parent
de71e0421b
commit
ffe9aaf2b9
19
ChangeLog
19
ChangeLog
@ -1,3 +1,22 @@
|
||||
2016-05-19 Joseph Myers <joseph@codesourcery.com>
|
||||
|
||||
[BZ #13304]
|
||||
* sysdeps/ieee754/ldbl-128ibm/s_fmal.c: Include <fenv.h>,
|
||||
<float.h>, <math_private.h> and <stdlib.h>.
|
||||
(add_split): New function.
|
||||
(mul_split): Likewise.
|
||||
(ext_val): New typedef.
|
||||
(store_ext_val): New function.
|
||||
(mul_ext_val): New function.
|
||||
(compare): New function.
|
||||
(add_split_ext): New function.
|
||||
(__fmal): After checking for Inf, NaN and zero, compute result as
|
||||
an exact sum of scaled double values in round-to-nearest before
|
||||
adding those up and adjusting for other rounding modes.
|
||||
* math/auto-libm-test-in: Remove xfail-rounding:ldbl-128ibm from
|
||||
tests of fma.
|
||||
* math/auto-libm-test-out: Regenerated.
|
||||
|
||||
2016-05-19 H.J. Lu <hongjiu.lu@intel.com>
|
||||
|
||||
[BZ #20119]
|
||||
|
||||
@ -2069,15 +2069,14 @@ fma -min -min 0 missing-errno
|
||||
fma -min -min -0 missing-errno
|
||||
|
||||
# Bug 6801: errno setting may be missing.
|
||||
# Bug 13304: results on directed rounding may be incorrect.
|
||||
fma max max min missing-errno xfail-rounding:ldbl-128ibm
|
||||
fma max max -min missing-errno xfail-rounding:ldbl-128ibm
|
||||
fma max -max min missing-errno xfail-rounding:ldbl-128ibm
|
||||
fma max -max -min missing-errno xfail-rounding:ldbl-128ibm
|
||||
fma -max max min missing-errno xfail-rounding:ldbl-128ibm
|
||||
fma -max max -min missing-errno xfail-rounding:ldbl-128ibm
|
||||
fma -max -max min missing-errno xfail-rounding:ldbl-128ibm
|
||||
fma -max -max -min missing-errno xfail-rounding:ldbl-128ibm
|
||||
fma max max min missing-errno
|
||||
fma max max -min missing-errno
|
||||
fma max -max min missing-errno
|
||||
fma max -max -min missing-errno
|
||||
fma -max max min missing-errno
|
||||
fma -max max -min missing-errno
|
||||
fma -max -max min missing-errno
|
||||
fma -max -max -min missing-errno
|
||||
|
||||
fma 0x1.7ff8p+13 0x1.000002p+0 0x1.ffffp-24
|
||||
fma 0x1.fffp+0 0x1.00001p+0 -0x1.fffp+0
|
||||
|
||||
12256
math/auto-libm-test-out
12256
math/auto-libm-test-out
File diff suppressed because it is too large
Load Diff
@ -17,25 +17,263 @@
|
||||
License along with the GNU C Library; if not, see
|
||||
<http://www.gnu.org/licenses/>. */
|
||||
|
||||
#include <fenv.h>
|
||||
#include <float.h>
|
||||
#include <math.h>
|
||||
#include <math_private.h>
|
||||
#include <math_ldbl_opt.h>
|
||||
#include <stdlib.h>
|
||||
|
||||
/* Calculate X + Y exactly and store the result in *HI + *LO. It is
|
||||
given that |X| >= |Y| and the values are small enough that no
|
||||
overflow occurs. */
|
||||
|
||||
static void
|
||||
add_split (double *hi, double *lo, double x, double y)
|
||||
{
|
||||
/* Apply Dekker's algorithm. */
|
||||
*hi = x + y;
|
||||
*lo = (x - *hi) + y;
|
||||
}
|
||||
|
||||
/* Calculate X * Y exactly and store the result in *HI + *LO. It is
|
||||
given that the values are small enough that no overflow occurs and
|
||||
large enough (or zero) that no underflow occurs. */
|
||||
|
||||
static void
|
||||
mul_split (double *hi, double *lo, double x, double y)
|
||||
{
|
||||
#ifdef __FP_FAST_FMA
|
||||
/* Fast built-in fused multiply-add. */
|
||||
*hi = x * y;
|
||||
*lo = __builtin_fma (x, y, -*hi);
|
||||
#else
|
||||
/* Apply Dekker's algorithm. */
|
||||
*hi = x * y;
|
||||
# define C ((1 << (DBL_MANT_DIG + 1) / 2) + 1)
|
||||
double x1 = x * C;
|
||||
double y1 = y * C;
|
||||
# undef C
|
||||
x1 = (x - x1) + x1;
|
||||
y1 = (y - y1) + y1;
|
||||
double x2 = x - x1;
|
||||
double y2 = y - y1;
|
||||
*lo = (((x1 * y1 - *hi) + x1 * y2) + x2 * y1) + x2 * y2;
|
||||
#endif
|
||||
}
|
||||
|
||||
/* Value with extended range, used in intermediate computations. */
|
||||
typedef struct
|
||||
{
|
||||
/* Value in [0.5, 1), as from frexp, or 0. */
|
||||
double val;
|
||||
/* Exponent of power of 2 it is multiplied by, or 0 for zero. */
|
||||
int exp;
|
||||
} ext_val;
|
||||
|
||||
/* Store D as an ext_val value. */
|
||||
|
||||
static void
|
||||
store_ext_val (ext_val *v, double d)
|
||||
{
|
||||
v->val = __frexp (d, &v->exp);
|
||||
}
|
||||
|
||||
/* Store X * Y as ext_val values *V0 and *V1. */
|
||||
|
||||
static void
|
||||
mul_ext_val (ext_val *v0, ext_val *v1, double x, double y)
|
||||
{
|
||||
int xexp, yexp;
|
||||
x = __frexp (x, &xexp);
|
||||
y = __frexp (y, &yexp);
|
||||
double hi, lo;
|
||||
mul_split (&hi, &lo, x, y);
|
||||
store_ext_val (v0, hi);
|
||||
if (hi != 0)
|
||||
v0->exp += xexp + yexp;
|
||||
store_ext_val (v1, lo);
|
||||
if (lo != 0)
|
||||
v1->exp += xexp + yexp;
|
||||
}
|
||||
|
||||
/* Compare absolute values of ext_val values pointed to by P and Q for
|
||||
qsort. */
|
||||
|
||||
static int
|
||||
compare (const void *p, const void *q)
|
||||
{
|
||||
const ext_val *pe = p;
|
||||
const ext_val *qe = q;
|
||||
if (pe->val == 0)
|
||||
return qe->val == 0 ? 0 : -1;
|
||||
else if (qe->val == 0)
|
||||
return 1;
|
||||
else if (pe->exp < qe->exp)
|
||||
return -1;
|
||||
else if (pe->exp > qe->exp)
|
||||
return 1;
|
||||
else
|
||||
{
|
||||
double pd = fabs (pe->val);
|
||||
double qd = fabs (qe->val);
|
||||
if (pd < qd)
|
||||
return -1;
|
||||
else if (pd == qd)
|
||||
return 0;
|
||||
else
|
||||
return 1;
|
||||
}
|
||||
}
|
||||
|
||||
/* Calculate *X + *Y exactly, storing the high part in *X (rounded to
|
||||
nearest) and the low part in *Y. It is given that |X| >= |Y|. */
|
||||
|
||||
static void
|
||||
add_split_ext (ext_val *x, ext_val *y)
|
||||
{
|
||||
int xexp = x->exp, yexp = y->exp;
|
||||
if (y->val == 0 || xexp - yexp > 53)
|
||||
return;
|
||||
double hi = x->val;
|
||||
double lo = __scalbn (y->val, yexp - xexp);
|
||||
add_split (&hi, &lo, hi, lo);
|
||||
store_ext_val (x, hi);
|
||||
if (hi != 0)
|
||||
x->exp += xexp;
|
||||
store_ext_val (y, lo);
|
||||
if (lo != 0)
|
||||
y->exp += xexp;
|
||||
}
|
||||
|
||||
long double
|
||||
__fmal (long double x, long double y, long double z)
|
||||
{
|
||||
/* An IBM long double 128 is really just 2 IEEE64 doubles, and in
|
||||
* the case of inf/nan only the first double counts. So we use the
|
||||
* (double) cast to avoid any data movement. */
|
||||
if ((isfinite ((double)x) && isfinite ((double)y)) && isinf ((double)z))
|
||||
return (z);
|
||||
double xhi, xlo, yhi, ylo, zhi, zlo;
|
||||
int64_t hx, hy, hz;
|
||||
int xexp, yexp, zexp;
|
||||
double scale_val;
|
||||
int scale_exp;
|
||||
ldbl_unpack (x, &xhi, &xlo);
|
||||
EXTRACT_WORDS64 (hx, xhi);
|
||||
xexp = (hx & 0x7ff0000000000000LL) >> 52;
|
||||
ldbl_unpack (y, &yhi, &ylo);
|
||||
EXTRACT_WORDS64 (hy, yhi);
|
||||
yexp = (hy & 0x7ff0000000000000LL) >> 52;
|
||||
ldbl_unpack (z, &zhi, &zlo);
|
||||
EXTRACT_WORDS64 (hz, zhi);
|
||||
zexp = (hz & 0x7ff0000000000000LL) >> 52;
|
||||
|
||||
/* If z is zero and x are y are nonzero, compute the result
|
||||
as x * y to avoid the wrong sign of a zero result if x * y
|
||||
underflows to 0. */
|
||||
/* If z is Inf or NaN, but x and y are finite, avoid any exceptions
|
||||
from computing x * y. */
|
||||
if (zexp == 0x7ff && xexp != 0x7ff && yexp != 0x7ff)
|
||||
return (z + x) + y;
|
||||
|
||||
/* If z is zero and x are y are nonzero, compute the result as x * y
|
||||
to avoid the wrong sign of a zero result if x * y underflows to
|
||||
0. */
|
||||
if (z == 0 && x != 0 && y != 0)
|
||||
return x * y;
|
||||
|
||||
/* If x or y or z is Inf/NaN, or if x * y is zero, compute as x * y
|
||||
+ z. */
|
||||
if (xexp == 0x7ff || yexp == 0x7ff || zexp == 0x7ff
|
||||
|| x == 0 || y == 0)
|
||||
return (x * y) + z;
|
||||
|
||||
{
|
||||
SET_RESTORE_ROUND (FE_TONEAREST);
|
||||
|
||||
ext_val vals[10];
|
||||
store_ext_val (&vals[0], zhi);
|
||||
store_ext_val (&vals[1], zlo);
|
||||
mul_ext_val (&vals[2], &vals[3], xhi, yhi);
|
||||
mul_ext_val (&vals[4], &vals[5], xhi, ylo);
|
||||
mul_ext_val (&vals[6], &vals[7], xlo, yhi);
|
||||
mul_ext_val (&vals[8], &vals[9], xlo, ylo);
|
||||
qsort (vals, 10, sizeof (ext_val), compare);
|
||||
/* Add up the values so that each element of VALS has absolute
|
||||
value at most equal to the last set bit of the next nonzero
|
||||
element. */
|
||||
for (size_t i = 0; i <= 8; i++)
|
||||
{
|
||||
add_split_ext (&vals[i + 1], &vals[i]);
|
||||
qsort (vals + i + 1, 9 - i, sizeof (ext_val), compare);
|
||||
}
|
||||
/* Add up the values in the other direction, so that each element
|
||||
of VALS has absolute value less than 5ulp of the next
|
||||
value. */
|
||||
size_t dstpos = 9;
|
||||
for (size_t i = 1; i <= 9; i++)
|
||||
{
|
||||
if (vals[dstpos].val == 0)
|
||||
{
|
||||
vals[dstpos] = vals[9 - i];
|
||||
vals[9 - i].val = 0;
|
||||
vals[9 - i].exp = 0;
|
||||
}
|
||||
else
|
||||
{
|
||||
add_split_ext (&vals[dstpos], &vals[9 - i]);
|
||||
if (vals[9 - i].val != 0)
|
||||
{
|
||||
if (9 - i < dstpos - 1)
|
||||
{
|
||||
vals[dstpos - 1] = vals[9 - i];
|
||||
vals[9 - i].val = 0;
|
||||
vals[9 - i].exp = 0;
|
||||
}
|
||||
dstpos--;
|
||||
}
|
||||
}
|
||||
}
|
||||
/* If the result is an exact zero, it results from adding two
|
||||
values with opposite signs; recompute in the original rounding
|
||||
mode. */
|
||||
if (vals[9].val == 0)
|
||||
goto zero_out;
|
||||
/* Adding the top three values will now give a result as accurate
|
||||
as the underlying long double arithmetic. */
|
||||
add_split_ext (&vals[9], &vals[8]);
|
||||
if (compare (&vals[8], &vals[7]) < 0)
|
||||
{
|
||||
ext_val tmp = vals[7];
|
||||
vals[7] = vals[8];
|
||||
vals[8] = tmp;
|
||||
}
|
||||
add_split_ext (&vals[8], &vals[7]);
|
||||
add_split_ext (&vals[9], &vals[8]);
|
||||
if (vals[9].exp > DBL_MAX_EXP || vals[9].exp < DBL_MIN_EXP)
|
||||
{
|
||||
/* Overflow or underflow, with the result depending on the
|
||||
original rounding mode, but not on the low part computed
|
||||
here. */
|
||||
scale_val = vals[9].val;
|
||||
scale_exp = vals[9].exp;
|
||||
goto scale_out;
|
||||
}
|
||||
double hi = __scalbn (vals[9].val, vals[9].exp);
|
||||
double lo = __scalbn (vals[8].val, vals[8].exp);
|
||||
/* It is possible that the low part became subnormal and was
|
||||
rounded so that the result is no longer canonical. */
|
||||
ldbl_canonicalize (&hi, &lo);
|
||||
long double ret = ldbl_pack (hi, lo);
|
||||
math_check_force_underflow (ret);
|
||||
return ret;
|
||||
}
|
||||
|
||||
scale_out:
|
||||
scale_val = math_opt_barrier (scale_val);
|
||||
scale_val = __scalbn (scale_val, scale_exp);
|
||||
if (fabs (scale_val) == DBL_MAX)
|
||||
return __copysignl (LDBL_MAX, scale_val);
|
||||
math_check_force_underflow (scale_val);
|
||||
return scale_val;
|
||||
|
||||
zero_out:;
|
||||
double zero = 0.0;
|
||||
zero = math_opt_barrier (zero);
|
||||
return zero - zero;
|
||||
}
|
||||
#if IS_IN (libm)
|
||||
long_double_symbol (libm, __fmal, fmal);
|
||||
|
||||
Loading…
Reference in New Issue
Block a user