Move math_check_force_underflow macros to separate math-underflow.h.

This patch continues cleaning up math_private.h by moving the
math_check_force_underflow set of macros to a separate header
math-underflow.h.

This header is included by the files that need it rather than from
math_private.h.  Moving these macros to a separate file removes the
math_private.h uses of macros from float.h, so the inclusion of
float.h in math_private.h is also removed; files that were depending
on that inclusion are fixed to include float.h directly.  The
inclusion of math-barriers.h from math_private.h will be removed in a
separate patch.

Tested for x86_64 and x86.  Also tested with build-many-glibcs.py that
installed stripped shared libraries are unchanged by this patch.

	* math/math-underflow.h: New file.
	* sysdeps/generic/math_private.h: Do not include <float.h>.
	(fabs_tg): Remove macro.  Moved to math-underflow.h.
	(min_of_type_f): Likewise.
	(min_of_type_): Likewise.
	(min_of_type_l): Likewise.
	(min_of_type_f128): Likewise.
	(min_of_type): Likewise.
	(math_check_force_underflow): Likewise.
	(math_check_force_underflow_nonneg): Likewise.
	(math_check_force_underflow_complex): Likewise.
	* math/e_exp2_template.c: Include <math-underflow.h>.
	* math/k_casinh_template.c: Likewise.
	* math/s_catan_template.c: Likewise.
	* math/s_catanh_template.c: Likewise.
	* math/s_ccosh_template.c: Likewise.
	* math/s_cexp_template.c: Likewise.
	* math/s_clog10_template.c: Likewise.
	* math/s_clog_template.c: Likewise.
	* math/s_csin_template.c: Likewise.
	* math/s_csinh_template.c: Likewise.
	* math/s_csqrt_template.c: Likewise.
	* math/s_ctan_template.c: Likewise.
	* math/s_ctanh_template.c: Likewise.
	* sysdeps/ieee754/dbl-64/e_asin.c: Likewise.
	* sysdeps/ieee754/dbl-64/e_atanh.c: Likewise.
	* sysdeps/ieee754/dbl-64/e_exp2.c: Likewise.
	* sysdeps/ieee754/dbl-64/e_gamma_r.c: Likewise.
	* sysdeps/ieee754/dbl-64/e_hypot.c: Likewise.
	* sysdeps/ieee754/dbl-64/e_j1.c: Likewise.
	* sysdeps/ieee754/dbl-64/e_jn.c: Likewise.
	* sysdeps/ieee754/dbl-64/e_pow.c: Likewise.
	* sysdeps/ieee754/dbl-64/e_sinh.c: Likewise.
	* sysdeps/ieee754/dbl-64/s_asinh.c: Likewise.
	* sysdeps/ieee754/dbl-64/s_atan.c: Likewise.
	* sysdeps/ieee754/dbl-64/s_erf.c: Likewise.
	* sysdeps/ieee754/dbl-64/s_expm1.c: Likewise.
	* sysdeps/ieee754/dbl-64/s_log1p.c: Likewise.
	* sysdeps/ieee754/dbl-64/s_sin.c: Likewise.
	* sysdeps/ieee754/dbl-64/s_sincos.c: Likewise.
	* sysdeps/ieee754/dbl-64/s_tan.c: Likewise.
	* sysdeps/ieee754/dbl-64/s_tanh.c: Likewise.
	* sysdeps/ieee754/flt-32/e_asinf.c: Likewise.
	* sysdeps/ieee754/flt-32/e_atanhf.c: Likewise.
	* sysdeps/ieee754/flt-32/e_gammaf_r.c: Likewise.
	* sysdeps/ieee754/flt-32/e_j1f.c: Likewise.
	* sysdeps/ieee754/flt-32/e_jnf.c: Likewise.
	* sysdeps/ieee754/flt-32/e_sinhf.c: Likewise.
	* sysdeps/ieee754/flt-32/k_sinf.c: Likewise.
	* sysdeps/ieee754/flt-32/k_tanf.c: Likewise.
	* sysdeps/ieee754/flt-32/s_asinhf.c: Likewise.
	* sysdeps/ieee754/flt-32/s_atanf.c: Likewise.
	* sysdeps/ieee754/flt-32/s_erff.c: Likewise.
	* sysdeps/ieee754/flt-32/s_expm1f.c: Likewise.
	* sysdeps/ieee754/flt-32/s_log1pf.c: Likewise.
	* sysdeps/ieee754/flt-32/s_tanhf.c: Likewise.
	* sysdeps/ieee754/ldbl-128/e_asinl.c: Likewise.
	* sysdeps/ieee754/ldbl-128/e_atanhl.c: Likewise.
	* sysdeps/ieee754/ldbl-128/e_expl.c: Likewise.
	* sysdeps/ieee754/ldbl-128/e_gammal_r.c: Likewise.
	* sysdeps/ieee754/ldbl-128/e_hypotl.c: Likewise.
	* sysdeps/ieee754/ldbl-128/e_j1l.c: Likewise.
	* sysdeps/ieee754/ldbl-128/e_jnl.c: Likewise.
	* sysdeps/ieee754/ldbl-128/e_sinhl.c: Likewise.
	* sysdeps/ieee754/ldbl-128/k_sincosl.c: Likewise.
	* sysdeps/ieee754/ldbl-128/k_sinl.c: Likewise.
	* sysdeps/ieee754/ldbl-128/k_tanl.c: Likewise.
	* sysdeps/ieee754/ldbl-128/s_asinhl.c: Likewise.
	* sysdeps/ieee754/ldbl-128/s_atanl.c: Likewise.
	* sysdeps/ieee754/ldbl-128/s_erfl.c: Likewise.
	* sysdeps/ieee754/ldbl-128/s_expm1l.c: Likewise.
	* sysdeps/ieee754/ldbl-128/s_log1pl.c: Likewise.
	* sysdeps/ieee754/ldbl-128/s_tanhl.c: Likewise.
	* sysdeps/ieee754/ldbl-128ibm/e_asinl.c: Likewise.
	* sysdeps/ieee754/ldbl-128ibm/e_atanhl.c: Likewise.
	* sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c: Likewise.
	* sysdeps/ieee754/ldbl-128ibm/e_hypotl.c: Likewise.
	* sysdeps/ieee754/ldbl-128ibm/e_j1l.c: Likewise.
	* sysdeps/ieee754/ldbl-128ibm/e_jnl.c: Likewise.
	* sysdeps/ieee754/ldbl-128ibm/e_powl.c: Likewise.
	* sysdeps/ieee754/ldbl-128ibm/e_sinhl.c: Likewise.
	* sysdeps/ieee754/ldbl-128ibm/k_sincosl.c: Likewise.
	* sysdeps/ieee754/ldbl-128ibm/k_sinl.c: Likewise.
	* sysdeps/ieee754/ldbl-128ibm/k_tanl.c: Likewise.
	* sysdeps/ieee754/ldbl-128ibm/s_asinhl.c: Likewise.
	* sysdeps/ieee754/ldbl-128ibm/s_atanl.c: Likewise.
	* sysdeps/ieee754/ldbl-128ibm/s_erfl.c: Likewise.
	* sysdeps/ieee754/ldbl-128ibm/s_fmal.c: Likewise.
	* sysdeps/ieee754/ldbl-128ibm/s_tanhl.c: Likewise.
	* sysdeps/ieee754/ldbl-96/e_asinl.c: Likewise.
	* sysdeps/ieee754/ldbl-96/e_atanhl.c: Likewise.
	* sysdeps/ieee754/ldbl-96/e_gammal_r.c: Likewise.
	* sysdeps/ieee754/ldbl-96/e_hypotl.c: Likewise.
	* sysdeps/ieee754/ldbl-96/e_j1l.c: Likewise.
	* sysdeps/ieee754/ldbl-96/e_jnl.c: Likewise.
	* sysdeps/ieee754/ldbl-96/e_sinhl.c: Likewise.
	* sysdeps/ieee754/ldbl-96/k_sinl.c: Likewise.
	* sysdeps/ieee754/ldbl-96/k_tanl.c: Likewise.
	* sysdeps/ieee754/ldbl-96/s_asinhl.c: Likewise.
	* sysdeps/ieee754/ldbl-96/s_erfl.c: Likewise.
	* sysdeps/ieee754/ldbl-96/s_tanhl.c: Likewise.
	* sysdeps/powerpc/fpu/e_hypot.c: Likewise.
	* sysdeps/x86/fpu/powl_helper.c: Likewise.
	* sysdeps/ieee754/dbl-64/s_nextup.c: Include <float.h>.
	* sysdeps/ieee754/flt-32/s_nextupf.c: Likewise.
	* sysdeps/ieee754/ldbl-128/s_nextupl.c: Likewise.
	* sysdeps/ieee754/ldbl-128ibm/s_nextupl.c: Likewise.
	* sysdeps/ieee754/ldbl-96/s_nextupl.c: Likewise.

