Project

General

Profile

Download (42.9 KB)

Task #6746 » fma.diff

Moritz Buhl, 06/02/2019 11:51 AM

View differences:

lib/libm/src/fpmath.h 1 Jun 2019 15:55:40 -0000
/*-
* SPDX-License-Identifier: BSD-2-Clause-FreeBSD
*
* Copyright (c) 2002, 2003 David Schultz <das@FreeBSD.ORG>
* 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 <sys/types.h>
#include <machine/ieee.h>
#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
lib/libm/src/math_private.h 1 Jun 2019 10:30:47 -0000
u_int32_t msw;
u_int32_t lsw;
} parts;
struct
{
u_int64_t w;
} xparts;
} ieee_double_shape_type;
#endif
......
u_int32_t lsw;
u_int32_t msw;
} parts;
struct
{
u_int64_t w;
} xparts;
} ieee_double_shape_type;
#endif
......
(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) \
......
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. */
lib/libm/src/s_fma.c 1 Jun 2019 11:20:54 -0000
/* $OpenBSD: s_fma.c,v 1.7 2016/09/12 19:47:02 guenther Exp $ */
/*-
* Copyright (c) 2005 David Schultz <das@FreeBSD.ORG>
* SPDX-License-Identifier: BSD-2-Clause-FreeBSD
*
* Copyright (c) 2005-2011 David Schultz <das@FreeBSD.ORG>
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
......
* SUCH DAMAGE.
*/
#include <sys/cdefs.h>
#include <fenv.h>
#include <float.h>
#include <math.h>
#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.
*
......
* 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;
......
* 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))
......
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);
lib/libm/src/s_fmaf.c 1 Jun 2019 10:31:18 -0000
/* $OpenBSD: s_fmaf.c,v 1.2 2012/12/05 23:20:04 deraadt Exp $ */
/*-
* Copyright (c) 2005 David Schultz <das@FreeBSD.ORG>
* SPDX-License-Identifier: BSD-2-Clause-FreeBSD
*
* Copyright (c) 2005-2011 David Schultz <das@FreeBSD.ORG>
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
......
* SUCH DAMAGE.
*/
#include <math.h>
#include <sys/cdefs.h>
#include <fenv.h>
#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);
}
lib/libm/src/s_fmal.c 1 Jun 2019 10:35:01 -0000
/* $OpenBSD: s_fmal.c,v 1.3 2013/11/12 19:00:38 martynas Exp $ */
/*-
* Copyright (c) 2005 David Schultz <das@FreeBSD.ORG>
* SPDX-License-Identifier: BSD-2-Clause-FreeBSD
*
* Copyright (c) 2005-2011 David Schultz <das@FreeBSD.ORG>
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
......
* SUCH DAMAGE.
*/
#include <sys/cdefs.h>
#include <fenv.h>
#include <float.h>
#include <math.h>
#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.
*
......
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;
......
* 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))
......
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));
}
regress/lib/libm/msun/Makefile 1 Jun 2019 11:20:01 -0000
TESTS =
TESTS += conj_test
TESTS += fenv_test
TESTS += fma_test
TESTS += lrint_test
PROGS= ${TESTS}
regress/lib/libm/msun/fma_test.c 1 Jun 2019 11:55:40 -0000
/*-
* Copyright (c) 2008 David Schultz <das@FreeBSD.org>
* 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 <sys/cdefs.h>
#include <sys/param.h>
#include <assert.h>
#include <fenv.h>
#include <float.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#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);
}
(2-2/4)