Index: lib/libm/src/fpmath.h =================================================================== RCS file: lib/libm/src/fpmath.h diff -N lib/libm/src/fpmath.h --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ lib/libm/src/fpmath.h 1 Jun 2019 15:55:40 -0000 @@ -0,0 +1,48 @@ +/*- + * SPDX-License-Identifier: BSD-2-Clause-FreeBSD + * + * Copyright (c) 2002, 2003 David Schultz + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS + * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + * + * $FreeBSD$ + */ + +#include +#include + +#define manl ext_fracl +#define manh ext_frach +#define exp ext_exp +#define sign ext_sign +#define junkl ext_padl +#define junkh ext_padh + +union IEEEl2bits { + long double e; + struct ieee_ext bits; +}; + +#define LDBL_TO_ARRAY32(p, a) EXT_TO_ARRAY32(p, a) +#define LDBL_MANH_SIZE DBL_FRACHBITS +#define LDBL_MANL_SIZE DBL_FRACLBITS Index: lib/libm/src/math_private.h =================================================================== RCS file: /cvs/src/lib/libm/src/math_private.h,v retrieving revision 1.18 diff -u -p -r1.18 math_private.h --- lib/libm/src/math_private.h 12 Sep 2016 19:47:02 -0000 1.18 +++ lib/libm/src/math_private.h 1 Jun 2019 10:30:47 -0000 @@ -236,6 +236,10 @@ typedef union u_int32_t msw; u_int32_t lsw; } parts; + struct + { + u_int64_t w; + } xparts; } ieee_double_shape_type; #endif @@ -250,6 +254,10 @@ typedef union u_int32_t lsw; u_int32_t msw; } parts; + struct + { + u_int64_t w; + } xparts; } ieee_double_shape_type; #endif @@ -264,6 +272,14 @@ do { \ (ix1) = ew_u.parts.lsw; \ } while (0) +/* Get a 64-bit int from a double. */ +#define EXTRACT_WORD64(ix,d) \ +do { \ + ieee_double_shape_type ew_u; \ + ew_u.value = (d); \ + (ix) = ew_u.xparts.w; \ +} while (0) + /* Get the more significant 32 bit int from a double. */ #define GET_HIGH_WORD(i,d) \ @@ -290,6 +306,14 @@ do { \ iw_u.parts.msw = (ix0); \ iw_u.parts.lsw = (ix1); \ (d) = iw_u.value; \ +} while (0) + +/* Set a double from a 64-bit int. */ +#define INSERT_WORD64(d,ix) \ +do { \ + ieee_double_shape_type iw_u; \ + iw_u.xparts.w = (ix); \ + (d) = iw_u.value; \ } while (0) /* Set the more significant 32 bits of a double from an int. */ Index: lib/libm/src/s_fma.c =================================================================== RCS file: /cvs/src/lib/libm/src/s_fma.c,v retrieving revision 1.7 diff -u -p -r1.7 s_fma.c --- lib/libm/src/s_fma.c 12 Sep 2016 19:47:02 -0000 1.7 +++ lib/libm/src/s_fma.c 1 Jun 2019 11:20:54 -0000 @@ -1,7 +1,7 @@ -/* $OpenBSD: s_fma.c,v 1.7 2016/09/12 19:47:02 guenther Exp $ */ - /*- - * Copyright (c) 2005 David Schultz + * SPDX-License-Identifier: BSD-2-Clause-FreeBSD + * + * Copyright (c) 2005-2011 David Schultz * All rights reserved. * * Redistribution and use in source and binary forms, with or without @@ -26,10 +26,138 @@ * SUCH DAMAGE. */ +#include + #include #include #include +#include "math_private.h" + +/* + * A struct dd represents a floating-point number with twice the precision + * of a double. We maintain the invariant that "hi" stores the 53 high-order + * bits of the result. + */ +struct dd { + double hi; + double lo; +}; + +/* + * Compute a+b exactly, returning the exact result in a struct dd. We assume + * that both a and b are finite, but make no assumptions about their relative + * magnitudes. + */ +static inline struct dd +dd_add(double a, double b) +{ + struct dd ret; + double s; + + ret.hi = a + b; + s = ret.hi - a; + ret.lo = (a - (ret.hi - s)) + (b - s); + return (ret); +} + +/* + * Compute a+b, with a small tweak: The least significant bit of the + * result is adjusted into a sticky bit summarizing all the bits that + * were lost to rounding. This adjustment negates the effects of double + * rounding when the result is added to another number with a higher + * exponent. For an explanation of round and sticky bits, see any reference + * on FPU design, e.g., + * + * J. Coonen. An Implementation Guide to a Proposed Standard for + * Floating-Point Arithmetic. Computer, vol. 13, no. 1, Jan 1980. + */ +static inline double +add_adjusted(double a, double b) +{ + struct dd sum; + uint64_t hibits, lobits; + + sum = dd_add(a, b); + if (sum.lo != 0) { + EXTRACT_WORD64(hibits, sum.hi); + if ((hibits & 1) == 0) { + /* hibits += (int)copysign(1.0, sum.hi * sum.lo) */ + EXTRACT_WORD64(lobits, sum.lo); + hibits += 1 - ((hibits ^ lobits) >> 62); + INSERT_WORD64(sum.hi, hibits); + } + } + return (sum.hi); +} + +/* + * Compute ldexp(a+b, scale) with a single rounding error. It is assumed + * that the result will be subnormal, and care is taken to ensure that + * double rounding does not occur. + */ +static inline double +add_and_denormalize(double a, double b, int scale) +{ + struct dd sum; + uint64_t hibits, lobits; + int bits_lost; + + sum = dd_add(a, b); + + /* + * If we are losing at least two bits of accuracy to denormalization, + * then the first lost bit becomes a round bit, and we adjust the + * lowest bit of sum.hi to make it a sticky bit summarizing all the + * bits in sum.lo. With the sticky bit adjusted, the hardware will + * break any ties in the correct direction. + * + * If we are losing only one bit to denormalization, however, we must + * break the ties manually. + */ + if (sum.lo != 0) { + EXTRACT_WORD64(hibits, sum.hi); + bits_lost = -((int)(hibits >> 52) & 0x7ff) - scale + 1; + if ((bits_lost != 1) ^ (int)(hibits & 1)) { + /* hibits += (int)copysign(1.0, sum.hi * sum.lo) */ + EXTRACT_WORD64(lobits, sum.lo); + hibits += 1 - (((hibits ^ lobits) >> 62) & 2); + INSERT_WORD64(sum.hi, hibits); + } + } + return (ldexp(sum.hi, scale)); +} + +/* + * Compute a*b exactly, returning the exact result in a struct dd. We assume + * that both a and b are normalized, so no underflow or overflow will occur. + * The current rounding mode must be round-to-nearest. + */ +static inline struct dd +dd_mul(double a, double b) +{ + static const double split = 0x1p27 + 1.0; + struct dd ret; + double ha, hb, la, lb, p, q; + + p = a * split; + ha = a - p; + ha += p; + la = a - ha; + + p = b * split; + hb = b - p; + hb += p; + lb = b - hb; + + p = ha * hb; + q = ha * lb + la * hb; + + ret.hi = p + q; + ret.lo = p - ret.hi + q + la * lb; + return (ret); +} + /* * Fused multiply-add: Compute x * y + z with a single rounding error. * @@ -47,14 +175,11 @@ * Hardware instructions should be used on architectures that support it, * since this implementation will likely be several times slower. */ -#if LDBL_MANT_DIG != 113 double fma(double x, double y, double z) { - static const double split = 0x1p27 + 1.0; - double xs, ys, zs; - double c, cc, hx, hy, p, q, tx, ty; - double r, rr, s; + double xs, ys, zs, adj; + struct dd xy, r; int oround; int ex, ey, ez; int spread; @@ -84,41 +209,6 @@ fma(double x, double y, double z) * will overflow, so we handle these cases specially. Rounding * modes other than FE_TONEAREST are painful. */ - if (spread > DBL_MANT_DIG * 2) { - fenv_t env; - feraiseexcept(FE_INEXACT); - switch(oround) { - case FE_TONEAREST: - return (x * y); - case FE_TOWARDZERO: - if ((x > 0.0) ^ (y < 0.0) ^ (z < 0.0)) - return (x * y); - feholdexcept(&env); - r = x * y; - if (!fetestexcept(FE_INEXACT)) - r = nextafter(r, 0); - feupdateenv(&env); - return (r); - case FE_DOWNWARD: - if (z > 0.0) - return (x * y); - feholdexcept(&env); - r = x * y; - if (!fetestexcept(FE_INEXACT)) - r = nextafter(r, -INFINITY); - feupdateenv(&env); - return (r); - default: /* FE_UPWARD */ - if (z < 0.0) - return (x * y); - feholdexcept(&env); - r = x * y; - if (!fetestexcept(FE_INEXACT)) - r = nextafter(r, INFINITY); - feupdateenv(&env); - return (r); - } - } if (spread < -DBL_MANT_DIG) { feraiseexcept(FE_INEXACT); if (!isnormal(z)) @@ -127,78 +217,75 @@ fma(double x, double y, double z) case FE_TONEAREST: return (z); case FE_TOWARDZERO: - if ((x > 0.0) ^ (y < 0.0) ^ (z < 0.0)) + if (x > 0.0 ^ y < 0.0 ^ z < 0.0) return (z); else return (nextafter(z, 0)); case FE_DOWNWARD: - if ((x > 0.0) ^ (y < 0.0)) + if (x > 0.0 ^ y < 0.0) return (z); else return (nextafter(z, -INFINITY)); default: /* FE_UPWARD */ - if ((x > 0.0) ^ (y < 0.0)) + if (x > 0.0 ^ y < 0.0) return (nextafter(z, INFINITY)); else return (z); } } + if (spread <= DBL_MANT_DIG * 2) + zs = ldexp(zs, -spread); + else + zs = copysign(DBL_MIN, zs); - /* - * Use Dekker's algorithm to perform the multiplication and - * subsequent addition in twice the machine precision. - * Arrange so that x * y = c + cc, and x * y + z = r + rr. - */ fesetround(FE_TONEAREST); + /* work around clang bug 8100 */ + volatile double vxs = xs; - p = xs * split; - hx = xs - p; - hx += p; - tx = xs - hx; - - p = ys * split; - hy = ys - p; - hy += p; - ty = ys - hy; - - p = hx * hy; - q = hx * ty + tx * hy; - c = p + q; - cc = p - c + q + tx * ty; - - zs = ldexp(zs, -spread); - r = c + zs; - s = r - c; - rr = (c - (r - s)) + (zs - s) + cc; + /* + * Basic approach for round-to-nearest: + * + * (xy.hi, xy.lo) = x * y (exact) + * (r.hi, r.lo) = xy.hi + z (exact) + * adj = xy.lo + r.lo (inexact; low bit is sticky) + * result = r.hi + adj (correctly rounded) + */ + xy = dd_mul(vxs, ys); + r = dd_add(xy.hi, zs); spread = ex + ey; - if (spread + ilogb(r) > -1023) { + + if (r.hi == 0.0) { + /* + * When the addends cancel to 0, ensure that the result has + * the correct sign. + */ fesetround(oround); - r = r + rr; - } else { + volatile double vzs = zs; /* XXX gcc CSE bug workaround */ + return (xy.hi + vzs + ldexp(xy.lo, spread)); + } + + if (oround != FE_TONEAREST) { /* - * The result is subnormal, so we round before scaling to - * avoid double rounding. + * There is no need to worry about double rounding in directed + * rounding modes. */ - p = ldexp(copysign(0x1p-1022, r), -spread); - c = r + p; - s = c - r; - cc = (r - (c - s)) + (p - s) + rr; fesetround(oround); - r = (c + cc) - p; + /* work around clang bug 8100 */ + volatile double vrlo = r.lo; + adj = vrlo + xy.lo; + return (ldexp(r.hi + adj, spread)); } - return (ldexp(r, spread)); -} -#else /* LDBL_MANT_DIG == 113 */ -/* - * 113 bits of precision is more than twice the precision of a double, - * so it is enough to represent the intermediate product exactly. - */ -double -fma(double x, double y, double z) -{ - return ((long double)x * y + z); + + adj = add_adjusted(r.lo, xy.lo); + if (spread + ilogb(r.hi) > -1023) + return (ldexp(r.hi + adj, spread)); + else + return (add_and_denormalize(r.hi, adj, spread)); } -#endif /* LDBL_MANT_DIG != 113 */ + +#if (LDBL_MANT_DIG == 53) +__weak_reference(fma, fmal); +#endif DEF_STD(fma); LDBL_MAYBE_UNUSED_CLONE(fma); Index: lib/libm/src/s_fmaf.c =================================================================== RCS file: /cvs/src/lib/libm/src/s_fmaf.c,v retrieving revision 1.2 diff -u -p -r1.2 s_fmaf.c --- lib/libm/src/s_fmaf.c 5 Dec 2012 23:20:04 -0000 1.2 +++ lib/libm/src/s_fmaf.c 1 Jun 2019 10:31:18 -0000 @@ -1,7 +1,7 @@ -/* $OpenBSD: s_fmaf.c,v 1.2 2012/12/05 23:20:04 deraadt Exp $ */ - /*- - * Copyright (c) 2005 David Schultz + * SPDX-License-Identifier: BSD-2-Clause-FreeBSD + * + * Copyright (c) 2005-2011 David Schultz * All rights reserved. * * Redistribution and use in source and binary forms, with or without @@ -26,22 +26,45 @@ * SUCH DAMAGE. */ -#include +#include + +#include + +#include "math.h" +#include "math_private.h" /* * Fused multiply-add: Compute x * y + z with a single rounding error. * * A double has more than twice as much precision than a float, so - * direct double-precision arithmetic suffices. - * - * XXX We are relying on the compiler to convert from double to float - * using the current rounding mode and with the appropriate - * side-effects. But on at least one platform (gcc 3.4.2/sparc64), - * this appears to be too much to ask for. The precision - * reduction should be done manually. + * direct double-precision arithmetic suffices, except where double + * rounding occurs. */ float fmaf(float x, float y, float z) { - return ((double)x * y + z); + double xy, result; + uint32_t hr, lr; + + xy = (double)x * y; + result = xy + z; + EXTRACT_WORDS(hr, lr, result); + /* Common case: The double precision result is fine. */ + if ((lr & 0x1fffffff) != 0x10000000 || /* not a halfway case */ + (hr & 0x7ff00000) == 0x7ff00000 || /* NaN */ + result - xy == z || /* exact */ + fegetround() != FE_TONEAREST) /* not round-to-nearest */ + return (result); + + /* + * If result is inexact, and exactly halfway between two float values, + * we need to adjust the low-order bit in the direction of the error. + */ + fesetround(FE_TOWARDZERO); + volatile double vxy = xy; /* XXX work around gcc CSE bug */ + double adjusted_result = vxy + z; + fesetround(FE_TONEAREST); + if (result == adjusted_result) + SET_LOW_WORD(adjusted_result, lr + 1); + return (adjusted_result); } Index: lib/libm/src/s_fmal.c =================================================================== RCS file: /cvs/src/lib/libm/src/s_fmal.c,v retrieving revision 1.3 diff -u -p -r1.3 s_fmal.c --- lib/libm/src/s_fmal.c 12 Nov 2013 19:00:38 -0000 1.3 +++ lib/libm/src/s_fmal.c 1 Jun 2019 10:35:01 -0000 @@ -1,7 +1,7 @@ -/* $OpenBSD: s_fmal.c,v 1.3 2013/11/12 19:00:38 martynas Exp $ */ - /*- - * Copyright (c) 2005 David Schultz + * SPDX-License-Identifier: BSD-2-Clause-FreeBSD + * + * Copyright (c) 2005-2011 David Schultz * All rights reserved. * * Redistribution and use in source and binary forms, with or without @@ -26,10 +26,134 @@ * SUCH DAMAGE. */ +#include + #include #include #include +#include "fpmath.h" + +/* + * A struct dd represents a floating-point number with twice the precision + * of a long double. We maintain the invariant that "hi" stores the high-order + * bits of the result. + */ +struct dd { + long double hi; + long double lo; +}; + +/* + * Compute a+b exactly, returning the exact result in a struct dd. We assume + * that both a and b are finite, but make no assumptions about their relative + * magnitudes. + */ +static inline struct dd +dd_add(long double a, long double b) +{ + struct dd ret; + long double s; + + ret.hi = a + b; + s = ret.hi - a; + ret.lo = (a - (ret.hi - s)) + (b - s); + return (ret); +} + +/* + * Compute a+b, with a small tweak: The least significant bit of the + * result is adjusted into a sticky bit summarizing all the bits that + * were lost to rounding. This adjustment negates the effects of double + * rounding when the result is added to another number with a higher + * exponent. For an explanation of round and sticky bits, see any reference + * on FPU design, e.g., + * + * J. Coonen. An Implementation Guide to a Proposed Standard for + * Floating-Point Arithmetic. Computer, vol. 13, no. 1, Jan 1980. + */ +static inline long double +add_adjusted(long double a, long double b) +{ + struct dd sum; + union IEEEl2bits u; + + sum = dd_add(a, b); + if (sum.lo != 0) { + u.e = sum.hi; + if ((u.bits.manl & 1) == 0) + sum.hi = nextafterl(sum.hi, INFINITY * sum.lo); + } + return (sum.hi); +} + +/* + * Compute ldexp(a+b, scale) with a single rounding error. It is assumed + * that the result will be subnormal, and care is taken to ensure that + * double rounding does not occur. + */ +static inline long double +add_and_denormalize(long double a, long double b, int scale) +{ + struct dd sum; + int bits_lost; + union IEEEl2bits u; + + sum = dd_add(a, b); + + /* + * If we are losing at least two bits of accuracy to denormalization, + * then the first lost bit becomes a round bit, and we adjust the + * lowest bit of sum.hi to make it a sticky bit summarizing all the + * bits in sum.lo. With the sticky bit adjusted, the hardware will + * break any ties in the correct direction. + * + * If we are losing only one bit to denormalization, however, we must + * break the ties manually. + */ + if (sum.lo != 0) { + u.e = sum.hi; + bits_lost = -u.bits.exp - scale + 1; + if ((bits_lost != 1) ^ (int)(u.bits.manl & 1)) + sum.hi = nextafterl(sum.hi, INFINITY * sum.lo); + } + return (ldexp(sum.hi, scale)); +} + +/* + * Compute a*b exactly, returning the exact result in a struct dd. We assume + * that both a and b are normalized, so no underflow or overflow will occur. + * The current rounding mode must be round-to-nearest. + */ +static inline struct dd +dd_mul(long double a, long double b) +{ +#if LDBL_MANT_DIG == 64 + static const long double split = 0x1p32L + 1.0; +#elif LDBL_MANT_DIG == 113 + static const long double split = 0x1p57L + 1.0; +#endif + struct dd ret; + long double ha, hb, la, lb, p, q; + + p = a * split; + ha = a - p; + ha += p; + la = a - ha; + + p = b * split; + hb = b - p; + hb += p; + lb = b - hb; + + p = ha * hb; + q = ha * lb + la * hb; + + ret.hi = p + q; + ret.lo = p - ret.hi + q + la * lb; + return (ret); +} + /* * Fused multiply-add: Compute x * y + z with a single rounding error. * @@ -42,14 +166,8 @@ long double fmal(long double x, long double y, long double z) { -#if LDBL_MANT_DIG == 64 - static const long double split = 0x1p32L + 1.0; -#elif LDBL_MANT_DIG == 113 - static const long double split = 0x1p57L + 1.0; -#endif - long double xs, ys, zs; - long double c, cc, hx, hy, p, q, tx, ty; - long double r, rr, s; + long double xs, ys, zs, adj; + struct dd xy, r; int oround; int ex, ey, ez; int spread; @@ -79,41 +197,6 @@ fmal(long double x, long double y, long * will overflow, so we handle these cases specially. Rounding * modes other than FE_TONEAREST are painful. */ - if (spread > LDBL_MANT_DIG * 2) { - fenv_t env; - feraiseexcept(FE_INEXACT); - switch(oround) { - case FE_TONEAREST: - return (x * y); - case FE_TOWARDZERO: - if ((x > 0.0) ^ (y < 0.0) ^ (z < 0.0)) - return (x * y); - feholdexcept(&env); - r = x * y; - if (!fetestexcept(FE_INEXACT)) - r = nextafterl(r, 0); - feupdateenv(&env); - return (r); - case FE_DOWNWARD: - if (z > 0.0) - return (x * y); - feholdexcept(&env); - r = x * y; - if (!fetestexcept(FE_INEXACT)) - r = nextafterl(r, -INFINITY); - feupdateenv(&env); - return (r); - default: /* FE_UPWARD */ - if (z < 0.0) - return (x * y); - feholdexcept(&env); - r = x * y; - if (!fetestexcept(FE_INEXACT)) - r = nextafterl(r, INFINITY); - feupdateenv(&env); - return (r); - } - } if (spread < -LDBL_MANT_DIG) { feraiseexcept(FE_INEXACT); if (!isnormal(z)) @@ -122,65 +205,69 @@ fmal(long double x, long double y, long case FE_TONEAREST: return (z); case FE_TOWARDZERO: - if ((x > 0.0) ^ (y < 0.0) ^ (z < 0.0)) + if (x > 0.0 ^ y < 0.0 ^ z < 0.0) return (z); else return (nextafterl(z, 0)); case FE_DOWNWARD: - if ((x > 0.0) ^ (y < 0.0)) + if (x > 0.0 ^ y < 0.0) return (z); else return (nextafterl(z, -INFINITY)); default: /* FE_UPWARD */ - if ((x > 0.0) ^ (y < 0.0)) + if (x > 0.0 ^ y < 0.0) return (nextafterl(z, INFINITY)); else return (z); } } + if (spread <= LDBL_MANT_DIG * 2) + zs = ldexpl(zs, -spread); + else + zs = copysignl(LDBL_MIN, zs); - /* - * Use Dekker's algorithm to perform the multiplication and - * subsequent addition in twice the machine precision. - * Arrange so that x * y = c + cc, and x * y + z = r + rr. - */ fesetround(FE_TONEAREST); + /* work around clang bug 8100 */ + volatile long double vxs = xs; - p = xs * split; - hx = xs - p; - hx += p; - tx = xs - hx; - - p = ys * split; - hy = ys - p; - hy += p; - ty = ys - hy; - - p = hx * hy; - q = hx * ty + tx * hy; - c = p + q; - cc = p - c + q + tx * ty; - - zs = ldexpl(zs, -spread); - r = c + zs; - s = r - c; - rr = (c - (r - s)) + (zs - s) + cc; + /* + * Basic approach for round-to-nearest: + * + * (xy.hi, xy.lo) = x * y (exact) + * (r.hi, r.lo) = xy.hi + z (exact) + * adj = xy.lo + r.lo (inexact; low bit is sticky) + * result = r.hi + adj (correctly rounded) + */ + xy = dd_mul(vxs, ys); + r = dd_add(xy.hi, zs); spread = ex + ey; - if (spread + ilogbl(r) > -16383) { + + if (r.hi == 0.0) { + /* + * When the addends cancel to 0, ensure that the result has + * the correct sign. + */ fesetround(oround); - r = r + rr; - } else { + volatile long double vzs = zs; /* XXX gcc CSE bug workaround */ + return (xy.hi + vzs + ldexpl(xy.lo, spread)); + } + + if (oround != FE_TONEAREST) { /* - * The result is subnormal, so we round before scaling to - * avoid double rounding. + * There is no need to worry about double rounding in directed + * rounding modes. */ - p = ldexpl(copysignl(0x1p-16382L, r), -spread); - c = r + p; - s = c - r; - cc = (r - (c - s)) + (p - s) + rr; fesetround(oround); - r = (c + cc) - p; + /* work around clang bug 8100 */ + volatile long double vrlo = r.lo; + adj = vrlo + xy.lo; + return (ldexpl(r.hi + adj, spread)); } - return (ldexpl(r, spread)); + + adj = add_adjusted(r.lo, xy.lo); + if (spread + ilogbl(r.hi) > -16383) + return (ldexpl(r.hi + adj, spread)); + else + return (add_and_denormalize(r.hi, adj, spread)); } Index: regress/lib/libm/msun/Makefile =================================================================== RCS file: /cvs/src/regress/lib/libm/msun/Makefile,v retrieving revision 1.1.1.1 diff -u -p -r1.1.1.1 Makefile --- regress/lib/libm/msun/Makefile 21 Feb 2019 16:14:03 -0000 1.1.1.1 +++ regress/lib/libm/msun/Makefile 1 Jun 2019 11:20:01 -0000 @@ -3,6 +3,7 @@ TESTS = TESTS += conj_test TESTS += fenv_test +TESTS += fma_test TESTS += lrint_test PROGS= ${TESTS} Index: regress/lib/libm/msun/fma_test.c =================================================================== RCS file: regress/lib/libm/msun/fma_test.c diff -N regress/lib/libm/msun/fma_test.c --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ regress/lib/libm/msun/fma_test.c 1 Jun 2019 11:55:40 -0000 @@ -0,0 +1,543 @@ +/*- + * Copyright (c) 2008 David Schultz + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS + * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + */ + +/* + * Tests for fma{,f,l}(). + */ + +#include + +#include +#include +#include +#include +#include +#include +#include + +#include "test-utils.h" + +#pragma STDC FENV_ACCESS ON + +/* + * Test that a function returns the correct value and sets the + * exception flags correctly. The exceptmask specifies which + * exceptions we should check. We need to be lenient for several + * reasons, but mainly because on some architectures it's impossible + * to raise FE_OVERFLOW without raising FE_INEXACT. + * + * These are macros instead of functions so that assert provides more + * meaningful error messages. + */ +#define test(func, x, y, z, result, exceptmask, excepts) do { \ + volatile long double _vx = (x), _vy = (y), _vz = (z); \ + assert(feclearexcept(FE_ALL_EXCEPT) == 0); \ + assert(fpequal((func)(_vx, _vy, _vz), (result))); \ + assert(((void)(func), fetestexcept(exceptmask) == (excepts))); \ +} while (0) + +#define testall(x, y, z, result, exceptmask, excepts) do { \ + test(fma, (double)(x), (double)(y), (double)(z), \ + (double)(result), (exceptmask), (excepts)); \ + test(fmaf, (float)(x), (float)(y), (float)(z), \ + (float)(result), (exceptmask), (excepts)); \ + test(fmal, (x), (y), (z), (result), (exceptmask), (excepts)); \ +} while (0) + +/* Test in all rounding modes. */ +#define testrnd(func, x, y, z, rn, ru, rd, rz, exceptmask, excepts) do { \ + fesetround(FE_TONEAREST); \ + test((func), (x), (y), (z), (rn), (exceptmask), (excepts)); \ + fesetround(FE_UPWARD); \ + test((func), (x), (y), (z), (ru), (exceptmask), (excepts)); \ + fesetround(FE_DOWNWARD); \ + test((func), (x), (y), (z), (rd), (exceptmask), (excepts)); \ + fesetround(FE_TOWARDZERO); \ + test((func), (x), (y), (z), (rz), (exceptmask), (excepts)); \ +} while (0) + +/* + * This is needed because clang constant-folds fma in ways that are incorrect + * in rounding modes other than FE_TONEAREST. + */ +static volatile double one = 1.0; + +static void +test_zeroes(void) +{ + const int rd = (fegetround() == FE_DOWNWARD); + + testall(0.0, 0.0, 0.0, 0.0, ALL_STD_EXCEPT, 0); + testall(1.0, 0.0, 0.0, 0.0, ALL_STD_EXCEPT, 0); + testall(0.0, 1.0, 0.0, 0.0, ALL_STD_EXCEPT, 0); + testall(0.0, 0.0, 1.0, 1.0, ALL_STD_EXCEPT, 0); + + testall(-0.0, 0.0, 0.0, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0); + testall(0.0, -0.0, 0.0, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0); + testall(-0.0, -0.0, 0.0, 0.0, ALL_STD_EXCEPT, 0); + testall(0.0, 0.0, -0.0, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0); + testall(-0.0, -0.0, -0.0, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0); + + testall(-0.0, 0.0, -0.0, -0.0, ALL_STD_EXCEPT, 0); + testall(0.0, -0.0, -0.0, -0.0, ALL_STD_EXCEPT, 0); + + testall(-one, one, one, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0); + testall(one, -one, one, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0); + testall(-one, -one, -one, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0); + + switch (fegetround()) { + case FE_TONEAREST: + case FE_TOWARDZERO: + test(fmaf, -FLT_MIN, FLT_MIN, 0.0, -0.0, + ALL_STD_EXCEPT, FE_INEXACT | FE_UNDERFLOW); + test(fma, -DBL_MIN, DBL_MIN, 0.0, -0.0, + ALL_STD_EXCEPT, FE_INEXACT | FE_UNDERFLOW); + test(fmal, -LDBL_MIN, LDBL_MIN, 0.0, -0.0, + ALL_STD_EXCEPT, FE_INEXACT | FE_UNDERFLOW); + } +} + +static void +test_infinities(void) +{ + + testall(INFINITY, 1.0, -1.0, INFINITY, ALL_STD_EXCEPT, 0); + testall(-1.0, INFINITY, 0.0, -INFINITY, ALL_STD_EXCEPT, 0); + testall(0.0, 0.0, INFINITY, INFINITY, ALL_STD_EXCEPT, 0); + testall(1.0, 1.0, INFINITY, INFINITY, ALL_STD_EXCEPT, 0); + testall(1.0, 1.0, -INFINITY, -INFINITY, ALL_STD_EXCEPT, 0); + + testall(INFINITY, -INFINITY, 1.0, -INFINITY, ALL_STD_EXCEPT, 0); + testall(INFINITY, INFINITY, 1.0, INFINITY, ALL_STD_EXCEPT, 0); + testall(-INFINITY, -INFINITY, INFINITY, INFINITY, ALL_STD_EXCEPT, 0); + + testall(0.0, INFINITY, 1.0, NAN, ALL_STD_EXCEPT, FE_INVALID); + testall(INFINITY, 0.0, -0.0, NAN, ALL_STD_EXCEPT, FE_INVALID); + + /* The invalid exception is optional in this case. */ + testall(INFINITY, 0.0, NAN, NAN, ALL_STD_EXCEPT & ~FE_INVALID, 0); + + testall(INFINITY, INFINITY, -INFINITY, NAN, + ALL_STD_EXCEPT, FE_INVALID); + testall(-INFINITY, INFINITY, INFINITY, NAN, + ALL_STD_EXCEPT, FE_INVALID); + testall(INFINITY, -1.0, INFINITY, NAN, + ALL_STD_EXCEPT, FE_INVALID); + + test(fmaf, FLT_MAX, FLT_MAX, -INFINITY, -INFINITY, ALL_STD_EXCEPT, 0); + test(fma, DBL_MAX, DBL_MAX, -INFINITY, -INFINITY, ALL_STD_EXCEPT, 0); + test(fmal, LDBL_MAX, LDBL_MAX, -INFINITY, -INFINITY, + ALL_STD_EXCEPT, 0); + test(fmaf, FLT_MAX, -FLT_MAX, INFINITY, INFINITY, ALL_STD_EXCEPT, 0); + test(fma, DBL_MAX, -DBL_MAX, INFINITY, INFINITY, ALL_STD_EXCEPT, 0); + test(fmal, LDBL_MAX, -LDBL_MAX, INFINITY, INFINITY, + ALL_STD_EXCEPT, 0); +} + +static void +test_nans(void) +{ + + testall(NAN, 0.0, 0.0, NAN, ALL_STD_EXCEPT, 0); + testall(1.0, NAN, 1.0, NAN, ALL_STD_EXCEPT, 0); + testall(1.0, -1.0, NAN, NAN, ALL_STD_EXCEPT, 0); + testall(0.0, 0.0, NAN, NAN, ALL_STD_EXCEPT, 0); + testall(NAN, NAN, NAN, NAN, ALL_STD_EXCEPT, 0); + + /* x*y should not raise an inexact/overflow/underflow if z is NaN. */ + testall(M_PI, M_PI, NAN, NAN, ALL_STD_EXCEPT, 0); + test(fmaf, FLT_MIN, FLT_MIN, NAN, NAN, ALL_STD_EXCEPT, 0); + test(fma, DBL_MIN, DBL_MIN, NAN, NAN, ALL_STD_EXCEPT, 0); + test(fmal, LDBL_MIN, LDBL_MIN, NAN, NAN, ALL_STD_EXCEPT, 0); + test(fmaf, FLT_MAX, FLT_MAX, NAN, NAN, ALL_STD_EXCEPT, 0); + test(fma, DBL_MAX, DBL_MAX, NAN, NAN, ALL_STD_EXCEPT, 0); + test(fmal, LDBL_MAX, LDBL_MAX, NAN, NAN, ALL_STD_EXCEPT, 0); +} + +/* + * Tests for cases where z is very small compared to x*y. + */ +static void +test_small_z(void) +{ + + /* x*y positive, z positive */ + if (fegetround() == FE_UPWARD) { + test(fmaf, one, one, 0x1.0p-100, 1.0 + FLT_EPSILON, + ALL_STD_EXCEPT, FE_INEXACT); + test(fma, one, one, 0x1.0p-200, 1.0 + DBL_EPSILON, + ALL_STD_EXCEPT, FE_INEXACT); + test(fmal, one, one, 0x1.0p-200, 1.0 + LDBL_EPSILON, + ALL_STD_EXCEPT, FE_INEXACT); + } else { + testall(0x1.0p100, one, 0x1.0p-100, 0x1.0p100, + ALL_STD_EXCEPT, FE_INEXACT); + } + + /* x*y negative, z negative */ + if (fegetround() == FE_DOWNWARD) { + test(fmaf, -one, one, -0x1.0p-100, -(1.0 + FLT_EPSILON), + ALL_STD_EXCEPT, FE_INEXACT); + test(fma, -one, one, -0x1.0p-200, -(1.0 + DBL_EPSILON), + ALL_STD_EXCEPT, FE_INEXACT); + test(fmal, -one, one, -0x1.0p-200, -(1.0 + LDBL_EPSILON), + ALL_STD_EXCEPT, FE_INEXACT); + } else { + testall(0x1.0p100, -one, -0x1.0p-100, -0x1.0p100, + ALL_STD_EXCEPT, FE_INEXACT); + } + + /* x*y positive, z negative */ + if (fegetround() == FE_DOWNWARD || fegetround() == FE_TOWARDZERO) { + test(fmaf, one, one, -0x1.0p-100, 1.0 - FLT_EPSILON / 2, + ALL_STD_EXCEPT, FE_INEXACT); + test(fma, one, one, -0x1.0p-200, 1.0 - DBL_EPSILON / 2, + ALL_STD_EXCEPT, FE_INEXACT); + test(fmal, one, one, -0x1.0p-200, 1.0 - LDBL_EPSILON / 2, + ALL_STD_EXCEPT, FE_INEXACT); + } else { + testall(0x1.0p100, one, -0x1.0p-100, 0x1.0p100, + ALL_STD_EXCEPT, FE_INEXACT); + } + + /* x*y negative, z positive */ + if (fegetround() == FE_UPWARD || fegetround() == FE_TOWARDZERO) { + test(fmaf, -one, one, 0x1.0p-100, -1.0 + FLT_EPSILON / 2, + ALL_STD_EXCEPT, FE_INEXACT); + test(fma, -one, one, 0x1.0p-200, -1.0 + DBL_EPSILON / 2, + ALL_STD_EXCEPT, FE_INEXACT); + test(fmal, -one, one, 0x1.0p-200, -1.0 + LDBL_EPSILON / 2, + ALL_STD_EXCEPT, FE_INEXACT); + } else { + testall(-0x1.0p100, one, 0x1.0p-100, -0x1.0p100, + ALL_STD_EXCEPT, FE_INEXACT); + } +} + +/* + * Tests for cases where z is very large compared to x*y. + */ +static void +test_big_z(void) +{ + + /* z positive, x*y positive */ + if (fegetround() == FE_UPWARD) { + test(fmaf, 0x1.0p-50, 0x1.0p-50, 1.0, 1.0 + FLT_EPSILON, + ALL_STD_EXCEPT, FE_INEXACT); + test(fma, 0x1.0p-100, 0x1.0p-100, 1.0, 1.0 + DBL_EPSILON, + ALL_STD_EXCEPT, FE_INEXACT); + test(fmal, 0x1.0p-100, 0x1.0p-100, 1.0, 1.0 + LDBL_EPSILON, + ALL_STD_EXCEPT, FE_INEXACT); + } else { + testall(-0x1.0p-50, -0x1.0p-50, 0x1.0p100, 0x1.0p100, + ALL_STD_EXCEPT, FE_INEXACT); + } + + /* z negative, x*y negative */ + if (fegetround() == FE_DOWNWARD) { + test(fmaf, -0x1.0p-50, 0x1.0p-50, -1.0, -(1.0 + FLT_EPSILON), + ALL_STD_EXCEPT, FE_INEXACT); + test(fma, -0x1.0p-100, 0x1.0p-100, -1.0, -(1.0 + DBL_EPSILON), + ALL_STD_EXCEPT, FE_INEXACT); + test(fmal, -0x1.0p-100, 0x1.0p-100, -1.0, -(1.0 + LDBL_EPSILON), + ALL_STD_EXCEPT, FE_INEXACT); + } else { + testall(0x1.0p-50, -0x1.0p-50, -0x1.0p100, -0x1.0p100, + ALL_STD_EXCEPT, FE_INEXACT); + } + + /* z negative, x*y positive */ + if (fegetround() == FE_UPWARD || fegetround() == FE_TOWARDZERO) { + test(fmaf, -0x1.0p-50, -0x1.0p-50, -1.0, + -1.0 + FLT_EPSILON / 2, ALL_STD_EXCEPT, FE_INEXACT); + test(fma, -0x1.0p-100, -0x1.0p-100, -1.0, + -1.0 + DBL_EPSILON / 2, ALL_STD_EXCEPT, FE_INEXACT); + test(fmal, -0x1.0p-100, -0x1.0p-100, -1.0, + -1.0 + LDBL_EPSILON / 2, ALL_STD_EXCEPT, FE_INEXACT); + } else { + testall(0x1.0p-50, 0x1.0p-50, -0x1.0p100, -0x1.0p100, + ALL_STD_EXCEPT, FE_INEXACT); + } + + /* z positive, x*y negative */ + if (fegetround() == FE_DOWNWARD || fegetround() == FE_TOWARDZERO) { + test(fmaf, 0x1.0p-50, -0x1.0p-50, 1.0, 1.0 - FLT_EPSILON / 2, + ALL_STD_EXCEPT, FE_INEXACT); + test(fma, 0x1.0p-100, -0x1.0p-100, 1.0, 1.0 - DBL_EPSILON / 2, + ALL_STD_EXCEPT, FE_INEXACT); + test(fmal, 0x1.0p-100, -0x1.0p-100, 1.0, 1.0 - LDBL_EPSILON / 2, + ALL_STD_EXCEPT, FE_INEXACT); + } else { + testall(-0x1.0p-50, 0x1.0p-50, 0x1.0p100, 0x1.0p100, + ALL_STD_EXCEPT, FE_INEXACT); + } +} + +static void +test_accuracy(void) +{ + + /* ilogb(x*y) - ilogb(z) = 20 */ + testrnd(fmaf, -0x1.c139d8p-51, -0x1.600e7ap32, 0x1.26558cp-38, + 0x1.34e48ap-18, 0x1.34e48cp-18, 0x1.34e48ap-18, 0x1.34e48ap-18, + ALL_STD_EXCEPT, FE_INEXACT); + testrnd(fma, -0x1.c139d7b84f1a3p-51, -0x1.600e7a2a16484p32, + 0x1.26558cac31580p-38, 0x1.34e48a78aae97p-18, + 0x1.34e48a78aae97p-18, 0x1.34e48a78aae96p-18, + 0x1.34e48a78aae96p-18, ALL_STD_EXCEPT, FE_INEXACT); +#if LDBL_MANT_DIG == 113 + testrnd(fmal, -0x1.c139d7b84f1a3079263afcc5bae3p-51L, + -0x1.600e7a2a164840edbe2e7d301a72p32L, + 0x1.26558cac315807eb07e448042101p-38L, + 0x1.34e48a78aae96c76ed36077dd387p-18L, + 0x1.34e48a78aae96c76ed36077dd388p-18L, + 0x1.34e48a78aae96c76ed36077dd387p-18L, + 0x1.34e48a78aae96c76ed36077dd387p-18L, + ALL_STD_EXCEPT, FE_INEXACT); +#elif LDBL_MANT_DIG == 64 + testrnd(fmal, -0x1.c139d7b84f1a307ap-51L, -0x1.600e7a2a164840eep32L, + 0x1.26558cac315807ecp-38L, 0x1.34e48a78aae96c78p-18L, + 0x1.34e48a78aae96c78p-18L, 0x1.34e48a78aae96c76p-18L, + 0x1.34e48a78aae96c76p-18L, ALL_STD_EXCEPT, FE_INEXACT); +#elif LDBL_MANT_DIG == 53 + testrnd(fmal, -0x1.c139d7b84f1a3p-51L, -0x1.600e7a2a16484p32L, + 0x1.26558cac31580p-38L, 0x1.34e48a78aae97p-18L, + 0x1.34e48a78aae97p-18L, 0x1.34e48a78aae96p-18L, + 0x1.34e48a78aae96p-18L, ALL_STD_EXCEPT, FE_INEXACT); +#endif + + /* ilogb(x*y) - ilogb(z) = -40 */ + testrnd(fmaf, 0x1.98210ap53, 0x1.9556acp-24, 0x1.d87da4p70, + 0x1.d87da4p70, 0x1.d87da6p70, 0x1.d87da4p70, 0x1.d87da4p70, + ALL_STD_EXCEPT, FE_INEXACT); + testrnd(fma, 0x1.98210ac83fe2bp53, 0x1.9556ac1475f0fp-24, + 0x1.d87da3aafc60ep70, 0x1.d87da3aafda40p70, + 0x1.d87da3aafda40p70, 0x1.d87da3aafda3fp70, + 0x1.d87da3aafda3fp70, ALL_STD_EXCEPT, FE_INEXACT); +#if LDBL_MANT_DIG == 113 + testrnd(fmal, 0x1.98210ac83fe2a8f65b6278b74cebp53L, + 0x1.9556ac1475f0f28968b61d0de65ap-24L, + 0x1.d87da3aafc60d830aa4c6d73b749p70L, + 0x1.d87da3aafda3f36a69eb86488224p70L, + 0x1.d87da3aafda3f36a69eb86488225p70L, + 0x1.d87da3aafda3f36a69eb86488224p70L, + 0x1.d87da3aafda3f36a69eb86488224p70L, + ALL_STD_EXCEPT, FE_INEXACT); +#elif LDBL_MANT_DIG == 64 + testrnd(fmal, 0x1.98210ac83fe2a8f6p53L, 0x1.9556ac1475f0f28ap-24L, + 0x1.d87da3aafc60d83p70L, 0x1.d87da3aafda3f36ap70L, + 0x1.d87da3aafda3f36ap70L, 0x1.d87da3aafda3f368p70L, + 0x1.d87da3aafda3f368p70L, ALL_STD_EXCEPT, FE_INEXACT); +#elif LDBL_MANT_DIG == 53 + testrnd(fmal, 0x1.98210ac83fe2bp53L, 0x1.9556ac1475f0fp-24L, + 0x1.d87da3aafc60ep70L, 0x1.d87da3aafda40p70L, + 0x1.d87da3aafda40p70L, 0x1.d87da3aafda3fp70L, + 0x1.d87da3aafda3fp70L, ALL_STD_EXCEPT, FE_INEXACT); +#endif + + /* ilogb(x*y) - ilogb(z) = 0 */ + testrnd(fmaf, 0x1.31ad02p+100, 0x1.2fbf7ap-42, -0x1.c3e106p+58, + -0x1.64c27cp+56, -0x1.64c27ap+56, -0x1.64c27cp+56, + -0x1.64c27ap+56, ALL_STD_EXCEPT, FE_INEXACT); + testrnd(fma, 0x1.31ad012ede8aap+100, 0x1.2fbf79c839067p-42, + -0x1.c3e106929056ep+58, -0x1.64c282b970a5fp+56, + -0x1.64c282b970a5ep+56, -0x1.64c282b970a5fp+56, + -0x1.64c282b970a5ep+56, ALL_STD_EXCEPT, FE_INEXACT); +#if LDBL_MANT_DIG == 113 + testrnd(fmal, 0x1.31ad012ede8aa282fa1c19376d16p+100L, + 0x1.2fbf79c839066f0f5c68f6d2e814p-42L, + -0x1.c3e106929056ec19de72bfe64215p+58L, + -0x1.64c282b970a612598fc025ca8cddp+56L, + -0x1.64c282b970a612598fc025ca8cddp+56L, + -0x1.64c282b970a612598fc025ca8cdep+56L, + -0x1.64c282b970a612598fc025ca8cddp+56L, + ALL_STD_EXCEPT, FE_INEXACT); +#elif LDBL_MANT_DIG == 64 + testrnd(fmal, 0x1.31ad012ede8aa4eap+100L, 0x1.2fbf79c839066aeap-42L, + -0x1.c3e106929056e61p+58L, -0x1.64c282b970a60298p+56L, + -0x1.64c282b970a60298p+56L, -0x1.64c282b970a6029ap+56L, + -0x1.64c282b970a60298p+56L, ALL_STD_EXCEPT, FE_INEXACT); +#elif LDBL_MANT_DIG == 53 + testrnd(fmal, 0x1.31ad012ede8aap+100L, 0x1.2fbf79c839067p-42L, + -0x1.c3e106929056ep+58L, -0x1.64c282b970a5fp+56L, + -0x1.64c282b970a5ep+56L, -0x1.64c282b970a5fp+56L, + -0x1.64c282b970a5ep+56L, ALL_STD_EXCEPT, FE_INEXACT); +#endif + + /* x*y (rounded) ~= -z */ + /* XXX spurious inexact exceptions */ + testrnd(fmaf, 0x1.bbffeep-30, -0x1.1d164cp-74, 0x1.ee7296p-104, + -0x1.c46ea8p-128, -0x1.c46ea8p-128, -0x1.c46ea8p-128, + -0x1.c46ea8p-128, ALL_STD_EXCEPT & ~FE_INEXACT, 0); + testrnd(fma, 0x1.bbffeea6fc7d6p-30, 0x1.1d164c6cbf078p-74, + -0x1.ee72993aff948p-104, -0x1.71f72ac7d9d8p-159, + -0x1.71f72ac7d9d8p-159, -0x1.71f72ac7d9d8p-159, + -0x1.71f72ac7d9d8p-159, ALL_STD_EXCEPT & ~FE_INEXACT, 0); +#if LDBL_MANT_DIG == 113 + testrnd(fmal, 0x1.bbffeea6fc7d65927d147f437675p-30L, + 0x1.1d164c6cbf078b7a22607d1cd6a2p-74L, + -0x1.ee72993aff94973876031bec0944p-104L, + 0x1.64e086175b3a2adc36e607058814p-217L, + 0x1.64e086175b3a2adc36e607058814p-217L, + 0x1.64e086175b3a2adc36e607058814p-217L, + 0x1.64e086175b3a2adc36e607058814p-217L, + ALL_STD_EXCEPT & ~FE_INEXACT, 0); +#elif LDBL_MANT_DIG == 64 + testrnd(fmal, 0x1.bbffeea6fc7d6592p-30L, 0x1.1d164c6cbf078b7ap-74L, + -0x1.ee72993aff949736p-104L, 0x1.af190e7a1ee6ad94p-168L, + 0x1.af190e7a1ee6ad94p-168L, 0x1.af190e7a1ee6ad94p-168L, + 0x1.af190e7a1ee6ad94p-168L, ALL_STD_EXCEPT & ~FE_INEXACT, 0); +#elif LDBL_MANT_DIG == 53 + testrnd(fmal, 0x1.bbffeea6fc7d6p-30L, 0x1.1d164c6cbf078p-74L, + -0x1.ee72993aff948p-104L, -0x1.71f72ac7d9d8p-159L, + -0x1.71f72ac7d9d8p-159L, -0x1.71f72ac7d9d8p-159L, + -0x1.71f72ac7d9d8p-159L, ALL_STD_EXCEPT & ~FE_INEXACT, 0); +#endif +} + +static void +test_double_rounding(void) +{ + + /* + * a = 0x1.8000000000001p0 + * b = 0x1.8000000000001p0 + * c = -0x0.0000000000000000000000000080...1p+1 + * a * b = 0x1.2000000000001800000000000080p+1 + * + * The correct behavior is to round DOWN to 0x1.2000000000001p+1 in + * round-to-nearest mode. An implementation that computes a*b+c in + * double+double precision, however, will get 0x1.20000000000018p+1, + * and then round UP. + */ + fesetround(FE_TONEAREST); + test(fma, 0x1.8000000000001p0, 0x1.8000000000001p0, + -0x1.0000000000001p-104, 0x1.2000000000001p+1, + ALL_STD_EXCEPT, FE_INEXACT); + fesetround(FE_DOWNWARD); + test(fma, 0x1.8000000000001p0, 0x1.8000000000001p0, + -0x1.0000000000001p-104, 0x1.2000000000001p+1, + ALL_STD_EXCEPT, FE_INEXACT); + fesetround(FE_UPWARD); + test(fma, 0x1.8000000000001p0, 0x1.8000000000001p0, + -0x1.0000000000001p-104, 0x1.2000000000002p+1, + ALL_STD_EXCEPT, FE_INEXACT); + + fesetround(FE_TONEAREST); + test(fmaf, 0x1.800002p+0, 0x1.800002p+0, -0x1.000002p-46, 0x1.200002p+1, + ALL_STD_EXCEPT, FE_INEXACT); + fesetround(FE_DOWNWARD); + test(fmaf, 0x1.800002p+0, 0x1.800002p+0, -0x1.000002p-46, 0x1.200002p+1, + ALL_STD_EXCEPT, FE_INEXACT); + fesetround(FE_UPWARD); + test(fmaf, 0x1.800002p+0, 0x1.800002p+0, -0x1.000002p-46, 0x1.200004p+1, + ALL_STD_EXCEPT, FE_INEXACT); + + fesetround(FE_TONEAREST); +#if LDBL_MANT_DIG == 64 + test(fmal, 0x1.4p+0L, 0x1.0000000000000004p+0L, 0x1p-128L, + 0x1.4000000000000006p+0L, ALL_STD_EXCEPT, FE_INEXACT); +#elif LDBL_MANT_DIG == 113 + test(fmal, 0x1.8000000000000000000000000001p+0L, + 0x1.8000000000000000000000000001p+0L, + -0x1.0000000000000000000000000001p-224L, + 0x1.2000000000000000000000000001p+1L, ALL_STD_EXCEPT, FE_INEXACT); +#endif + +} + +int +main(void) +{ + int rmodes[] = { FE_TONEAREST, FE_UPWARD, FE_DOWNWARD, FE_TOWARDZERO }; + unsigned i, j; + +#if defined(__i386__) + printf("1..0 # SKIP all testcases fail on i386\n"); + exit(0); +#endif + + j = 1; + + printf("1..19\n"); + + for (i = 0; i < nitems(rmodes); i++, j++) { + printf("rmode = %d\n", rmodes[i]); + fesetround(rmodes[i]); + test_zeroes(); + printf("ok %d - fma zeroes\n", j); + } + + for (i = 0; i < nitems(rmodes); i++, j++) { +#if defined(__amd64__) + printf("ok %d # SKIP testcase fails assertion on " + "amd64\n", j); + continue; +#else + printf("rmode = %d\n", rmodes[i]); + fesetround(rmodes[i]); + test_infinities(); + printf("ok %d - fma infinities\n", j); +#endif + } + + fesetround(FE_TONEAREST); + test_nans(); + printf("ok %d - fma NaNs\n", j); + j++; + + for (i = 0; i < nitems(rmodes); i++, j++) { + printf("rmode = %d\n", rmodes[i]); + fesetround(rmodes[i]); + test_small_z(); + printf("ok %d - fma small z\n", j); + } + + for (i = 0; i < nitems(rmodes); i++, j++) { + printf("rmode = %d\n", rmodes[i]); + fesetround(rmodes[i]); + test_big_z(); + printf("ok %d - fma big z\n", j); + } + + fesetround(FE_TONEAREST); + test_accuracy(); + printf("ok %d - fma accuracy\n", j); + j++; + + test_double_rounding(); + printf("ok %d - fma double rounding\n", j); + j++; + + /* + * TODO: + * - Tests for subnormals + * - Cancellation tests (e.g., z = (double)x*y, but x*y is inexact) + */ + + return (0); +}