[PATCH] math: let architectures signal underflow after narrowing
Matt Turner <[email protected]>
| Newsgroups | gmane.comp.lib.glibc.alpha |
|---|---|
| Message-ID | <[email protected]> |
IEEE 754 determines tininess after rounding from the result rounded to
the precision of the destination with an unbounded exponent range. Alpha
determines it from the delivered result instead, and the two differ
where a value is tiny but reaches the smallest normal once rounded,
since the spacing below the smallest normal is twice that of the binade
the unbounded rounding lands in. Alpha signals no underflow for such a
result, and the narrowing functions cannot recover it: no trap is taken,
so the kernel emulation never sees the operation, and the final
narrowing is a single conversion instruction that has already decided
the question.
Add a CHECK_NARROW_TINY hook, in a sysdeps header that generic code
defines to do nothing, and use it in the round-to-odd narrowing macros
and in fmaf. Round-to-odd already gives the right answer wherever
tininess follows IEEE 754, so nothing changes for architectures that do.
The narrowing macros bound the round-to-odd value only as a temporary,
so add, subtract, multiply and divide now name it; fused multiply-add
already did. NARROW_MIN_NORMAL selects the smallest normal of the result
type rather than taking it as a macro argument, which would have to be
threaded through every narrowing function. On architectures where the
hook is empty the preprocessor discards both arguments, so neither is
expanded.
fmaf is not a narrowing function but converts its double result to float
on return, and so has the same problem; route its returns through one
helper.
On x86_64 the generated code is unchanged: s_fdiv and s_fmaf are
identical with and without this, in both the static and shared builds,
once debug information is stripped.
On Alpha this fixes the twenty remaining math failures, all of the form
Exception "Underflow" not set, leaving math with no failures.
---
math/math-narrow.h | 49 ++++++++++++++++-----
sysdeps/alpha/math-narrow-tininess.h | 61 ++++++++++++++++++++++++++
sysdeps/generic/math-narrow-tininess.h | 29 ++++++++++++
sysdeps/ieee754/dbl-64/s_fmaf.c | 24 ++++++++--
4 files changed, 148 insertions(+), 15 deletions(-)
create mode 100644 sysdeps/alpha/math-narrow-tininess.h
create mode 100644 sysdeps/generic/math-narrow-tininess.h
diff --git ./math/math-narrow.h ./math/math-narrow.h
index bb319cc2ca..ef79ef2380 100644
--- ./math/math-narrow.h
+++ ./math/math-narrow.h
@@ -28,8 +28,20 @@
#include <math_private.h>
#include <fenv_private.h>
#include <math-narrow-alias.h>
+#include <math-narrow-tininess.h>
+#include <float.h>
#include <stdbool.h>
+/* The smallest positive normal value of TYPE, for CHECK_NARROW_TINY.
+ Selected at compile time; the generic CHECK_NARROW_TINY ignores its
+ arguments, so this expands to nothing on architectures that determine
+ tininess as IEEE 754 describes it. */
+#define NARROW_MIN_NORMAL(TYPE) \
+ __builtin_choose_expr \
+ (__builtin_types_compatible_p (TYPE, float), FLT_MIN, \
+ __builtin_choose_expr \
+ (__builtin_types_compatible_p (TYPE, double), DBL_MIN, LDBL_MIN))
+
/* Carry out a computation using round-to-odd. The computation is
EXPR; the union type in which to store the result is UNION and the
subfield of the "ieee" field of that union with the low part of the
@@ -97,8 +109,12 @@
if ((X) == -(Y)) \
ret = (TYPE) ((X) + (Y)); \
else \
- ret = (TYPE) ROUND_TO_ODD (math_opt_barrier (X) + (Y), \
- UNION, SUFFIX, MANTISSA, false); \
+ { \
+ __typeof (X) w = ROUND_TO_ODD (math_opt_barrier (X) + (Y), \
+ UNION, SUFFIX, MANTISSA, false); \
+ ret = (TYPE) w; \
+ CHECK_NARROW_TINY (ret, w, NARROW_MIN_NORMAL (TYPE)); \
+ } \
\
CHECK_NARROW_ADD (ret, (X), (Y)); \
return ret; \
@@ -155,8 +171,12 @@
if ((X) == (Y)) \
ret = (TYPE) ((X) - (Y)); \
else \
- ret = (TYPE) ROUND_TO_ODD (math_opt_barrier (X) - (Y), \
- UNION, SUFFIX, MANTISSA, false); \
+ { \
+ __typeof (X) w = ROUND_TO_ODD (math_opt_barrier (X) - (Y), \
+ UNION, SUFFIX, MANTISSA, false); \
+ ret = (TYPE) w; \
+ CHECK_NARROW_TINY (ret, w, NARROW_MIN_NORMAL (TYPE)); \
+ } \
\
CHECK_NARROW_SUB (ret, (X), (Y)); \
return ret; \
@@ -209,9 +229,11 @@
{ \
TYPE ret; \
\
- ret = (TYPE) ROUND_TO_ODD (math_opt_barrier (X) * (Y), \
- UNION, SUFFIX, MANTISSA, \
- CLEAR_UNDERFLOW); \
+ __typeof (X) w = ROUND_TO_ODD (math_opt_barrier (X) * (Y), \
+ UNION, SUFFIX, MANTISSA, \
+ CLEAR_UNDERFLOW); \
+ ret = (TYPE) w; \
+ CHECK_NARROW_TINY (ret, w, NARROW_MIN_NORMAL (TYPE)); \
\
CHECK_NARROW_MUL (ret, (X), (Y)); \
return ret; \
@@ -264,9 +286,11 @@
{ \
TYPE ret; \
\
- ret = (TYPE) ROUND_TO_ODD (math_opt_barrier (X) / (Y), \
- UNION, SUFFIX, MANTISSA, \
- CLEAR_UNDERFLOW); \
+ __typeof (X) w = ROUND_TO_ODD (math_opt_barrier (X) / (Y), \
+ UNION, SUFFIX, MANTISSA, \
+ CLEAR_UNDERFLOW); \
+ ret = (TYPE) w; \
+ CHECK_NARROW_TINY (ret, w, NARROW_MIN_NORMAL (TYPE)); \
\
CHECK_NARROW_DIV (ret, (X), (Y)); \
return ret; \
@@ -371,7 +395,10 @@
if (tmp == 0) \
ret = (TYPE) (math_opt_barrier (X) * (Y) + (Z)); \
else \
- ret = (TYPE) tmp; \
+ { \
+ ret = (TYPE) tmp; \
+ CHECK_NARROW_TINY (ret, tmp, NARROW_MIN_NORMAL (TYPE)); \
+ } \
\
CHECK_NARROW_FMA (ret, (X), (Y), (Z)); \
return ret; \
diff --git ./sysdeps/alpha/math-narrow-tininess.h ./sysdeps/alpha/math-narrow-tininess.h
new file mode 100644
index 0000000000..66dbda9688
--- /dev/null
+++ ./sysdeps/alpha/math-narrow-tininess.h
@@ -0,0 +1,61 @@
+/* Underflow signalling for narrowing functions. Alpha version.
+ Copyright (C) 2026 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+
+ The GNU C Library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 of the License, or (at your option) any later version.
+
+ The GNU C Library is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public
+ License along with the GNU C Library; if not, see
+ <https://www.gnu.org/licenses/>. */
+
+#ifndef _MATH_NARROW_TININESS_H
+#define _MATH_NARROW_TININESS_H 1
+
+#include <fenv.h>
+
+/* IEEE 754 determines tininess after rounding from the result rounded to
+ the precision of the destination with an unbounded exponent range, but
+ Alpha determines it from the delivered result. The two differ where
+ the value is tiny but reaches the smallest normal once rounded, since
+ the spacing below the smallest normal is twice that of the binade the
+ unbounded rounding lands in. Alpha signals no underflow for such a
+ result. For example fdiv (-0x4p-128, 0x1.000002p+0) in FE_DOWNWARD has
+ the exact quotient -0x1.fffffc0000080p-127, which rounds with an
+ unbounded exponent range to -0x1.fffffep-127 and so is tiny, but is
+ delivered as -0x1p-126.
+
+ Signal underflow for that case. RET is the narrowed result, W the
+ round-to-odd value it came from and MIN_NORMAL the smallest positive
+ normal value of the narrower type. Rounding W with the exponent range
+ unbounded is what decides whether underflow is due: dividing by
+ MIN_NORMAL is exact and moves the value into the normal range, so the
+ conversion that follows rounds the significand alone.
+
+ The comparisons are the quiet ones, as RET is a NaN whenever an
+ argument was. */
+#ifdef FE_UNDERFLOW
+# define CHECK_NARROW_TINY(RET, W, MIN_NORMAL) \
+ do \
+ { \
+ if (((RET) == (MIN_NORMAL) || (RET) == -(MIN_NORMAL)) \
+ && isless (W, MIN_NORMAL) && isgreater (W, -(MIN_NORMAL))) \
+ { \
+ __typeof (RET) __scaled = (__typeof (RET)) ((W) / (MIN_NORMAL)); \
+ if (__scaled > -1 && __scaled < 1) \
+ __feraiseexcept (FE_UNDERFLOW); \
+ } \
+ } \
+ while (0)
+#else
+# define CHECK_NARROW_TINY(RET, W, MIN_NORMAL) do { } while (0)
+#endif
+
+#endif /* math-narrow-tininess.h */
diff --git ./sysdeps/generic/math-narrow-tininess.h ./sysdeps/generic/math-narrow-tininess.h
new file mode 100644
index 0000000000..9fc1e62a02
--- /dev/null
+++ ./sysdeps/generic/math-narrow-tininess.h
@@ -0,0 +1,29 @@
+/* Underflow signalling for narrowing functions. Generic version.
+ Copyright (C) 2026 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+
+ The GNU C Library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 of the License, or (at your option) any later version.
+
+ The GNU C Library is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public
+ License along with the GNU C Library; if not, see
+ <https://www.gnu.org/licenses/>. */
+
+#ifndef _MATH_NARROW_TININESS_H
+#define _MATH_NARROW_TININESS_H 1
+
+/* Narrowing a round-to-odd value signals underflow by itself wherever
+ tininess is determined as IEEE 754 describes it, so there is nothing
+ to do here. RET is the narrowed result, W the round-to-odd value it
+ came from and MIN_NORMAL the smallest positive normal value of the
+ narrower type. */
+#define CHECK_NARROW_TINY(RET, W, MIN_NORMAL) do { } while (0)
+
+#endif /* math-narrow-tininess.h */
diff --git ./sysdeps/ieee754/dbl-64/s_fmaf.c ./sysdeps/ieee754/dbl-64/s_fmaf.c
index 639c273486..f3b206d40e 100644
--- ./sysdeps/ieee754/dbl-64/s_fmaf.c
+++ ./sysdeps/ieee754/dbl-64/s_fmaf.c
@@ -21,8 +21,24 @@
#include <fenv.h>
#include <libm-alias-float.h>
#include <math-use-builtins.h>
+#include <math-narrow-tininess.h>
+#include <float.h>
#include "math_config.h"
+#if ! USE_FMAF_BUILTIN
+/* Narrow the double result to float. The conversion is what raises
+ underflow, so architectures that do not determine tininess as IEEE 754
+ describes it compensate here; elsewhere CHECK_NARROW_TINY does
+ nothing and this is a plain conversion. */
+static inline float
+narrow_fmaf_result (double result)
+{
+ float ret = (float) result;
+ CHECK_NARROW_TINY (ret, result, FLT_MIN);
+ return ret;
+}
+#endif
+
float
__fmaf (float x, float y, float z)
{
@@ -37,17 +53,17 @@ __fmaf (float x, float y, float z)
/* If not exact or at round to even boundary, the result is correct in
all rounding modes. */
if (__glibc_likely ((u & 0xfffffff) != 0))
- return result;
+ return narrow_fmaf_result (result);
/* Also check if the double result appears exact when it might not be and
thus it will not set the underflow flag if denormal. */
if ((u & 0x10000000) == 0
&& ((u >> MANTISSA_WIDTH) & 0x7ff) > EXPONENT_BIAS - 126)
- return result;
+ return narrow_fmaf_result (result);
/* Return if result is exact in all rounding modes. */
if (result - xy == z && result - z == xy)
- return result;
+ return narrow_fmaf_result (result);
/* This is where 'double-rouding' might return a wrong value, and thus
needs adjusting the low-order bits in the direction of the error. */
@@ -62,7 +78,7 @@ __fmaf (float x, float y, float z)
u++;
else
u--;
- return asdouble (u);
+ return narrow_fmaf_result (asdouble (u));
#endif /* ! USE_FMAF_BUILTIN */
}
#ifndef __fmaf
--
2.54.0