From: Joseph Myers <joseph@codesourcery.com>
Date: Thu, 10 May 2018 00:53:04 +0000 (+0000)
Subject: Move math_check_force_underflow macros to separate math-underflow.h.
X-Git-Tag: glibc-2.28~291
X-Git-Url: http://git.ipfire.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8f5b00d375dbd7f5e15e57b24fec3bd5a4b1e98d;p=thirdparty%2Fglibc.git

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.
---

diff --git a/ChangeLog b/ChangeLog
index 7ab225e320e..b2b66020e81 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 94a21054c2c..30ed4af49bb 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 f7b679e55fb..451f05c2791 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 00000000000..c5e5c79768d
--- /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 c3894c5a0a5..19f4e11cb47 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 fdeeb589e15..245a1a6412a 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 6840c4819e3..1664f0690a3 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 e099e42a4c3..5fdab1e5ac3 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 82a52faebfe..a947502e94a 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 8aa1f74c8b1..1248cac5932 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 d6af614177f..5f95e92ee00 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 6466ae80c96..95f28c2b669 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 dcf14c75370..b1153906b90 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 af898d8f810..ab11c252990 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 9527ab4ccc5..bd1292239fe 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 f6faf718958..9cd0941431b 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 768cbb9ec18..6bf56945a65 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 608afb5a8c3..df9c8786b27 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 a944abeee89..045cbbb88d6 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 a3750f9f311..2744549cbdd 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 a89075169b4..a2c33cc4edb 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 d0e387e1a21..734f3ca64a5 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 6ef0fbea1c2..9181b22bb8d 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 542d03a7e36..96d5b23ccc6 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 49c24fb4894..c4e34211ac0 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 5577ccebd92..192ff8594d9 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 ff5718c82ec..38db092d045 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 48dfca33176..5f820604f8c 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 3e136e7ee07..8ef74f2196e 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 b7cf5ce3214..c01f1be93b3 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 db13a5793bc..d37a2b66579 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 ba1dbe27b66..b369ac9f5ba 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 c7460371e44..1d8d44befe0 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 05a25eb060a..04ff8b6717b 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 321bf441422..673a97102de 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 55eb1449c57..e03073b8bd6 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 598000a8cd7..43f11ad1f05 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 d6405104489..8b23add3479 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 8aadc4c31a9..887b5f7d1ed 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 a4414ce27ec..cd15ed7d4be 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 17b3663dde7..20f7db81ea7 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 a195d594667..dcf3c35358f 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 9f0e55860f0..228ece20752 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 27ecdd56fa7..0812b54dca4 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 03a4cfd6a55..f61b0f5ef64 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 1d00abdbf1c..693934d59ab 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 03589707784..3b04ce77518 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 009084ee19e..9d2ef899756 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 7ba8c0b79de..87ec7ba21f4 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 c2b96cfa2b7..cc3d63d35d6 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 53c2f1fada7..6021facdfb1 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 4c8c2e22e6a..13b7683e179 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 c4c61cbed4d..735dc9b8d02 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 95b9fe06776..e2730b8fc20 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 dd82b3aec83..7bafd4ae296 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 bff3a5e88cc..e6f46f52e9c 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 635b4e7e429..7739eec2919 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 cce5b6e4fd8..39e7cf3084a 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 36668dd0a32..6eb33ae375e 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 a9f38c7f0ab..9c19fb1bf2f 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 e79023c69ad..8da794d80cb 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 92a482658d3..a733ee76958 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 e05368df538..022ccf4c45b 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 88e91c702dd..c0abfe2926c 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 ea0d29c4931..66881af01a5 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 b8b2ffeba1f..f3181dbc8e6 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 5628c176809..20827e6d098 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 8f62f96c503..fc309dab0a2 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 1d42fcecef8..0d6e313aa83 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 b576f42030a..25c286b8ff7 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 0ce702491e4..84ea7ee0f5d 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 7f3b04fde68..842f77b7ed7 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 6130609f2a7..5126900f964 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 dde5720f5d9..71b3addfba0 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 b6c7011501b..f59ad4e113b 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 67d9d24ce77..f869fb068cf 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 4374220b1b8..ba95337c313 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 1590c6db958..46d1d7b52a8 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 232e00c345f..3927fca25c3 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 21ba14bd836..d4977e54146 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 0560d820ae4..32cf36c65d5 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 7b761b0afa8..5302fee5221 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 3821fd75cd9..08aeb4c0d47 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 7743fe882fa..9b532a3a255 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 e6457a1c1cf..3504862402e 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 85982ca1586..806906a58a4 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 af3487f661c..c1f43ff3d75 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 5e5fd100c58..fc7a5c55dc6 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 560ef226f46..f664e30c987 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 6d0cb9ef677..581615d563f 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 fe249b9126e..394921f5640 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 64f7c913331..a4b39783bcc 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 30a3c2ac0e6..2549f71d194 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 f8641d5ce41..9b5151baa2b 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 50a08ed90bb..2b9ae1f6770 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 0f8974085b5..1e42df70a7e 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 a1ca1fecee0..5dff32ce73f 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 b0578fffb5a..b1b3e0637b7 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 58352c0e54c..b55db6d2590 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 d270d308908..651eedd7920 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