fix mangitude of applied float error

6 files changed, 103 insertions, 69 deletions

diff --git a/src/tools/miri/src/intrinsics/mod.rs b/src/tools/miri/src/intrinsics/mod.rs
index b5e81460773..5e46768b0e6 100644
--- a/src/tools/miri/src/intrinsics/mod.rs
+++ b/src/tools/miri/src/intrinsics/mod.rs
@@ -10,7 +10,7 @@ use rustc_abi::Size;
 use rustc_apfloat::ieee::{IeeeFloat, Semantics};
 use rustc_apfloat::{self, Float, Round};
 use rustc_middle::mir;
-use rustc_middle::ty::{self, FloatTy, ScalarInt};
+use rustc_middle::ty::{self, FloatTy};
 use rustc_span::{Symbol, sym};
 
 use self::atomic::EvalContextExt as _;
@@ -230,7 +230,7 @@ pub trait EvalContextExt<'tcx>: crate::MiriInterpCxExt<'tcx> {
                     let res = apply_random_float_error_ulp(
                         this,
                         res,
-                        2, // log2(4)
+                        4,
                     );
 
                     // Clamp the result to the guaranteed range of this function according to the C standard,
@@ -274,7 +274,7 @@ pub trait EvalContextExt<'tcx>: crate::MiriInterpCxExt<'tcx> {
                     let res = apply_random_float_error_ulp(
                         this,
                         res,
-                        2, // log2(4)
+                        4,
                     );
 
                     // Clamp the result to the guaranteed range of this function according to the C standard,
@@ -336,9 +336,7 @@ pub trait EvalContextExt<'tcx>: crate::MiriInterpCxExt<'tcx> {
 
                     // Apply a relative error of 4ULP to introduce some non-determinism
                     // simulating imprecise implementations and optimizations.
-                    apply_random_float_error_ulp(
-                        this, res, 2, // log2(4)
-                    )
+                    apply_random_float_error_ulp(this, res, 4)
                 });
                 let res = this.adjust_nan(res, &[f1, f2]);
                 this.write_scalar(res, dest)?;
@@ -354,9 +352,7 @@ pub trait EvalContextExt<'tcx>: crate::MiriInterpCxExt<'tcx> {
 
                     // Apply a relative error of 4ULP to introduce some non-determinism
                     // simulating imprecise implementations and optimizations.
-                    apply_random_float_error_ulp(
-                        this, res, 2, // log2(4)
-                    )
+                    apply_random_float_error_ulp(this, res, 4)
                 });
                 let res = this.adjust_nan(res, &[f1, f2]);
                 this.write_scalar(res, dest)?;
@@ -373,9 +369,7 @@ pub trait EvalContextExt<'tcx>: crate::MiriInterpCxExt<'tcx> {
 
                     // Apply a relative error of 4ULP to introduce some non-determinism
                     // simulating imprecise implementations and optimizations.
-                    apply_random_float_error_ulp(
-                        this, res, 2, // log2(4)
-                    )
+                    apply_random_float_error_ulp(this, res, 4)
                 });
                 let res = this.adjust_nan(res, &[f]);
                 this.write_scalar(res, dest)?;
@@ -391,9 +385,7 @@ pub trait EvalContextExt<'tcx>: crate::MiriInterpCxExt<'tcx> {
 
                     // Apply a relative error of 4ULP to introduce some non-determinism
                     // simulating imprecise implementations and optimizations.
-                    apply_random_float_error_ulp(
-                        this, res, 2, // log2(4)
-                    )
+                    apply_random_float_error_ulp(this, res, 4)
                 });
                 let res = this.adjust_nan(res, &[f]);
                 this.write_scalar(res, dest)?;
@@ -448,7 +440,7 @@ pub trait EvalContextExt<'tcx>: crate::MiriInterpCxExt<'tcx> {
                 }
                 // Apply a relative error of 4ULP to simulate non-deterministic precision loss
                 // due to optimizations.
-                let res = apply_random_float_error_to_imm(this, res, 2 /* log2(4) */)?;
+                let res = crate::math::apply_random_float_error_to_imm(this, res, 4)?;
                 this.write_immediate(*res, dest)?;
             }
 
@@ -486,31 +478,6 @@ pub trait EvalContextExt<'tcx>: crate::MiriInterpCxExt<'tcx> {
     }
 }
 
