Start removing force_eval uses

src/math/asinh.rs

@@ -26,10 +26,6 @@ pub fn asinh(mut x: f64) -> f64 {
     } else if e >= 0x3ff - 26 {
         /* |x| >= 0x1p-26, up to 1.6ulp error in [0.125,0.5] */
         x = log1p(x + x * x / (sqrt(x * x + 1.0) + 1.0));
-    } else {
-        /* |x| < 0x1p-26, raise inexact if x != 0 */
-        let x1p120 = f64::from_bits(0x4770000000000000);
-        force_eval!(x + x1p120);
     }
 
     if sign {

src/math/asinhf.rs

@@ -25,10 +25,6 @@ pub fn asinhf(mut x: f32) -> f32 {
     } else if i >= 0x3f800000 - (12 << 23) {
         /* |x| >= 0x1p-12, up to 1.6ulp error in [0.125,0.5] */
         x = log1pf(x + x * x / (sqrtf(x * x + 1.0) + 1.0));
-    } else {
-        /* |x| < 0x1p-12, raise inexact if x!=0 */
-        let x1p120 = f32::from_bits(0x7b800000);
-        force_eval!(x + x1p120);
     }
 
     if sign {

src/math/atan.rs

@@ -83,11 +83,6 @@ pub fn atan(x: f64) -> f64 {
         /* |x| < 0.4375 */
         if ix < 0x3e40_0000 {
             /* |x| < 2^-27 */
-            if ix < 0x0010_0000 {
-                /* raise underflow for subnormal x */
-                force_eval!(x as f32);
-            }
-
             return x;
         }
 

src/math/atanf.rs

@@ -63,10 +63,6 @@ pub fn atanf(mut x: f32) -> f32 {
         /* |x| < 0.4375 */
         if ix < 0x39800000 {
             /* |x| < 2**-12 */
-            if ix < 0x00800000 {
-                /* raise underflow for subnormal x */
-                force_eval!(x * x);
-            }
             return x;
         }
         -1

src/math/atanh.rs

@@ -15,12 +15,7 @@ pub fn atanh(x: f64) -> f64 {
     let mut y = f64::from_bits(u & 0x7fff_ffff_ffff_ffff);
 
     if e < 0x3ff - 1 {
-        if e < 0x3ff - 32 {
-            /* handle underflow */
-            if e == 0 {
-                force_eval!(y as f32);
-            }
-        } else {
+        if e >= 0x3ff - 32 {
             /* |x| < 0.5, up to 1.7ulp error */
             y = 0.5 * log1p(2.0 * y + 2.0 * y * y / (1.0 - y));
         }

src/math/atanhf.rs

@@ -15,12 +15,7 @@ pub fn atanhf(mut x: f32) -> f32 {
     x = f32::from_bits(u);
 
     if u < 0x3f800000 - (1 << 23) {
-        if u < 0x3f800000 - (32 << 23) {
-            /* handle underflow */
-            if u < (1 << 23) {
-                force_eval!((x * x) as f32);
-            }
-        } else {
+        if u >= 0x3f800000 - (32 << 23) {
             /* |x| < 0.5, up to 1.7ulp error */
             x = 0.5 * log1pf(2.0 * x + 2.0 * x * x / (1.0 - x));
         }

src/math/cosf.rs

@@ -28,18 +28,13 @@ const C4_PIO2: f64 = 4. * FRAC_PI_2; /* 0x401921FB, 0x54442D18 */
 pub fn cosf(x: f32) -> f32 {
     let x64 = x as f64;
 
-    let x1p120 = f32::from_bits(0x7b800000); // 0x1p120f === 2 ^ 120
-
     let mut ix = x.to_bits();
     let sign = (ix >> 31) != 0;
     ix &= 0x7fffffff;
 
     if ix <= 0x3f490fda {
         /* |x| ~<= pi/4 */
         if ix < 0x39800000 {
-            /* |x| < 2**-12 */
-            /* raise inexact if x != 0 */
-            force_eval!(x + x1p120);
             return 1.;
         }
         return k_cosf(x64);

src/math/cosh.rs

@@ -18,8 +18,6 @@ pub fn cosh(mut x: f64) -> f64 {
     /* |x| < log(2) */
     if w < 0x3fe62e42 {
         if w < 0x3ff00000 - (26 << 20) {
-            let x1p120 = f64::from_bits(0x4770000000000000);
-            force_eval!(x + x1p120);
             return 1.;
         }
         let t = expm1(x); // exponential minus 1

src/math/coshf.rs

@@ -9,8 +9,6 @@ use super::k_expo2f;
 /// Angles are specified in radians.
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
 pub fn coshf(mut x: f32) -> f32 {
-    let x1p120 = f32::from_bits(0x7b800000); // 0x1p120f === 2 ^ 120
-
     /* |x| */
     let mut ix = x.to_bits();
     ix &= 0x7fffffff;
@@ -20,7 +18,6 @@ pub fn coshf(mut x: f32) -> f32 {
     /* |x| < log(2) */
     if w < 0x3f317217 {
         if w < (0x3f800000 - (12 << 23)) {
-            force_eval!(x + x1p120);
             return 1.;
         }
         let t = expm1f(x);

src/math/exp.rs

@@ -84,7 +84,6 @@ const P5: f64 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
 pub fn exp(mut x: f64) -> f64 {
     let x1p1023 = f64::from_bits(0x7fe0000000000000); // 0x1p1023 === 2 ^ 1023
-    let x1p_149 = f64::from_bits(0x36a0000000000000); // 0x1p-149 === 2 ^ -149
 
     let hi: f64;
     let lo: f64;
@@ -106,13 +105,10 @@ pub fn exp(mut x: f64) -> f64 {
             return x;
         }
         if x > 709.782712893383973096 {
-            /* overflow if x!=inf */
             x *= x1p1023;
             return x;
         }
         if x < -708.39641853226410622 {
-            /* underflow if x!=-inf */
-            force_eval!((-x1p_149 / x) as f32);
             if x < -745.13321910194110842 {
                 return 0.;
             }
@@ -137,8 +133,6 @@ pub fn exp(mut x: f64) -> f64 {
         hi = x;
         lo = 0.;
     } else {
-        /* inexact if x!=0 */
-        force_eval!(x1p1023 + x);
         return 1. + x;
     }
 

src/math/exp2.rs

@@ -336,7 +336,6 @@ pub fn exp2(mut x: f64) -> f64 {
     // union {double f; uint64_t i;} u = {x};
     // union {uint32_t u; int32_t i;} k;
     let x1p1023 = f64::from_bits(0x7fe0000000000000);
-    let x1p52 = f64::from_bits(0x4330000000000000);
     let _0x1p_149 = f64::from_bits(0xb6a0000000000000);
 
     /* Filter out exceptional cases. */
@@ -356,10 +355,6 @@ pub fn exp2(mut x: f64) -> f64 {
         }
         if ui >> 63 != 0 {
             /* x <= -1022 */
-            /* underflow */
-            if x <= -1075.0 || x - x1p52 + x1p52 != x {
-                force_eval!((_0x1p_149 / x) as f32);
-            }
             if x <= -1075.0 {
                 return 0.0;
             }

src/math/exp2f.rs

@@ -102,9 +102,6 @@ pub fn exp2f(mut x: f32) -> f32 {
         }
         if ui >= 0x80000000 {
             /* x < -126 */
-            if ui >= 0xc3160000 || (ui & 0x0000ffff != 0) {
-                force_eval!(f32::from_bits(0x80000001) / x);
-            }
             if ui >= 0xc3160000 {
                 /* x <= -150 */
                 return 0.0;

src/math/expf.rs

@@ -33,7 +33,6 @@ const P2: f32 = -2.7667332906e-3; /* -0xb55215.0p-32 */
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
 pub fn expf(mut x: f32) -> f32 {
     let x1p127 = f32::from_bits(0x7f000000); // 0x1p127f === 2 ^ 127
-    let x1p_126 = f32::from_bits(0x800000); // 0x1p-126f === 2 ^ -126  /*original 0x1p-149f    ??????????? */
     let mut hx = x.to_bits();
     let sign = (hx >> 31) as i32; /* sign bit of x */
     let signb: bool = sign != 0;
@@ -53,8 +52,6 @@ pub fn expf(mut x: f32) -> f32 {
             return x;
         }
         if signb {
-            /* underflow */
-            force_eval!(-x1p_126 / x);
             if hx >= 0x42cff1b5 {
                 /* x <= -103.972084f */
                 return 0.;
@@ -84,8 +81,6 @@ pub fn expf(mut x: f32) -> f32 {
         hi = x;
         lo = 0.;
     } else {
-        /* raise inexact */
-        force_eval!(x1p127 + x);
         return 1. + x;
     }
 

src/math/expm1.rs

@@ -82,9 +82,6 @@ pub fn expm1(mut x: f64) -> f64 {
         c = (hi - x) - lo;
     } else if hx < 0x3c900000 {
         /* |x| < 2**-54, return x */
-        if hx < 0x00100000 {
-            force_eval!(x);
-        }
         return x;
     } else {
         c = 0.0;

src/math/expm1f.rs

@@ -84,9 +84,6 @@ pub fn expm1f(mut x: f32) -> f32 {
         c = (hi - x) - lo;
     } else if hx < 0x33000000 {
         /* when |x|<2**-25, return x */
-        if hx < 0x00800000 {
-            force_eval!(x * x);
-        }
         return x;
     } else {
         k = 0;

src/math/floor.rs

@@ -48,7 +48,12 @@ pub fn floor(x: f64) -> f64 {
     };
     /* special case because of non-nearest rounding modes */
     if e < 0x3ff {
-        force_eval!(y);
+        // FIXME: Using the black_box via force_eval! here prevents optimizations that prove this
+        // code is free of panics. Perhaps that means that the read_volatile is an insufficient
+        // optimization barrier. But we have no tests for that.
+        unsafe {
+            core::ptr::read_volatile(&y);
+        }
         return if (ui >> 63) != 0 { -1. } else { 0. };
     }
     if y > 0. {

src/math/floorf.rs

@@ -24,13 +24,11 @@ pub fn floorf(x: f32) -> f32 {
         if (ui & m) == 0 {
             return x;
         }
-        force_eval!(x + f32::from_bits(0x7b800000));
         if ui >> 31 != 0 {
             ui += m;
         }
         ui &= !m;
     } else {
-        force_eval!(x + f32::from_bits(0x7b800000));
         if ui >> 31 == 0 {
             ui = 0;
         } else if ui << 1 != 0 {

src/math/fmaf.rs

@@ -26,7 +26,6 @@
  */
 
 use core::f32;
-use core::ptr::read_volatile;
 
 use super::fenv::{
     feclearexcept, fegetround, feraiseexcept, fetestexcept, FE_INEXACT, FE_TONEAREST, FE_UNDERFLOW,
@@ -74,7 +73,7 @@ pub fn fmaf(x: f32, y: f32, mut z: f32) -> f32 {
         if e < 0x3ff - 126 && e >= 0x3ff - 149 && fetestexcept(FE_INEXACT) != 0 {
             feclearexcept(FE_INEXACT);
             // prevent `xy + vz` from being CSE'd with `xy + z` above
-            let vz: f32 = unsafe { read_volatile(&z) };
+            let vz = force_eval!(z);
             result = xy + vz as f64;
             if fetestexcept(FE_INEXACT) != 0 {
                 feraiseexcept(FE_UNDERFLOW);

src/math/ilogb.rs

@@ -9,7 +9,6 @@ pub fn ilogb(x: f64) -> i32 {
     if e == 0 {
         i <<= 12;
         if i == 0 {
-            force_eval!(0.0 / 0.0);
             return FP_ILOGB0;
         }
         /* subnormal x */
@@ -20,7 +19,6 @@ pub fn ilogb(x: f64) -> i32 {
         }
         e
     } else if e == 0x7ff {
-        force_eval!(0.0 / 0.0);
         if (i << 12) != 0 {
             FP_ILOGBNAN
         } else {

src/math/ilogbf.rs

@@ -9,7 +9,6 @@ pub fn ilogbf(x: f32) -> i32 {
     if e == 0 {
         i <<= 9;
         if i == 0 {
-            force_eval!(0.0 / 0.0);
             return FP_ILOGB0;
         }
         /* subnormal x */
@@ -20,7 +19,6 @@ pub fn ilogbf(x: f32) -> i32 {
         }
         e
     } else if e == 0xff {
-        force_eval!(0.0 / 0.0);
         if (i << 9) != 0 {
             FP_ILOGBNAN
         } else {

src/math/log1p.rs

@@ -95,10 +95,6 @@ pub fn log1p(x: f64) -> f64 {
         }
         if hx << 1 < 0x3ca00000 << 1 {
             /* |x| < 2**-53 */
-            /* underflow if subnormal */
-            if (hx & 0x7ff00000) == 0 {
-                force_eval!(x as f32);
-            }
             return x;
         }
         if hx <= 0xbfd2bec4 {

src/math/log1pf.rs

@@ -50,10 +50,6 @@ pub fn log1pf(x: f32) -> f32 {
         }
         if ix << 1 < 0x33800000 << 1 {
             /* |x| < 2**-24 */
-            /* underflow if subnormal */
-            if (ix & 0x7f800000) == 0 {
-                force_eval!(x * x);
-            }
             return x;
         }
         if ix <= 0xbe95f619 {

src/math/mod.rs

@@ -1,6 +1,6 @@
 macro_rules! force_eval {
     ($e:expr) => {
-        unsafe { ::core::ptr::read_volatile(&$e) }
+        core::hint::black_box($e)
     };
 }
 

src/math/nextafter.rs

@@ -23,15 +23,6 @@ pub fn nextafter(x: f64, y: f64) -> f64 {
         ux_i += 1;
     }
 
-    let e = ux_i >> 52 & 0x7ff;
-    // raise overflow if ux.f is infinite and x is finite
-    if e == 0x7ff {
-        force_eval!(x + x);
-    }
     let ux_f = f64::from_bits(ux_i);
-    // raise underflow if ux.f is subnormal or zero
-    if e == 0 {
-        force_eval!(x * x + ux_f * ux_f);
-    }
     ux_f
 }

src/math/nextafterf.rs

@@ -23,15 +23,6 @@ pub fn nextafterf(x: f32, y: f32) -> f32 {
         ux_i += 1;
     }
 
-    let e = ux_i & 0x7f80_0000_u32;
-    // raise overflow if ux_f is infinite and x is finite
-    if e == 0x7f80_0000_u32 {
-        force_eval!(x + x);
-    }
     let ux_f = f32::from_bits(ux_i);
-    // raise underflow if ux_f is subnormal or zero
-    if e == 0 {
-        force_eval!(x * x + ux_f * ux_f);
-    }
     ux_f
 }

src/math/sin.rs

@@ -42,21 +42,13 @@ use super::{k_cos, k_sin, rem_pio2};
 //      TRIG(x) returns trig(x) nearly rounded
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
 pub fn sin(x: f64) -> f64 {
-    let x1p120 = f64::from_bits(0x4770000000000000); // 0x1p120f === 2 ^ 120
-
     /* High word of x. */
     let ix = (f64::to_bits(x) >> 32) as u32 & 0x7fffffff;
 
     /* |x| ~< pi/4 */
     if ix <= 0x3fe921fb {
         if ix < 0x3e500000 {
             /* |x| < 2**-26 */
-            /* raise inexact if x != 0 and underflow if subnormal*/
-            if ix < 0x00100000 {
-                force_eval!(x / x1p120);
-            } else {
-                force_eval!(x + x1p120);
-            }
             return x;
         }
         return k_sin(x, 0.0, 0);

src/math/sincos.rs

@@ -25,13 +25,6 @@ pub fn sincos(x: f64) -> (f64, f64) {
     if ix <= 0x3fe921fb {
         /* if |x| < 2**-27 * sqrt(2) */
         if ix < 0x3e46a09e {
-            /* raise inexact if x!=0 and underflow if subnormal */
-            let x1p120 = f64::from_bits(0x4770000000000000); // 0x1p120 == 2^120
-            if ix < 0x00100000 {
-                force_eval!(x / x1p120);
-            } else {
-                force_eval!(x + x1p120);
-            }
             return (x, 1.0);
         }
         return (k_sin(x, 0.0, 0), k_cos(x, 0.0));