100 files changed, 287 insertions, 53 deletions

diff --git a/ChangeLog b/ChangeLog
index 7ab225e320..b2b66020e8 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,114 @@
+2018-05-10  Joseph Myers  <joseph@codesourcery.com>
+
+	* math/math-underflow.h: New file.
+	* sysdeps/generic/math_private.h: Do not include <float.h>.
+	(fabs_tg): Remove macro.  Moved to math-underflow.h.
+	(min_of_type_f): Likewise.
+	(min_of_type_): Likewise.
+	(min_of_type_l): Likewise.
+	(min_of_type_f128): Likewise.
+	(min_of_type): Likewise.
+	(math_check_force_underflow): Likewise.
+	(math_check_force_underflow_nonneg): Likewise.
+	(math_check_force_underflow_complex): Likewise.
+	* math/e_exp2_template.c: Include <math-underflow.h>.
+	* math/k_casinh_template.c: Likewise.
+	* math/s_catan_template.c: Likewise.
+	* math/s_catanh_template.c: Likewise.
+	* math/s_ccosh_template.c: Likewise.
+	* math/s_cexp_template.c: Likewise.
+	* math/s_clog10_template.c: Likewise.
+	* math/s_clog_template.c: Likewise.
+	* math/s_csin_template.c: Likewise.
+	* math/s_csinh_template.c: Likewise.
+	* math/s_csqrt_template.c: Likewise.
+	* math/s_ctan_template.c: Likewise.
+	* math/s_ctanh_template.c: Likewise.
+	* sysdeps/ieee754/dbl-64/e_asin.c: Likewise.
+	* sysdeps/ieee754/dbl-64/e_atanh.c: Likewise.
+	* sysdeps/ieee754/dbl-64/e_exp2.c: Likewise.
+	* sysdeps/ieee754/dbl-64/e_gamma_r.c: Likewise.
+	* sysdeps/ieee754/dbl-64/e_hypot.c: Likewise.
+	* sysdeps/ieee754/dbl-64/e_j1.c: Likewise.
+	* sysdeps/ieee754/dbl-64/e_jn.c: Likewise.
+	* sysdeps/ieee754/dbl-64/e_pow.c: Likewise.
+	* sysdeps/ieee754/dbl-64/e_sinh.c: Likewise.
+	* sysdeps/ieee754/dbl-64/s_asinh.c: Likewise.
+	* sysdeps/ieee754/dbl-64/s_atan.c: Likewise.
+	* sysdeps/ieee754/dbl-64/s_erf.c: Likewise.
+	* sysdeps/ieee754/dbl-64/s_expm1.c: Likewise.
+	* sysdeps/ieee754/dbl-64/s_log1p.c: Likewise.
+	* sysdeps/ieee754/dbl-64/s_sin.c: Likewise.
+	* sysdeps/ieee754/dbl-64/s_sincos.c: Likewise.
+	* sysdeps/ieee754/dbl-64/s_tan.c: Likewise.
+	* sysdeps/ieee754/dbl-64/s_tanh.c: Likewise.
+	* sysdeps/ieee754/flt-32/e_asinf.c: Likewise.
+	* sysdeps/ieee754/flt-32/e_atanhf.c: Likewise.
+	* sysdeps/ieee754/flt-32/e_gammaf_r.c: Likewise.
+	* sysdeps/ieee754/flt-32/e_j1f.c: Likewise.
+	* sysdeps/ieee754/flt-32/e_jnf.c: Likewise.
+	* sysdeps/ieee754/flt-32/e_sinhf.c: Likewise.
+	* sysdeps/ieee754/flt-32/k_sinf.c: Likewise.
+	* sysdeps/ieee754/flt-32/k_tanf.c: Likewise.
+	* sysdeps/ieee754/flt-32/s_asinhf.c: Likewise.
+	* sysdeps/ieee754/flt-32/s_atanf.c: Likewise.
+	* sysdeps/ieee754/flt-32/s_erff.c: Likewise.
+	* sysdeps/ieee754/flt-32/s_expm1f.c: Likewise.
+	* sysdeps/ieee754/flt-32/s_log1pf.c: Likewise.
+	* sysdeps/ieee754/flt-32/s_tanhf.c: Likewise.
+	* sysdeps/ieee754/ldbl-128/e_asinl.c: Likewise.
+	* sysdeps/ieee754/ldbl-128/e_atanhl.c: Likewise.
+	* sysdeps/ieee754/ldbl-128/e_expl.c: Likewise.
+	* sysdeps/ieee754/ldbl-128/e_gammal_r.c: Likewise.
+	* sysdeps/ieee754/ldbl-128/e_hypotl.c: Likewise.
+	* sysdeps/ieee754/ldbl-128/e_j1l.c: Likewise.
+	* sysdeps/ieee754/ldbl-128/e_jnl.c: Likewise.
+	* sysdeps/ieee754/ldbl-128/e_sinhl.c: Likewise.
+	* sysdeps/ieee754/ldbl-128/k_sincosl.c: Likewise.
+	* sysdeps/ieee754/ldbl-128/k_sinl.c: Likewise.
+	* sysdeps/ieee754/ldbl-128/k_tanl.c: Likewise.
+	* sysdeps/ieee754/ldbl-128/s_asinhl.c: Likewise.
+	* sysdeps/ieee754/ldbl-128/s_atanl.c: Likewise.
+	* sysdeps/ieee754/ldbl-128/s_erfl.c: Likewise.
+	* sysdeps/ieee754/ldbl-128/s_expm1l.c: Likewise.
+	* sysdeps/ieee754/ldbl-128/s_log1pl.c: Likewise.
+	* sysdeps/ieee754/ldbl-128/s_tanhl.c: Likewise.
+	* sysdeps/ieee754/ldbl-128ibm/e_asinl.c: Likewise.
+	* sysdeps/ieee754/ldbl-128ibm/e_atanhl.c: Likewise.
+	* sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c: Likewise.
+	* sysdeps/ieee754/ldbl-128ibm/e_hypotl.c: Likewise.
+	* sysdeps/ieee754/ldbl-128ibm/e_j1l.c: Likewise.
+	* sysdeps/ieee754/ldbl-128ibm/e_jnl.c: Likewise.
+	* sysdeps/ieee754/ldbl-128ibm/e_powl.c: Likewise.
+	* sysdeps/ieee754/ldbl-128ibm/e_sinhl.c: Likewise.
+	* sysdeps/ieee754/ldbl-128ibm/k_sincosl.c: Likewise.
+	* sysdeps/ieee754/ldbl-128ibm/k_sinl.c: Likewise.
+	* sysdeps/ieee754/ldbl-128ibm/k_tanl.c: Likewise.
+	* sysdeps/ieee754/ldbl-128ibm/s_asinhl.c: Likewise.
+	* sysdeps/ieee754/ldbl-128ibm/s_atanl.c: Likewise.
+	* sysdeps/ieee754/ldbl-128ibm/s_erfl.c: Likewise.
+	* sysdeps/ieee754/ldbl-128ibm/s_fmal.c: Likewise.
+	* sysdeps/ieee754/ldbl-128ibm/s_tanhl.c: Likewise.
+	* sysdeps/ieee754/ldbl-96/e_asinl.c: Likewise.
+	* sysdeps/ieee754/ldbl-96/e_atanhl.c: Likewise.
+	* sysdeps/ieee754/ldbl-96/e_gammal_r.c: Likewise.
+	* sysdeps/ieee754/ldbl-96/e_hypotl.c: Likewise.
+	* sysdeps/ieee754/ldbl-96/e_j1l.c: Likewise.
+	* sysdeps/ieee754/ldbl-96/e_jnl.c: Likewise.
+	* sysdeps/ieee754/ldbl-96/e_sinhl.c: Likewise.
+	* sysdeps/ieee754/ldbl-96/k_sinl.c: Likewise.
+	* sysdeps/ieee754/ldbl-96/k_tanl.c: Likewise.
+	* sysdeps/ieee754/ldbl-96/s_asinhl.c: Likewise.
+	* sysdeps/ieee754/ldbl-96/s_erfl.c: Likewise.
+	* sysdeps/ieee754/ldbl-96/s_tanhl.c: Likewise.
+	* sysdeps/powerpc/fpu/e_hypot.c: Likewise.
+	* sysdeps/x86/fpu/powl_helper.c: Likewise.
+	* sysdeps/ieee754/dbl-64/s_nextup.c: Include <float.h>.
+	* sysdeps/ieee754/flt-32/s_nextupf.c: Likewise.
+	* sysdeps/ieee754/ldbl-128/s_nextupl.c: Likewise.
+	* sysdeps/ieee754/ldbl-128ibm/s_nextupl.c: Likewise.
+	* sysdeps/ieee754/ldbl-96/s_nextupl.c: Likewise.
+
 2018-05-09  Joseph Myers  <joseph@codesourcery.com>
 
 	* sysdeps/generic/math-barriers.h: New file.