-/// Applies a random ULP floating point error to `val` and returns the new value.
-/// So if you want an X ULP error, `ulp_exponent` should be log2(X).
-///
-/// Will fail if `val` is not a floating point number.
-fn apply_random_float_error_to_imm<'tcx>(
-    ecx: &mut MiriInterpCx<'tcx>,
-    val: ImmTy<'tcx>,
-    ulp_exponent: u32,
-) -> InterpResult<'tcx, ImmTy<'tcx>> {
-    let scalar = val.to_scalar_int()?;
-    let res: ScalarInt = match val.layout.ty.kind() {
-        ty::Float(FloatTy::F16) =>
-            apply_random_float_error_ulp(ecx, scalar.to_f16(), ulp_exponent).into(),
-        ty::Float(FloatTy::F32) =>
-            apply_random_float_error_ulp(ecx, scalar.to_f32(), ulp_exponent).into(),
-        ty::Float(FloatTy::F64) =>
-            apply_random_float_error_ulp(ecx, scalar.to_f64(), ulp_exponent).into(),
-        ty::Float(FloatTy::F128) =>
-            apply_random_float_error_ulp(ecx, scalar.to_f128(), ulp_exponent).into(),
-        _ => bug!("intrinsic called with non-float input type"),
-    };
-
-    interp_ok(ImmTy::from_scalar_int(res, val.layout))
-}
-
 /// For the intrinsics:
 /// - sinf32, sinf64
 /// - cosf32, cosf64
diff --git a/src/tools/miri/src/machine.rs b/src/tools/miri/src/machine.rs
index 7802a58f112..e0ee9541948 100644
--- a/src/tools/miri/src/machine.rs
+++ b/src/tools/miri/src/machine.rs
@@ -1278,7 +1278,7 @@ impl<'tcx> Machine<'tcx> for MiriMachine<'tcx> {
         ecx: &mut InterpCx<'tcx, Self>,
         val: ImmTy<'tcx>,
     ) -> InterpResult<'tcx, ImmTy<'tcx>> {
-        crate::math::apply_random_float_error_to_imm(ecx, val, 2 /* log2(4) */)
+        crate::math::apply_random_float_error_to_imm(ecx, val, 4)
     }
 
     #[inline(always)]
diff --git a/src/tools/miri/src/math.rs b/src/tools/miri/src/math.rs
index dd994b03f5d..604d6c0fd73 100644
--- a/src/tools/miri/src/math.rs
+++ b/src/tools/miri/src/math.rs
@@ -6,10 +6,6 @@ use rustc_middle::ty::{self, FloatTy, ScalarInt};
 use crate::*;
 
 /// Disturbes a floating-point result by a relative error in the range (-2^scale, 2^scale).
