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 #include #include +#include +#include #include +/* 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 + . */ + +#ifndef _MATH_NARROW_TININESS_H +#define _MATH_NARROW_TININESS_H 1 + +#include + +/* 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 + . */ + +#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 #include #include +#include +#include #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