diff --git a/math/e_exp2_template.c b/math/e_exp2_template.c
index 94a21054c2..30ed4af49b 100644
--- a/math/e_exp2_template.c
+++ b/math/e_exp2_template.c
@@ -18,6 +18,7 @@
 
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <float.h>
 
 #define declare_mgen_finite_alias_x(from, to) \
diff --git a/math/k_casinh_template.c b/math/k_casinh_template.c
index f7b679e55f..451f05c279 100644
--- a/math/k_casinh_template.c
+++ b/math/k_casinh_template.c
@@ -21,6 +21,7 @@
 #include <complex.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <float.h>
 
 /* Return the complex inverse hyperbolic sine of finite nonzero Z,
diff --git a/math/math-underflow.h b/math/math-underflow.h
new file mode 100644
index 0000000000..c5e5c79768
--- /dev/null
+++ b/math/math-underflow.h
@@ -0,0 +1,79 @@
+/* Check for underflow and force underflow exceptions.
+   Copyright (C) 2015-2018 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
+   <http://www.gnu.org/licenses/>.  */
+
+#ifndef _MATH_UNDERFLOW_H
+#define _MATH_UNDERFLOW_H	1
+
+#include <float.h>
+#include <math.h>
+
+#include <math-barriers.h>
+
+#define fabs_tg(x) __MATH_TG ((x), (__typeof (x)) __builtin_fabs, (x))
+
+/* These must be function-like macros because some __MATH_TG
+   implementations macro-expand the function-name argument before
+   concatenating a suffix to it.  */
+#define min_of_type_f() FLT_MIN
+#define min_of_type_() DBL_MIN
+#define min_of_type_l() LDBL_MIN
+#define min_of_type_f128() FLT128_MIN
+
+#define min_of_type(x) __MATH_TG ((x), (__typeof (x)) min_of_type_, ())
+
+/* If X (which is not a NaN) is subnormal, force an underflow
+   exception.  */
+#define math_check_force_underflow(x)				\
+  do								\
+    {								\
+      __typeof (x) force_underflow_tmp = (x);			\
+      if (fabs_tg (force_underflow_tmp)				\
+	  < min_of_type (force_underflow_tmp))			\
+	{							\
+	  __typeof (force_underflow_tmp) force_underflow_tmp2	\
+	    = force_underflow_tmp * force_underflow_tmp;	\
+	  math_force_eval (force_underflow_tmp2);		\
+	}							\
+    }								\
+  while (0)
+/* Likewise, but X is also known to be nonnegative.  */
+#define math_check_force_underflow_nonneg(x)			\
+  do								\
+    {								\
+      __typeof (x) force_underflow_tmp = (x);			\
+      if (force_underflow_tmp					\
+	  < min_of_type (force_underflow_tmp))			\
+	{							\
+	  __typeof (force_underflow_tmp) force_underflow_tmp2	\
+	    = force_underflow_tmp * force_underflow_tmp;	\
+	  math_force_eval (force_underflow_tmp2);		\
+	}							\
+    }								\
+  while (0)
+/* Likewise, for both real and imaginary parts of a complex
+   result.  */
+#define math_check_force_underflow_complex(x)				\
+  do									\
+    {									\
+      __typeof (x) force_underflow_complex_tmp = (x);			\
+      math_check_force_underflow (__real__ force_underflow_complex_tmp); \
+      math_check_force_underflow (__imag__ force_underflow_complex_tmp); \
+    }									\
+  while (0)
+
+#endif /* math-underflow.h */
diff --git a/math/s_catan_template.c b/math/s_catan_template.c
index c3894c5a0a..19f4e11cb4 100644
--- a/math/s_catan_template.c
+++ b/math/s_catan_template.c
@@ -20,6 +20,7 @@
 #include <complex.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <float.h>
 
 CFLOAT
diff --git a/math/s_catanh_template.c b/math/s_catanh_template.c
index fdeeb589e1..245a1a6412 100644
--- a/math/s_catanh_template.c
+++ b/math/s_catanh_template.c
@@ -20,6 +20,7 @@
 #include <complex.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <float.h>
 
 CFLOAT
diff --git a/math/s_ccosh_template.c b/math/s_ccosh_template.c
index 6840c4819e..1664f0690a 100644
--- a/math/s_ccosh_template.c
+++ b/math/s_ccosh_template.c
@@ -21,6 +21,7 @@
 #include <fenv.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <float.h>
 
 CFLOAT
diff --git a/math/s_cexp_template.c b/math/s_cexp_template.c
index e099e42a4c..5fdab1e5ac 100644
--- a/math/s_cexp_template.c
+++ b/math/s_cexp_template.c
@@ -21,6 +21,7 @@
 #include <fenv.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <float.h>
 
 CFLOAT
diff --git a/math/s_clog10_template.c b/math/s_clog10_template.c
index 82a52faebf..a947502e94 100644
--- a/math/s_clog10_template.c
+++ b/math/s_clog10_template.c
@@ -20,6 +20,7 @@
 #include <complex.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <float.h>
 
 /* log_10 (2).  */
diff --git a/math/s_clog_template.c b/math/s_clog_template.c
index 8aa1f74c8b..1248cac593 100644
--- a/math/s_clog_template.c
+++ b/math/s_clog_template.c
@@ -20,6 +20,7 @@
 #include <complex.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <float.h>
 
 CFLOAT
diff --git a/math/s_csin_template.c b/math/s_csin_template.c
index d6af614177..5f95e92ee0 100644
--- a/math/s_csin_template.c
+++ b/math/s_csin_template.c
@@ -21,6 +21,7 @@
 #include <fenv.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <float.h>
 
 CFLOAT
diff --git a/math/s_csinh_template.c b/math/s_csinh_template.c
index 6466ae80c9..95f28c2b66 100644
--- a/math/s_csinh_template.c
+++ b/math/s_csinh_template.c
@@ -21,6 +21,7 @@
 #include <fenv.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <float.h>
 
 CFLOAT
diff --git a/math/s_csqrt_template.c b/math/s_csqrt_template.c
index dcf14c7537..b1153906b9 100644
--- a/math/s_csqrt_template.c
+++ b/math/s_csqrt_template.c
@@ -21,6 +21,7 @@
 #include <complex.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <float.h>
 
 CFLOAT
diff --git a/math/s_ctan_template.c b/math/s_ctan_template.c
index af898d8f81..ab11c25299 100644
--- a/math/s_ctan_template.c
+++ b/math/s_ctan_template.c
@@ -21,6 +21,7 @@
 #include <fenv.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <float.h>
 
 CFLOAT
diff --git a/math/s_ctanh_template.c b/math/s_ctanh_template.c
index 9527ab4ccc..bd1292239f 100644
--- a/math/s_ctanh_template.c
+++ b/math/s_ctanh_template.c
@@ -21,6 +21,7 @@
 #include <fenv.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <float.h>
 
 CFLOAT
diff --git a/sysdeps/generic/math_private.h b/sysdeps/generic/math_private.h
index f6faf71895..9cd0941431 100644
--- a/sysdeps/generic/math_private.h
+++ b/sysdeps/generic/math_private.h
@@ -20,7 +20,6 @@
 #include <stdint.h>
 #include <sys/types.h>
 #include <fenv.h>
-#include <float.h>
 #include <get-rounding-mode.h>
 
 /* Gather machine dependent _Floatn support.  */
@@ -265,58 +264,6 @@ extern void __docos (double __x, double __dx, double __v[]);
 
 #include <math-barriers.h>
 