-///
-/// For a 2^N ULP error, you can use an `err_scale` of `-(F::PRECISION - 1 - N)`.
-/// In other words, a 1 ULP (absolute) error is the same as a `2^-(F::PRECISION-1)` relative error.
-/// (Subtracting 1 compensates for the integer bit.)
 pub(crate) fn apply_random_float_error<F: rustc_apfloat::Float>(
     ecx: &mut crate::MiriInterpCx<'_>,
     val: F,
@@ -17,11 +13,13 @@ pub(crate) fn apply_random_float_error<F: rustc_apfloat::Float>(
 ) -> F {
     if !ecx.machine.float_nondet
         || matches!(ecx.machine.float_rounding_error, FloatRoundingErrorMode::None)
+        // relative errors don't do anything to zeros... avoid messing up the sign
+        || val.is_zero()
     {
         return val;
     }
-
     let rng = ecx.machine.rng.get_mut();
+
     // Generate a random integer in the range [0, 2^PREC).
     // (When read as binary, the position of the first `1` determines the exponent,
     // and the remaining bits fill the mantissa. `PREC` is one plus the size of the mantissa,
@@ -37,43 +35,66 @@ pub(crate) fn apply_random_float_error<F: rustc_apfloat::Float>(
     let err = r.scalbn(err_scale.strict_sub(F::PRECISION.try_into().unwrap()));
     // give it a random sign
     let err = if rng.random() { -err } else { err };
-    // multiple the value with (1+err)
-    (val * (F::from_u128(1).value + err).value).value
+    // Compute `val*(1+err)`, distributed out as `val + val*err` to avoid the imprecise addition
+    // error being amplified by multiplication.
+    (val + (val * err).value).value
 }
 
-/// [`apply_random_float_error`] gives instructions to apply a 2^N ULP error.
-/// This function implements these instructions such that applying a 2^N ULP error is less error prone.
-/// So for a 2^N ULP error, you would pass N as the `ulp_exponent` argument.
+/// Applies an error of `[-N, +N]` ULP to the given value.
 pub(crate) fn apply_random_float_error_ulp<F: rustc_apfloat::Float>(
     ecx: &mut crate::MiriInterpCx<'_>,
     val: F,
-    ulp_exponent: u32,
+    max_error: u32,
 ) -> F {
-    let n = i32::try_from(ulp_exponent)
-        .expect("`err_scale_for_ulp`: exponent is too large to create an error scale");
-    // we know this fits
-    let prec = i32::try_from(F::PRECISION).unwrap();
-    let err_scale = -(prec - n - 1);
-    apply_random_float_error(ecx, val, err_scale)
+    // We could try to be clever and reuse `apply_random_float_error`, but that is hard to get right
+    // (see <https://github.com/rust-lang/miri/pull/4558#discussion_r2316838085> for why) so we
+    // implement the logic directly instead.
+    if !ecx.machine.float_nondet
+        || matches!(ecx.machine.float_rounding_error, FloatRoundingErrorMode::None)
+        // FIXME: also disturb zeros? That requires a lot more cases in `fixed_float_value`
+        // and might make the std test suite quite unhappy.
+        || val.is_zero()
+    {
+        return val;
+    }
+    let rng = ecx.machine.rng.get_mut();
+
+    let max_error = i64::from(max_error);
+    let error = match ecx.machine.float_rounding_error {
+        FloatRoundingErrorMode::Random => rng.random_range(-max_error..=max_error),
+        FloatRoundingErrorMode::Max =>
+            if rng.random() {
+                max_error
+            } else {
+                -max_error
+            },
+        FloatRoundingErrorMode::None => unreachable!(),
+    };
+    // If upwards ULP and downwards ULP differ, we take the average.
+    let ulp = (((val.next_up().value - val).value + (val - val.next_down().value).value).value
+        / F::from_u128(2).value)
+        .value;
+    // Shift the value by N times the ULP
+    (val + (ulp * F::from_i128(error.into()).value).value).value
 }
 
-/// Applies a random 16ULP floating point error to `val` and returns the new value.
+/// Applies an error of `[-N, +N]` ULP to the given value.
 /// Will fail if `val` is not a floating point number.
 pub(crate) fn apply_random_float_error_to_imm<'tcx>(
     ecx: &mut MiriInterpCx<'tcx>,
     val: ImmTy<'tcx>,
-    ulp_exponent: u32,
+    max_error: u32,
 ) -> InterpResult<'tcx, ImmTy<'tcx>> {
     let scalar = val.to_scalar_int()?;
     let res: ScalarInt = match val.layout.ty.kind() {
         ty::Float(FloatTy::F16) =>
-            apply_random_float_error_ulp(ecx, scalar.to_f16(), ulp_exponent).into(),
+            apply_random_float_error_ulp(ecx, scalar.to_f16(), max_error).into(),
         ty::Float(FloatTy::F32) =>
-            apply_random_float_error_ulp(ecx, scalar.to_f32(), ulp_exponent).into(),
+            apply_random_float_error_ulp(ecx, scalar.to_f32(), max_error).into(),
         ty::Float(FloatTy::F64) =>
-            apply_random_float_error_ulp(ecx, scalar.to_f64(), ulp_exponent).into(),
+            apply_random_float_error_ulp(ecx, scalar.to_f64(), max_error).into(),
         ty::Float(FloatTy::F128) =>
-            apply_random_float_error_ulp(ecx, scalar.to_f128(), ulp_exponent).into(),
+            apply_random_float_error_ulp(ecx, scalar.to_f128(), max_error).into(),
         _ => bug!("intrinsic called with non-float input type"),
     };
 
diff --git a/src/tools/miri/tests/pass/float.rs b/src/tools/miri/tests/pass/float.rs
index fe7316c6680..2be262d76a4 100644
--- a/src/tools/miri/tests/pass/float.rs
+++ b/src/tools/miri/tests/pass/float.rs
@@ -39,9 +39,8 @@ macro_rules! assert_approx_eq {
     }};
 
     ($a:expr, $b: expr) => {
-        // accept up to 12ULP (4ULP for host floats and 4ULP for miri artificial error and 4 for any additional effects
-        // due to having multiple error sources.
-        assert_approx_eq!($a, $b, 12);
+        // accept up to 8ULP (4ULP for host floats and 4ULP for miri artificial error).
+        assert_approx_eq!($a, $b, 8);
     };
 }
 
