diff --git a/doc/release/upcoming_changes/19479.compatibility.rst b/doc/release/upcoming_changes/19479.compatibility.rst
new file mode 100644
index 000000000000..83533a305749
--- /dev/null
+++ b/doc/release/upcoming_changes/19479.compatibility.rst
@@ -0,0 +1,7 @@
+Distutils forces strict floating point model on clang
+-----------------------------------------------------
+NumPy now sets the ``-ftrapping-math`` option on clang to enforce correct
+floating point error handling for universal functions.
+Clang defaults to non-IEEE and C99 conform behaviour otherwise.
+This change (using the equivalent but newer ``-ffp-exception-behavior=strict``)
+was attempted in NumPy 1.21, but was effectively never used.
diff --git a/numpy/core/src/umath/loops_exponent_log.dispatch.c.src b/numpy/core/src/umath/loops_exponent_log.dispatch.c.src
index 41e0bf37b6f8..9970ad2ea994 100644
--- a/numpy/core/src/umath/loops_exponent_log.dispatch.c.src
+++ b/numpy/core/src/umath/loops_exponent_log.dispatch.c.src
@@ -149,8 +149,8 @@ fma_get_exponent(__m256 x)
     __m256 normal_mask = _mm256_cmp_ps(x, _mm256_set1_ps(FLT_MIN), _CMP_GE_OQ);
 
     /*
-     * It is necessary for temp1 to be volatile, a bug in clang optimizes it out which leads
-     * to an overflow warning in some cases. See https://github.com/numpy/numpy/issues/18005
+     * The volatile is probably unnecessary now since we compile clang with
+     * `-ftrapping-math`: https://github.com/numpy/numpy/issues/18005
      */
     volatile __m256 temp1 = _mm256_blendv_ps(x, _mm256_set1_ps(0.0f), normal_mask);
     __m256 temp = _mm256_mul_ps(temp1, two_power_100);
@@ -180,8 +180,8 @@ fma_get_mantissa(__m256 x)
     __m256 normal_mask = _mm256_cmp_ps(x, _mm256_set1_ps(FLT_MIN), _CMP_GE_OQ);
 
     /*
-     * It is necessary for temp1 to be volatile, a bug in clang optimizes it out which leads
-     * to an overflow warning in some cases. See https://github.com/numpy/numpy/issues/18005
+     * The volatile is probably unnecessary now since we compile clang with
+     * `-ftrapping-math`: https://github.com/numpy/numpy/issues/18005
      */
     volatile __m256 temp1 = _mm256_blendv_ps(x, _mm256_set1_ps(0.0f), normal_mask);
     __m256 temp = _mm256_mul_ps(temp1, two_power_100);
diff --git a/numpy/distutils/ccompiler.py b/numpy/distutils/ccompiler.py
index 061f4862dc19..6d063ee4e0b8 100644
--- a/numpy/distutils/ccompiler.py
+++ b/numpy/distutils/ccompiler.py
@@ -388,7 +388,8 @@ def CCompiler_customize_cmd(self, cmd, ignore=()):
     if hasattr(self, 'compiler') and 'clang' in self.compiler[0]:
         # clang defaults to a non-strict floating error point model.
         # Since NumPy and most Python libs give warnings for these, override:
-        self.compiler.append('-ffp-exception-behavior=strict')
+        self.compiler.append('-ftrapping-math')
+        self.compiler_so.append('-ftrapping-math')
 
     def allow(attr):
         return getattr(cmd, attr, None) is not None and attr not in ignore