-#define fabs_tg(x) __MATH_TG ((x), (__typeof (x)) __builtin_fabs, (x))
-
-/* These must be function-like macros because some __MATH_TG
-   implementations macro-expand the function-name argument before
-   concatenating a suffix to it.  */
-#define min_of_type_f() FLT_MIN
-#define min_of_type_() DBL_MIN
-#define min_of_type_l() LDBL_MIN
-#define min_of_type_f128() FLT128_MIN
-
-#define min_of_type(x) __MATH_TG ((x), (__typeof (x)) min_of_type_, ())
-
-/* If X (which is not a NaN) is subnormal, force an underflow
-   exception.  */
-#define math_check_force_underflow(x)				\
-  do								\
-    {								\
-      __typeof (x) force_underflow_tmp = (x);			\
-      if (fabs_tg (force_underflow_tmp)				\
-	  < min_of_type (force_underflow_tmp))			\
-	{							\
-	  __typeof (force_underflow_tmp) force_underflow_tmp2	\
-	    = force_underflow_tmp * force_underflow_tmp;	\
-	  math_force_eval (force_underflow_tmp2);		\
-	}							\
-    }								\
-  while (0)
-/* Likewise, but X is also known to be nonnegative.  */
-#define math_check_force_underflow_nonneg(x)			\
-  do								\
-    {								\
-      __typeof (x) force_underflow_tmp = (x);			\
-      if (force_underflow_tmp					\
-	  < min_of_type (force_underflow_tmp))			\
-	{							\
-	  __typeof (force_underflow_tmp) force_underflow_tmp2	\
-	    = force_underflow_tmp * force_underflow_tmp;	\
-	  math_force_eval (force_underflow_tmp2);		\
-	}							\
-    }								\
-  while (0)
-/* Likewise, for both real and imaginary parts of a complex
-   result.  */
-#define math_check_force_underflow_complex(x)				\
-  do									\
-    {									\
-      __typeof (x) force_underflow_complex_tmp = (x);			\
-      math_check_force_underflow (__real__ force_underflow_complex_tmp); \
-      math_check_force_underflow (__imag__ force_underflow_complex_tmp); \
-    }									\
-  while (0)
-
 /* The standards only specify one variant of the fenv.h interfaces.
    But at least for some architectures we can be more efficient if we
    know what operations are going to be performed.  Therefore we
diff --git a/sysdeps/ieee754/dbl-64/e_asin.c b/sysdeps/ieee754/dbl-64/e_asin.c
index 768cbb9ec1..6bf56945a6 100644
--- a/sysdeps/ieee754/dbl-64/e_asin.c
+++ b/sysdeps/ieee754/dbl-64/e_asin.c
@@ -42,6 +42,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 #ifndef SECTION
 # define SECTION
diff --git a/sysdeps/ieee754/dbl-64/e_atanh.c b/sysdeps/ieee754/dbl-64/e_atanh.c
index 608afb5a8c..df9c8786b2 100644
--- a/sysdeps/ieee754/dbl-64/e_atanh.c
+++ b/sysdeps/ieee754/dbl-64/e_atanh.c
@@ -39,6 +39,7 @@
 #include <inttypes.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 static const double huge = 1e300;
 
diff --git a/sysdeps/ieee754/dbl-64/e_exp2.c b/sysdeps/ieee754/dbl-64/e_exp2.c
index a944abeee8..045cbbb88d 100644
--- a/sysdeps/ieee754/dbl-64/e_exp2.c
+++ b/sysdeps/ieee754/dbl-64/e_exp2.c
@@ -30,6 +30,7 @@
 #include <fenv.h>
 #include <inttypes.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 #include "t_exp2.h"
 
diff --git a/sysdeps/ieee754/dbl-64/e_gamma_r.c b/sysdeps/ieee754/dbl-64/e_gamma_r.c
index a3750f9f31..2744549cbd 100644
--- a/sysdeps/ieee754/dbl-64/e_gamma_r.c
+++ b/sysdeps/ieee754/dbl-64/e_gamma_r.c
@@ -20,6 +20,7 @@
 #include <math.h>
 #include <math-narrow-eval.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <float.h>
 
 /* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's
diff --git a/sysdeps/ieee754/dbl-64/e_hypot.c b/sysdeps/ieee754/dbl-64/e_hypot.c
index a89075169b..a2c33cc4ed 100644
--- a/sysdeps/ieee754/dbl-64/e_hypot.c
+++ b/sysdeps/ieee754/dbl-64/e_hypot.c
@@ -44,6 +44,7 @@
 
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 double
 __ieee754_hypot (double x, double y)
diff --git a/sysdeps/ieee754/dbl-64/e_j1.c b/sysdeps/ieee754/dbl-64/e_j1.c
index d0e387e1a2..734f3ca64a 100644
--- a/sysdeps/ieee754/dbl-64/e_j1.c
+++ b/sysdeps/ieee754/dbl-64/e_j1.c
@@ -63,6 +63,7 @@
 #include <math.h>
 #include <math-narrow-eval.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 static double pone (double), qone (double);
 
diff --git a/sysdeps/ieee754/dbl-64/e_jn.c b/sysdeps/ieee754/dbl-64/e_jn.c
index 6ef0fbea1c..9181b22bb8 100644
--- a/sysdeps/ieee754/dbl-64/e_jn.c
+++ b/sysdeps/ieee754/dbl-64/e_jn.c
@@ -41,6 +41,7 @@
 #include <math.h>
 #include <math-narrow-eval.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 static const double
   invsqrtpi = 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
diff --git a/sysdeps/ieee754/dbl-64/e_pow.c b/sysdeps/ieee754/dbl-64/e_pow.c
index 542d03a7e3..96d5b23ccc 100644
--- a/sysdeps/ieee754/dbl-64/e_pow.c
+++ b/sysdeps/ieee754/dbl-64/e_pow.c
@@ -38,6 +38,7 @@
 #include "MathLib.h"
 #include "upow.tbl"
 #include <math_private.h>
+#include <math-underflow.h>
 #include <fenv.h>
 
 #ifndef SECTION
diff --git a/sysdeps/ieee754/dbl-64/e_sinh.c b/sysdeps/ieee754/dbl-64/e_sinh.c
index 49c24fb489..c4e34211ac 100644
--- a/sysdeps/ieee754/dbl-64/e_sinh.c
+++ b/sysdeps/ieee754/dbl-64/e_sinh.c
@@ -36,6 +36,7 @@ static char rcsid[] = "$NetBSD: e_sinh.c,v 1.7 1995/05/10 20:46:13 jtc Exp $";
 #include <math.h>
 #include <math-narrow-eval.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 static const double one = 1.0, shuge = 1.0e307;
 
diff --git a/sysdeps/ieee754/dbl-64/s_asinh.c b/sysdeps/ieee754/dbl-64/s_asinh.c
index 5577ccebd9..192ff8594d 100644
--- a/sysdeps/ieee754/dbl-64/s_asinh.c
+++ b/sysdeps/ieee754/dbl-64/s_asinh.c
@@ -24,6 +24,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <libm-alias-double.h>
 
 static const double
diff --git a/sysdeps/ieee754/dbl-64/s_atan.c b/sysdeps/ieee754/dbl-64/s_atan.c
index ff5718c82e..38db092d04 100644
--- a/sysdeps/ieee754/dbl-64/s_atan.c
+++ b/sysdeps/ieee754/dbl-64/s_atan.c
@@ -46,6 +46,7 @@
 #include <libm-alias-double.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <stap-probe.h>
 
 void __mpatan (mp_no *, mp_no *, int);	/* see definition in mpatan.c */
diff --git a/sysdeps/ieee754/dbl-64/s_erf.c b/sysdeps/ieee754/dbl-64/s_erf.c
index 48dfca3317..5f820604f8 100644
--- a/sysdeps/ieee754/dbl-64/s_erf.c
+++ b/sysdeps/ieee754/dbl-64/s_erf.c
@@ -117,6 +117,7 @@ static char rcsid[] = "$NetBSD: s_erf.c,v 1.8 1995/05/10 20:47:05 jtc Exp $";
 #include <math.h>
 #include <math-narrow-eval.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <libm-alias-double.h>
 #include <fix-int-fp-convert-zero.h>
 
diff --git a/sysdeps/ieee754/dbl-64/s_expm1.c b/sysdeps/ieee754/dbl-64/s_expm1.c
index 3e136e7ee0..8ef74f2196 100644
--- a/sysdeps/ieee754/dbl-64/s_expm1.c
+++ b/sysdeps/ieee754/dbl-64/s_expm1.c
@@ -112,6 +112,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <libm-alias-double.h>
 #define one Q[0]
 static const double
diff --git a/sysdeps/ieee754/dbl-64/s_log1p.c b/sysdeps/ieee754/dbl-64/s_log1p.c
index b7cf5ce321..c01f1be93b 100644
--- a/sysdeps/ieee754/dbl-64/s_log1p.c
+++ b/sysdeps/ieee754/dbl-64/s_log1p.c
@@ -81,6 +81,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <libc-diag.h>
 
 static const double