@@ -176,6 +175,7 @@ fn assert_eq_msg<T: PartialEq + Debug>(x: T, y: T, msg: impl Display) {
 }
 
 /// Check that floats have bitwise equality
+#[track_caller]
 fn assert_biteq<F: Float>(a: F, b: F, msg: impl Display) {
     let ab = a.to_bits();
     let bb = b.to_bits();
@@ -189,6 +189,7 @@ fn assert_biteq<F: Float>(a: F, b: F, msg: impl Display) {
 }
 
 /// Check that two floats have equality
+#[track_caller]
 fn assert_feq<F: Float>(a: F, b: F, msg: impl Display) {
     let ab = a.to_bits();
     let bb = b.to_bits();
diff --git a/src/tools/miri/tests/pass/float_extra_rounding_error.rs b/src/tools/miri/tests/pass/float_extra_rounding_error.rs
index 73512399f82..24d7cf2ccee 100644
--- a/src/tools/miri/tests/pass/float_extra_rounding_error.rs
+++ b/src/tools/miri/tests/pass/float_extra_rounding_error.rs
@@ -9,7 +9,7 @@ use std::hint::black_box;
 
 fn main() {
     let expected = cfg_select! {
-        random => 13, // FIXME: why is it 13?
+        random => 9, // -4 ..= +4 ULP error
         max => 2,
         none => 1,
     };
@@ -20,4 +20,12 @@ fn main() {
         values.insert(val.to_bits());
     }
     assert_eq!(values.len(), expected);
+
+    if !cfg!(none) {
+        // Ensure the smallest and biggest value are 8 ULP apart.
+        // We can just subtract the raw bit representations for this.
+        let min = *values.iter().min().unwrap();
+        let max = *values.iter().max().unwrap();
+        assert_eq!(min.abs_diff(max), 8);
+    }
 }
diff --git a/src/tools/miri/tests/pass/shims/x86/rounding-error.rs b/src/tools/miri/tests/pass/shims/x86/rounding-error.rs
new file mode 100644
index 00000000000..bf56111b2e4
--- /dev/null
+++ b/src/tools/miri/tests/pass/shims/x86/rounding-error.rs
@@ -0,0 +1,37 @@
+// We're testing x86 target specific features
+//@only-target: x86_64 i686
+
+//! rsqrt and rcp SSE/AVX operations are approximate. We use that as license to treat them as
+//! non-deterministic. Ensure that we do indeed see random results within the expected error bounds.
+
+#[cfg(target_arch = "x86")]
+use std::arch::x86::*;
+#[cfg(target_arch = "x86_64")]
+use std::arch::x86_64::*;
+use std::collections::HashSet;
+
+fn main() {
+    let mut vals = HashSet::new();
+    for _ in 0..50 {
+        unsafe {
+            // Compute the inverse square root of 4.0, four times.
+            let a = _mm_setr_ps(4.0, 4.0, 4.0, 4.0);
+            let exact = 0.5;
+            let r = _mm_rsqrt_ps(a);
+            let r: [f32; 4] = std::mem::transmute(r);
+            // Check the results.
+            for r in r {
+                vals.insert(r.to_bits());
+                // Ensure the relative error is less than 2^-12.
+                let rel_error = (r - exact) / exact;
+                let log_error = rel_error.abs().log2();
+                assert!(
+                    rel_error == 0.0 || log_error < -12.0,
+                    "got an error of {rel_error} = 2^{log_error}"
+                );
+            }
+        }
+    }
+    // Ensure we saw a bunch of different results.
+    assert!(vals.len() >= 50);
+}