diff --git a/sysdeps/ieee754/dbl-64/s_nextup.c b/sysdeps/ieee754/dbl-64/s_nextup.c
index db13a5793b..d37a2b6657 100644
--- a/sysdeps/ieee754/dbl-64/s_nextup.c
+++ b/sysdeps/ieee754/dbl-64/s_nextup.c
@@ -16,6 +16,7 @@
    License along with the GNU C Library; if not, see
    <http://www.gnu.org/licenses/>.  */
 
+#include <float.h>
 #include <math.h>
 #include <math_private.h>
 #include <libm-alias-double.h>
diff --git a/sysdeps/ieee754/dbl-64/s_sin.c b/sysdeps/ieee754/dbl-64/s_sin.c
index ba1dbe27b6..b369ac9f5b 100644
--- a/sysdeps/ieee754/dbl-64/s_sin.c
+++ b/sysdeps/ieee754/dbl-64/s_sin.c
@@ -41,6 +41,7 @@
 #include "MathLib.h"
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <libm-alias-double.h>
 #include <fenv.h>
 
diff --git a/sysdeps/ieee754/dbl-64/s_sincos.c b/sysdeps/ieee754/dbl-64/s_sincos.c
index c7460371e4..1d8d44befe 100644
--- a/sysdeps/ieee754/dbl-64/s_sincos.c
+++ b/sysdeps/ieee754/dbl-64/s_sincos.c
@@ -21,6 +21,7 @@
 #include <math.h>
 
 #include <math_private.h>
+#include <math-underflow.h>
 #include <libm-alias-double.h>
 
 #define IN_SINCOS
diff --git a/sysdeps/ieee754/dbl-64/s_tan.c b/sysdeps/ieee754/dbl-64/s_tan.c
index 05a25eb060..04ff8b6717 100644
--- a/sysdeps/ieee754/dbl-64/s_tan.c
+++ b/sysdeps/ieee754/dbl-64/s_tan.c
@@ -41,6 +41,7 @@
 #include "MathLib.h"
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <libm-alias-double.h>
 #include <fenv.h>
 #include <stap-probe.h>
diff --git a/sysdeps/ieee754/dbl-64/s_tanh.c b/sysdeps/ieee754/dbl-64/s_tanh.c
index 321bf44142..673a97102d 100644
--- a/sysdeps/ieee754/dbl-64/s_tanh.c
+++ b/sysdeps/ieee754/dbl-64/s_tanh.c
@@ -41,6 +41,7 @@ static char rcsid[] = "$NetBSD: s_tanh.c,v 1.7 1995/05/10 20:48:22 jtc Exp $";
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <libm-alias-double.h>
 
 static const double one = 1.0, two = 2.0, tiny = 1.0e-300;
diff --git a/sysdeps/ieee754/flt-32/e_asinf.c b/sysdeps/ieee754/flt-32/e_asinf.c
index 55eb1449c5..e03073b8bd 100644
--- a/sysdeps/ieee754/flt-32/e_asinf.c
+++ b/sysdeps/ieee754/flt-32/e_asinf.c
@@ -42,6 +42,7 @@ static char rcsid[] = "$NetBSD: e_asinf.c,v 1.5 1995/05/12 04:57:25 jtc Exp $";
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 static const float
 one =  1.0000000000e+00, /* 0x3F800000 */
diff --git a/sysdeps/ieee754/flt-32/e_atanhf.c b/sysdeps/ieee754/flt-32/e_atanhf.c
index 598000a8cd..43f11ad1f0 100644
--- a/sysdeps/ieee754/flt-32/e_atanhf.c
+++ b/sysdeps/ieee754/flt-32/e_atanhf.c
@@ -39,6 +39,7 @@
 #include <inttypes.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 static const float huge = 1e30;
 
diff --git a/sysdeps/ieee754/flt-32/e_gammaf_r.c b/sysdeps/ieee754/flt-32/e_gammaf_r.c
index d640510448..8b23add347 100644
--- a/sysdeps/ieee754/flt-32/e_gammaf_r.c
+++ b/sysdeps/ieee754/flt-32/e_gammaf_r.c
@@ -20,6 +20,7 @@
 #include <math.h>
 #include <math-narrow-eval.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <float.h>
 
 /* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's
diff --git a/sysdeps/ieee754/flt-32/e_j1f.c b/sysdeps/ieee754/flt-32/e_j1f.c
index 8aadc4c31a..887b5f7d1e 100644
--- a/sysdeps/ieee754/flt-32/e_j1f.c
+++ b/sysdeps/ieee754/flt-32/e_j1f.c
@@ -18,6 +18,7 @@
 #include <math.h>
 #include <math-narrow-eval.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 static float ponef(float), qonef(float);
 
diff --git a/sysdeps/ieee754/flt-32/e_jnf.c b/sysdeps/ieee754/flt-32/e_jnf.c
index a4414ce27e..cd15ed7d4b 100644
--- a/sysdeps/ieee754/flt-32/e_jnf.c
+++ b/sysdeps/ieee754/flt-32/e_jnf.c
@@ -18,6 +18,7 @@
 #include <math.h>
 #include <math-narrow-eval.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 static const float
 two   =  2.0000000000e+00, /* 0x40000000 */
diff --git a/sysdeps/ieee754/flt-32/e_sinhf.c b/sysdeps/ieee754/flt-32/e_sinhf.c
index 17b3663dde..20f7db81ea 100644
--- a/sysdeps/ieee754/flt-32/e_sinhf.c
+++ b/sysdeps/ieee754/flt-32/e_sinhf.c
@@ -17,6 +17,7 @@
 #include <math.h>
 #include <math-narrow-eval.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 static const float one = 1.0, shuge = 1.0e37;
 
diff --git a/sysdeps/ieee754/flt-32/k_sinf.c b/sysdeps/ieee754/flt-32/k_sinf.c
index a195d59466..dcf3c35358 100644
--- a/sysdeps/ieee754/flt-32/k_sinf.c
+++ b/sysdeps/ieee754/flt-32/k_sinf.c
@@ -20,6 +20,7 @@ static char rcsid[] = "$NetBSD: k_sinf.c,v 1.4 1995/05/10 20:46:33 jtc Exp $";
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 static const float
 half =  5.0000000000e-01,/* 0x3f000000 */
diff --git a/sysdeps/ieee754/flt-32/k_tanf.c b/sysdeps/ieee754/flt-32/k_tanf.c
index 9f0e55860f..228ece2075 100644
--- a/sysdeps/ieee754/flt-32/k_tanf.c
+++ b/sysdeps/ieee754/flt-32/k_tanf.c
@@ -20,6 +20,7 @@ static char rcsid[] = "$NetBSD: k_tanf.c,v 1.4 1995/05/10 20:46:39 jtc Exp $";
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 static const float
 one   =  1.0000000000e+00, /* 0x3f800000 */
 pio4  =  7.8539812565e-01, /* 0x3f490fda */
diff --git a/sysdeps/ieee754/flt-32/s_asinhf.c b/sysdeps/ieee754/flt-32/s_asinhf.c
index 27ecdd56fa..0812b54dca 100644
--- a/sysdeps/ieee754/flt-32/s_asinhf.c
+++ b/sysdeps/ieee754/flt-32/s_asinhf.c
@@ -16,6 +16,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <libm-alias-float.h>
 
 static const float
diff --git a/sysdeps/ieee754/flt-32/s_atanf.c b/sysdeps/ieee754/flt-32/s_atanf.c
index 03a4cfd6a5..f61b0f5ef6 100644
--- a/sysdeps/ieee754/flt-32/s_atanf.c
+++ b/sysdeps/ieee754/flt-32/s_atanf.c
@@ -20,6 +20,7 @@ static char rcsid[] = "$NetBSD: s_atanf.c,v 1.4 1995/05/10 20:46:47 jtc Exp $";
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <libm-alias-float.h>
 
 static const float atanhi[] = {
diff --git a/sysdeps/ieee754/flt-32/s_erff.c b/sysdeps/ieee754/flt-32/s_erff.c
index 1d00abdbf1..693934d59a 100644
--- a/sysdeps/ieee754/flt-32/s_erff.c
+++ b/sysdeps/ieee754/flt-32/s_erff.c
@@ -22,6 +22,7 @@ static char rcsid[] = "$NetBSD: s_erff.c,v 1.4 1995/05/10 20:47:07 jtc Exp $";
 #include <math.h>
 #include <math-narrow-eval.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <libm-alias-float.h>
 #include <fix-int-fp-convert-zero.h>
 
diff --git a/sysdeps/ieee754/flt-32/s_expm1f.c b/sysdeps/ieee754/flt-32/s_expm1f.c
index 0358970778..3b04ce7751 100644
--- a/sysdeps/ieee754/flt-32/s_expm1f.c
+++ b/sysdeps/ieee754/flt-32/s_expm1f.c
@@ -17,6 +17,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <libm-alias-float.h>
 
 static const float huge = 1.0e+30;
diff --git a/sysdeps/ieee754/flt-32/s_log1pf.c b/sysdeps/ieee754/flt-32/s_log1pf.c
index 009084ee19..9d2ef89975 100644
--- a/sysdeps/ieee754/flt-32/s_log1pf.c
+++ b/sysdeps/ieee754/flt-32/s_log1pf.c
@@ -16,6 +16,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <libc-diag.h>
 
 static const float
diff --git a/sysdeps/ieee754/flt-32/s_nextupf.c b/sysdeps/ieee754/flt-32/s_nextupf.c
index 7ba8c0b79d..87ec7ba21f 100644
--- a/sysdeps/ieee754/flt-32/s_nextupf.c
+++ b/sysdeps/ieee754/flt-32/s_nextupf.c
@@ -16,6 +16,7 @@
    License along with the GNU C Library; if not, see
    <http://www.gnu.org/licenses/>.  */
 
+#include <float.h>
 #include <math.h>
 #include <math_private.h>
 #include <libm-alias-float.h>
diff --git a/sysdeps/ieee754/flt-32/s_tanhf.c b/sysdeps/ieee754/flt-32/s_tanhf.c
index c2b96cfa2b..cc3d63d35d 100644
--- a/sysdeps/ieee754/flt-32/s_tanhf.c
+++ b/sysdeps/ieee754/flt-32/s_tanhf.c
@@ -20,6 +20,7 @@ static char rcsid[] = "$NetBSD: s_tanhf.c,v 1.4 1995/05/10 20:48:24 jtc Exp $";
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <libm-alias-float.h>
 
 static const float one=1.0, two=2.0, tiny = 1.0e-30;
diff --git a/sysdeps/ieee754/ldbl-128/e_asinl.c b/sysdeps/ieee754/ldbl-128/e_asinl.c
index 53c2f1fada..6021facdfb 100644
--- a/sysdeps/ieee754/ldbl-128/e_asinl.c
+++ b/sysdeps/ieee754/ldbl-128/e_asinl.c
@@ -62,6 +62,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 static const _Float128
   one = 1,
diff --git a/sysdeps/ieee754/ldbl-128/e_atanhl.c b/sysdeps/ieee754/ldbl-128/e_atanhl.c
index 4c8c2e22e6..13b7683e17 100644
--- a/sysdeps/ieee754/ldbl-128/e_atanhl.c
+++ b/sysdeps/ieee754/ldbl-128/e_atanhl.c
@@ -35,6 +35,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 static const _Float128 one = 1, huge = L(1e4900);
 
diff --git a/sysdeps/ieee754/ldbl-128/e_expl.c b/sysdeps/ieee754/ldbl-128/e_expl.c
index c4c61cbed4..735dc9b8d0 100644
--- a/sysdeps/ieee754/ldbl-128/e_expl.c
+++ b/sysdeps/ieee754/ldbl-128/e_expl.c
@@ -65,6 +65,7 @@
 #include <fenv.h>
 #include <inttypes.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <stdlib.h>
 #include "t_expl.h"
 
diff --git a/sysdeps/ieee754/ldbl-128/e_gammal_r.c b/sysdeps/ieee754/ldbl-128/e_gammal_r.c
index 95b9fe0677..e2730b8fc2 100644
--- a/sysdeps/ieee754/ldbl-128/e_gammal_r.c
+++ b/sysdeps/ieee754/ldbl-128/e_gammal_r.c
@@ -20,6 +20,7 @@
 
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <float.h>
 
 /* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's
diff --git a/sysdeps/ieee754/ldbl-128/e_hypotl.c b/sysdeps/ieee754/ldbl-128/e_hypotl.c
index dd82b3aec8..7bafd4ae29 100644
--- a/sysdeps/ieee754/ldbl-128/e_hypotl.c
+++ b/sysdeps/ieee754/ldbl-128/e_hypotl.c
@@ -47,6 +47,7 @@
 
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 _Float128
 __ieee754_hypotl(_Float128 x, _Float128 y)
diff --git a/sysdeps/ieee754/ldbl-128/e_j1l.c b/sysdeps/ieee754/ldbl-128/e_j1l.c
index bff3a5e88c..e6f46f52e9 100644
--- a/sysdeps/ieee754/ldbl-128/e_j1l.c
+++ b/sysdeps/ieee754/ldbl-128/e_j1l.c
@@ -98,6 +98,7 @@
 #include <errno.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <float.h>
 
 /* 1 / sqrt(pi) */
diff --git a/sysdeps/ieee754/ldbl-128/e_jnl.c b/sysdeps/ieee754/ldbl-128/e_jnl.c
index 635b4e7e42..7739eec291 100644
--- a/sysdeps/ieee754/ldbl-128/e_jnl.c
+++ b/sysdeps/ieee754/ldbl-128/e_jnl.c
@@ -60,6 +60,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 static const _Float128
   invsqrtpi = L(5.6418958354775628694807945156077258584405E-1),
diff --git a/sysdeps/ieee754/ldbl-128/e_sinhl.c b/sysdeps/ieee754/ldbl-128/e_sinhl.c
index cce5b6e4fd..39e7cf3084 100644
--- a/sysdeps/ieee754/ldbl-128/e_sinhl.c
+++ b/sysdeps/ieee754/ldbl-128/e_sinhl.c
@@ -56,6 +56,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 static const _Float128 one = 1.0, shuge = L(1.0e4931),
 ovf_thresh = L(1.1357216553474703894801348310092223067821E4);
diff --git a/sysdeps/ieee754/ldbl-128/k_sincosl.c b/sysdeps/ieee754/ldbl-128/k_sincosl.c
index 36668dd0a3..6eb33ae375 100644
--- a/sysdeps/ieee754/ldbl-128/k_sincosl.c
+++ b/sysdeps/ieee754/ldbl-128/k_sincosl.c
@@ -20,6 +20,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 static const _Float128 c[] = {
 #define ONE c[0]
diff --git a/sysdeps/ieee754/ldbl-128/k_sinl.c b/sysdeps/ieee754/ldbl-128/k_sinl.c
index a9f38c7f0a..9c19fb1bf2 100644
--- a/sysdeps/ieee754/ldbl-128/k_sinl.c
+++ b/sysdeps/ieee754/ldbl-128/k_sinl.c
@@ -20,6 +20,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 static const _Float128 c[] = {
 #define ONE c[0]
diff --git a/sysdeps/ieee754/ldbl-128/k_tanl.c b/sysdeps/ieee754/ldbl-128/k_tanl.c
index e79023c69a..8da794d80c 100644
--- a/sysdeps/ieee754/ldbl-128/k_tanl.c
+++ b/sysdeps/ieee754/ldbl-128/k_tanl.c
@@ -59,6 +59,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <libc-diag.h>
 
 static const _Float128
diff --git a/sysdeps/ieee754/ldbl-128/s_asinhl.c b/sysdeps/ieee754/ldbl-128/s_asinhl.c
index 92a482658d..a733ee7695 100644
--- a/sysdeps/ieee754/ldbl-128/s_asinhl.c
+++ b/sysdeps/ieee754/ldbl-128/s_asinhl.c
@@ -32,6 +32,7 @@ static char rcsid[] = "$NetBSD: $";
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <libm-alias-ldouble.h>
 
 static const _Float128
diff --git a/sysdeps/ieee754/ldbl-128/s_atanl.c b/sysdeps/ieee754/ldbl-128/s_atanl.c
index e05368df53..022ccf4c45 100644
--- a/sysdeps/ieee754/ldbl-128/s_atanl.c
+++ b/sysdeps/ieee754/ldbl-128/s_atanl.c
@@ -62,6 +62,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <libm-alias-ldouble.h>
 
 /* arctan(k/8), k = 0, ..., 82 */
diff --git a/sysdeps/ieee754/ldbl-128/s_erfl.c b/sysdeps/ieee754/ldbl-128/s_erfl.c
index 88e91c702d..c0abfe2926 100644
--- a/sysdeps/ieee754/ldbl-128/s_erfl.c
+++ b/sysdeps/ieee754/ldbl-128/s_erfl.c
@@ -100,6 +100,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <libm-alias-ldouble.h>
 
 /* Evaluate P[n] x^n  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
diff --git a/sysdeps/ieee754/ldbl-128/s_expm1l.c b/sysdeps/ieee754/ldbl-128/s_expm1l.c
index ea0d29c493..66881af01a 100644
--- a/sysdeps/ieee754/ldbl-128/s_expm1l.c
+++ b/sysdeps/ieee754/ldbl-128/s_expm1l.c
@@ -57,6 +57,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <libm-alias-ldouble.h>
 
 /* exp(x) - 1 = x + 0.5 x^2 + x^3 P(x)/Q(x)
diff --git a/sysdeps/ieee754/ldbl-128/s_log1pl.c b/sysdeps/ieee754/ldbl-128/s_log1pl.c
index b8b2ffeba1..f3181dbc8e 100644
--- a/sysdeps/ieee754/ldbl-128/s_log1pl.c
+++ b/sysdeps/ieee754/ldbl-128/s_log1pl.c
@@ -56,6 +56,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 /* Coefficients for log(1+x) = x - x^2 / 2 + x^3 P(x)/Q(x)
  * 1/sqrt(2) <= 1+x < sqrt(2)
diff --git a/sysdeps/ieee754/ldbl-128/s_nextupl.c b/sysdeps/ieee754/ldbl-128/s_nextupl.c
index 5628c17680..20827e6d09 100644
--- a/sysdeps/ieee754/ldbl-128/s_nextupl.c
+++ b/sysdeps/ieee754/ldbl-128/s_nextupl.c
@@ -16,6 +16,7 @@
    License along with the GNU C Library; if not, see
    <http://www.gnu.org/licenses/>.  */
 
+#include <float.h>
 #include <math.h>
 #include <math_private.h>
 #include <libm-alias-ldouble.h>
diff --git a/sysdeps/ieee754/ldbl-128/s_tanhl.c b/sysdeps/ieee754/ldbl-128/s_tanhl.c
index 8f62f96c50..fc309dab0a 100644
--- a/sysdeps/ieee754/ldbl-128/s_tanhl.c
+++ b/sysdeps/ieee754/ldbl-128/s_tanhl.c
@@ -44,6 +44,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <libm-alias-ldouble.h>
 
 static const _Float128 one = 1.0, two = 2.0, tiny = L(1.0e-4900);
diff --git a/sysdeps/ieee754/ldbl-128ibm/e_asinl.c b/sysdeps/ieee754/ldbl-128ibm/e_asinl.c
index 1d42fcecef..0d6e313aa8 100644
--- a/sysdeps/ieee754/ldbl-128ibm/e_asinl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/e_asinl.c
@@ -62,6 +62,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 static const long double
   one = 1.0L,
diff --git a/sysdeps/ieee754/ldbl-128ibm/e_atanhl.c b/sysdeps/ieee754/ldbl-128ibm/e_atanhl.c
index b576f42030..25c286b8ff 100644
--- a/sysdeps/ieee754/ldbl-128ibm/e_atanhl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/e_atanhl.c
@@ -31,6 +31,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 static const long double one = 1.0L, huge = 1e300L;
 
diff --git a/sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c b/sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c
index 0ce702491e..84ea7ee0f5 100644
--- a/sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c
+++ b/sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c
@@ -20,6 +20,7 @@
 
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <float.h>
 
 /* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's
diff --git a/sysdeps/ieee754/ldbl-128ibm/e_hypotl.c b/sysdeps/ieee754/ldbl-128ibm/e_hypotl.c
index 7f3b04fde6..842f77b7ed 100644
--- a/sysdeps/ieee754/ldbl-128ibm/e_hypotl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/e_hypotl.c
@@ -44,6 +44,7 @@
 
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 long double
 __ieee754_hypotl(long double x, long double y)
diff --git a/sysdeps/ieee754/ldbl-128ibm/e_j1l.c b/sysdeps/ieee754/ldbl-128ibm/e_j1l.c
index 6130609f2a..5126900f96 100644
--- a/sysdeps/ieee754/ldbl-128ibm/e_j1l.c
+++ b/sysdeps/ieee754/ldbl-128ibm/e_j1l.c
@@ -21,6 +21,7 @@
 #include <errno.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <float.h>
 
 /* 1 / sqrt(pi) */
diff --git a/sysdeps/ieee754/ldbl-128ibm/e_jnl.c b/sysdeps/ieee754/ldbl-128ibm/e_jnl.c
index dde5720f5d..71b3addfba 100644
--- a/sysdeps/ieee754/ldbl-128ibm/e_jnl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/e_jnl.c
@@ -60,6 +60,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 static const long double
   invsqrtpi = 5.6418958354775628694807945156077258584405E-1L,
diff --git a/sysdeps/ieee754/ldbl-128ibm/e_powl.c b/sysdeps/ieee754/ldbl-128ibm/e_powl.c
index b6c7011501..f59ad4e113 100644
--- a/sysdeps/ieee754/ldbl-128ibm/e_powl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/e_powl.c
@@ -66,6 +66,7 @@
 
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 static const long double bp[] = {
   1.0L,
diff --git a/sysdeps/ieee754/ldbl-128ibm/e_sinhl.c b/sysdeps/ieee754/ldbl-128ibm/e_sinhl.c
index 67d9d24ce7..f869fb068c 100644
--- a/sysdeps/ieee754/ldbl-128ibm/e_sinhl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/e_sinhl.c
@@ -31,6 +31,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 static const long double one = 1.0, shuge = 1.0e307;
 
diff --git a/sysdeps/ieee754/ldbl-128ibm/k_sincosl.c b/sysdeps/ieee754/ldbl-128ibm/k_sincosl.c
index 4374220b1b..ba95337c31 100644
--- a/sysdeps/ieee754/ldbl-128ibm/k_sincosl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/k_sincosl.c
@@ -20,6 +20,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 static const long double c[] = {
 #define ONE c[0]
diff --git a/sysdeps/ieee754/ldbl-128ibm/k_sinl.c b/sysdeps/ieee754/ldbl-128ibm/k_sinl.c
index 1590c6db95..46d1d7b52a 100644
--- a/sysdeps/ieee754/ldbl-128ibm/k_sinl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/k_sinl.c
@@ -20,6 +20,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 static const long double c[] = {
 #define ONE c[0]
diff --git a/sysdeps/ieee754/ldbl-128ibm/k_tanl.c b/sysdeps/ieee754/ldbl-128ibm/k_tanl.c
index 232e00c345..3927fca25c 100644
--- a/sysdeps/ieee754/ldbl-128ibm/k_tanl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/k_tanl.c
@@ -59,6 +59,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <libc-diag.h>
 
 static const long double
diff --git a/sysdeps/ieee754/ldbl-128ibm/s_asinhl.c b/sysdeps/ieee754/ldbl-128ibm/s_asinhl.c
index 21ba14bd83..d4977e5414 100644
--- a/sysdeps/ieee754/ldbl-128ibm/s_asinhl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/s_asinhl.c
@@ -28,6 +28,7 @@ static char rcsid[] = "$NetBSD: s_asinh.c,v 1.9 1995/05/12 04:57:37 jtc Exp $";
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <math_ldbl_opt.h>
 
 static const long double
diff --git a/sysdeps/ieee754/ldbl-128ibm/s_atanl.c b/sysdeps/ieee754/ldbl-128ibm/s_atanl.c
index 0560d820ae..32cf36c65d 100644
--- a/sysdeps/ieee754/ldbl-128ibm/s_atanl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/s_atanl.c
@@ -62,6 +62,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <math_ldbl_opt.h>
 
 /* arctan(k/8), k = 0, ..., 82 */
diff --git a/sysdeps/ieee754/ldbl-128ibm/s_erfl.c b/sysdeps/ieee754/ldbl-128ibm/s_erfl.c
index 7b761b0afa..5302fee522 100644
--- a/sysdeps/ieee754/ldbl-128ibm/s_erfl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/s_erfl.c
@@ -105,6 +105,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <math_ldbl_opt.h>
 #include <fix-int-fp-convert-zero.h>
 
diff --git a/sysdeps/ieee754/ldbl-128ibm/s_fmal.c b/sysdeps/ieee754/ldbl-128ibm/s_fmal.c
index 3821fd75cd..08aeb4c0d4 100644
--- a/sysdeps/ieee754/ldbl-128ibm/s_fmal.c
+++ b/sysdeps/ieee754/ldbl-128ibm/s_fmal.c
@@ -21,6 +21,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <math_ldbl_opt.h>
 #include <mul_split.h>
 #include <stdlib.h>
diff --git a/sysdeps/ieee754/ldbl-128ibm/s_nextupl.c b/sysdeps/ieee754/ldbl-128ibm/s_nextupl.c
index 7743fe882f..9b532a3a25 100644
--- a/sysdeps/ieee754/ldbl-128ibm/s_nextupl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/s_nextupl.c
@@ -16,6 +16,7 @@
    License along with the GNU C Library; if not, see
    <http://www.gnu.org/licenses/>.  */
 
+#include <float.h>
 #include <math.h>
 #include <math_private.h>
 #include <math_ldbl_opt.h>
diff --git a/sysdeps/ieee754/ldbl-128ibm/s_tanhl.c b/sysdeps/ieee754/ldbl-128ibm/s_tanhl.c
index e6457a1c1c..3504862402 100644
--- a/sysdeps/ieee754/ldbl-128ibm/s_tanhl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/s_tanhl.c
@@ -41,6 +41,7 @@ static char rcsid[] = "$NetBSD: s_tanh.c,v 1.7 1995/05/10 20:48:22 jtc Exp $";
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <math_ldbl_opt.h>
 
 static const long double one=1.0L, two=2.0L, tiny = 1.0e-300L;
diff --git a/sysdeps/ieee754/ldbl-96/e_asinl.c b/sysdeps/ieee754/ldbl-96/e_asinl.c
index 85982ca158..806906a58a 100644
--- a/sysdeps/ieee754/ldbl-96/e_asinl.c
+++ b/sysdeps/ieee754/ldbl-96/e_asinl.c
@@ -61,6 +61,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 static const long double
   one = 1.0L,
diff --git a/sysdeps/ieee754/ldbl-96/e_atanhl.c b/sysdeps/ieee754/ldbl-96/e_atanhl.c
index af3487f661..c1f43ff3d7 100644
--- a/sysdeps/ieee754/ldbl-96/e_atanhl.c
+++ b/sysdeps/ieee754/ldbl-96/e_atanhl.c
@@ -35,6 +35,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 static const long double one = 1.0, huge = 1e4900L;
 
diff --git a/sysdeps/ieee754/ldbl-96/e_gammal_r.c b/sysdeps/ieee754/ldbl-96/e_gammal_r.c
index 5e5fd100c5..fc7a5c55dc 100644
--- a/sysdeps/ieee754/ldbl-96/e_gammal_r.c
+++ b/sysdeps/ieee754/ldbl-96/e_gammal_r.c
@@ -19,6 +19,7 @@
 
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <float.h>
 
 /* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's
diff --git a/sysdeps/ieee754/ldbl-96/e_hypotl.c b/sysdeps/ieee754/ldbl-96/e_hypotl.c
index 560ef226f4..f664e30c98 100644
--- a/sysdeps/ieee754/ldbl-96/e_hypotl.c
+++ b/sysdeps/ieee754/ldbl-96/e_hypotl.c
@@ -48,6 +48,7 @@
 
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 long double __ieee754_hypotl(long double x, long double y)
 {
diff --git a/sysdeps/ieee754/ldbl-96/e_j1l.c b/sysdeps/ieee754/ldbl-96/e_j1l.c
index 6d0cb9ef67..581615d563 100644
--- a/sysdeps/ieee754/ldbl-96/e_j1l.c
+++ b/sysdeps/ieee754/ldbl-96/e_j1l.c
@@ -75,6 +75,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 static long double pone (long double), qone (long double);
 
diff --git a/sysdeps/ieee754/ldbl-96/e_jnl.c b/sysdeps/ieee754/ldbl-96/e_jnl.c
index fe249b9126..394921f564 100644
--- a/sysdeps/ieee754/ldbl-96/e_jnl.c
+++ b/sysdeps/ieee754/ldbl-96/e_jnl.c
@@ -60,6 +60,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 static const long double
   invsqrtpi = 5.64189583547756286948079e-1L, two = 2.0e0L, one = 1.0e0L;
diff --git a/sysdeps/ieee754/ldbl-96/e_sinhl.c b/sysdeps/ieee754/ldbl-96/e_sinhl.c
index 64f7c91333..a4b39783bc 100644
--- a/sysdeps/ieee754/ldbl-96/e_sinhl.c
+++ b/sysdeps/ieee754/ldbl-96/e_sinhl.c
@@ -39,6 +39,7 @@ static char rcsid[] = "$NetBSD: $";
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 static const long double one = 1.0, shuge = 1.0e4931L;
 
diff --git a/sysdeps/ieee754/ldbl-96/k_sinl.c b/sysdeps/ieee754/ldbl-96/k_sinl.c
index 30a3c2ac0e..2549f71d19 100644
--- a/sysdeps/ieee754/ldbl-96/k_sinl.c
+++ b/sysdeps/ieee754/ldbl-96/k_sinl.c
@@ -23,6 +23,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 
 /* The polynomials have not been optimized for extended-precision and
    may contain more terms than needed.  */
diff --git a/sysdeps/ieee754/ldbl-96/k_tanl.c b/sysdeps/ieee754/ldbl-96/k_tanl.c
index f8641d5ce4..9b5151baa2 100644
--- a/sysdeps/ieee754/ldbl-96/k_tanl.c
+++ b/sysdeps/ieee754/ldbl-96/k_tanl.c
@@ -59,6 +59,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <libc-diag.h>
 
 static const long double
diff --git a/sysdeps/ieee754/ldbl-96/s_asinhl.c b/sysdeps/ieee754/ldbl-96/s_asinhl.c
index 50a08ed90b..2b9ae1f677 100644
--- a/sysdeps/ieee754/ldbl-96/s_asinhl.c
+++ b/sysdeps/ieee754/ldbl-96/s_asinhl.c
@@ -32,6 +32,7 @@ static char rcsid[] = "$NetBSD: $";
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <libm-alias-ldouble.h>
 
 static const long double
diff --git a/sysdeps/ieee754/ldbl-96/s_erfl.c b/sysdeps/ieee754/ldbl-96/s_erfl.c
index 0f8974085b..1e42df70a7 100644
--- a/sysdeps/ieee754/ldbl-96/s_erfl.c
+++ b/sysdeps/ieee754/ldbl-96/s_erfl.c
@@ -108,6 +108,7 @@
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <libm-alias-ldouble.h>
 
 static const long double
diff --git a/sysdeps/ieee754/ldbl-96/s_nextupl.c b/sysdeps/ieee754/ldbl-96/s_nextupl.c
index a1ca1fecee..5dff32ce73 100644
--- a/sysdeps/ieee754/ldbl-96/s_nextupl.c
+++ b/sysdeps/ieee754/ldbl-96/s_nextupl.c
@@ -16,6 +16,7 @@
    License along with the GNU C Library; if not, see
    <http://www.gnu.org/licenses/>.  */
 
+#include <float.h>
 #include <math.h>
 #include <math_private.h>
 #include <libm-alias-ldouble.h>
diff --git a/sysdeps/ieee754/ldbl-96/s_tanhl.c b/sysdeps/ieee754/ldbl-96/s_tanhl.c
index b0578fffb5..b1b3e0637b 100644
--- a/sysdeps/ieee754/ldbl-96/s_tanhl.c
+++ b/sysdeps/ieee754/ldbl-96/s_tanhl.c
@@ -45,6 +45,7 @@ static char rcsid[] = "$NetBSD: $";
 #include <float.h>
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <libm-alias-ldouble.h>
 
 static const long double one=1.0, two=2.0, tiny = 1.0e-4900L;
diff --git a/sysdeps/powerpc/fpu/e_hypot.c b/sysdeps/powerpc/fpu/e_hypot.c
index 58352c0e54..b55db6d259 100644
--- a/sysdeps/powerpc/fpu/e_hypot.c
+++ b/sysdeps/powerpc/fpu/e_hypot.c
@@ -19,6 +19,7 @@
 
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <stdint.h>
 
 static const double two60   = 1.152921504606847e+18;
diff --git a/sysdeps/x86/fpu/powl_helper.c b/sysdeps/x86/fpu/powl_helper.c
index d270d30890..651eedd792 100644
--- a/sysdeps/x86/fpu/powl_helper.c
+++ b/sysdeps/x86/fpu/powl_helper.c
@@ -18,6 +18,7 @@
 
 #include <math.h>
 #include <math_private.h>
+#include <math-underflow.h>
 #include <stdbool.h>
 
 /* High parts and low parts of -log (k/16), for integer k from 12 to