apply random float error to most floating-point operations

4 files changed, 403 insertions, 84 deletions

diff --git a/src/tools/miri/src/intrinsics/mod.rs b/src/tools/miri/src/intrinsics/mod.rs
index ec4fdfe0bac..377f3a902ed 100644
--- a/src/tools/miri/src/intrinsics/mod.rs
+++ b/src/tools/miri/src/intrinsics/mod.rs
@@ -7,12 +7,13 @@ use rand::Rng;
 use rustc_abi::Size;
 use rustc_apfloat::{Float, Round};
 use rustc_middle::mir;
-use rustc_middle::ty::{self, FloatTy};
+use rustc_middle::ty::{self, FloatTy, ScalarInt};
 use rustc_span::{Symbol, sym};
 
 use self::atomic::EvalContextExt as _;
 use self::helpers::{ToHost, ToSoft, check_intrinsic_arg_count};
 use self::simd::EvalContextExt as _;
+use crate::math::apply_random_float_error_ulp;
 use crate::*;
 
 impl<'tcx> EvalContextExt<'tcx> for crate::MiriInterpCx<'tcx> {}
@@ -206,10 +207,26 @@ pub trait EvalContextExt<'tcx>: crate::MiriInterpCxExt<'tcx> {
                 this.write_scalar(res, dest)?;
             }
 
+            "sqrtf32" => {
+                let [f] = check_intrinsic_arg_count(args)?;
+                let f = this.read_scalar(f)?.to_f32()?;
+                // Sqrt is specified to be fully precise.
+                let res = math::sqrt(f);
+                let res = this.adjust_nan(res, &[f]);
+                this.write_scalar(res, dest)?;
+            }
+            "sqrtf64" => {
+                let [f] = check_intrinsic_arg_count(args)?;
+                let f = this.read_scalar(f)?.to_f64()?;
+                // Sqrt is specified to be fully precise.
+                let res = math::sqrt(f);
+                let res = this.adjust_nan(res, &[f]);
+                this.write_scalar(res, dest)?;
+            }
+
             #[rustfmt::skip]
             | "sinf32"
             | "cosf32"
-            | "sqrtf32"
             | "expf32"
             | "exp2f32"
             | "logf32"
@@ -218,26 +235,33 @@ pub trait EvalContextExt<'tcx>: crate::MiriInterpCxExt<'tcx> {
             => {
                 let [f] = check_intrinsic_arg_count(args)?;
                 let f = this.read_scalar(f)?.to_f32()?;
-                // Using host floats except for sqrt (but it's fine, these operations do not have
+                // Using host floats (but it's fine, these operations do not have
                 // guaranteed precision).
+                let host = f.to_host();
                 let res = match intrinsic_name {
-                    "sinf32" => f.to_host().sin().to_soft(),
-                    "cosf32" => f.to_host().cos().to_soft(),
-                    "sqrtf32" => math::sqrt(f),
-                    "expf32" => f.to_host().exp().to_soft(),
-                    "exp2f32" => f.to_host().exp2().to_soft(),
-                    "logf32" => f.to_host().ln().to_soft(),
-                    "log10f32" => f.to_host().log10().to_soft(),
-                    "log2f32" => f.to_host().log2().to_soft(),
+                    "sinf32" => host.sin(),
+                    "cosf32" => host.cos(),
+                    "expf32" => host.exp(),
+                    "exp2f32" => host.exp2(),
+                    "logf32" => host.ln(),
+                    "log10f32" => host.log10(),
+                    "log2f32" => host.log2(),
                     _ => bug!(),
                 };
+                let res = res.to_soft();
+                // Apply a relative error of 16ULP to introduce some non-determinism
+                // simulating imprecise implementations and optimizations.
+                let res = apply_random_float_error_ulp(
+                    this,
+                    res,
+                    4
+                );
                 let res = this.adjust_nan(res, &[f]);
                 this.write_scalar(res, dest)?;
             }
             #[rustfmt::skip]
             | "sinf64"
             | "cosf64"
-            | "sqrtf64"
             | "expf64"
             | "exp2f64"
             | "logf64"
@@ -246,19 +270,27 @@ pub trait EvalContextExt<'tcx>: crate::MiriInterpCxExt<'tcx> {
             => {
                 let [f] = check_intrinsic_arg_count(args)?;
                 let f = this.read_scalar(f)?.to_f64()?;
-                // Using host floats except for sqrt (but it's fine, these operations do not have
+                // Using host floats (but it's fine, these operations do not have
                 // guaranteed precision).
+                let host = f.to_host();
                 let res = match intrinsic_name {
-                    "sinf64" => f.to_host().sin().to_soft(),
-                    "cosf64" => f.to_host().cos().to_soft(),
-                    "sqrtf64" => math::sqrt(f),
-                    "expf64" => f.to_host().exp().to_soft(),
-                    "exp2f64" => f.to_host().exp2().to_soft(),
-                    "logf64" => f.to_host().ln().to_soft(),
-                    "log10f64" => f.to_host().log10().to_soft(),
-                    "log2f64" => f.to_host().log2().to_soft(),
+                    "sinf64" => host.sin(),
+                    "cosf64" => host.cos(),
+                    "expf64" => host.exp(),
+                    "exp2f64" => host.exp2(),
+                    "logf64" => host.ln(),
+                    "log10f64" => host.log10(),
+                    "log2f64" => host.log2(),
                     _ => bug!(),
                 };
+                let res = res.to_soft();
+                // Apply a relative error of 16ULP to introduce some non-determinism
+                // simulating imprecise implementations and optimizations.
+                let res = apply_random_float_error_ulp(
+                    this,
+                    res,
+                    4
+                );
                 let res = this.adjust_nan(res, &[f]);
                 this.write_scalar(res, dest)?;
             }
@@ -316,6 +348,7 @@ pub trait EvalContextExt<'tcx>: crate::MiriInterpCxExt<'tcx> {
             }
 
             "powf32" => {
+                // FIXME: apply random relative error but without altering behaviour of powf
                 let [f1, f2] = check_intrinsic_arg_count(args)?;
                 let f1 = this.read_scalar(f1)?.to_f32()?;
                 let f2 = this.read_scalar(f2)?.to_f32()?;
@@ -325,6 +358,7 @@ pub trait EvalContextExt<'tcx>: crate::MiriInterpCxExt<'tcx> {
                 this.write_scalar(res, dest)?;
             }
             "powf64" => {
+                // FIXME: apply random relative error but without altering behaviour of powf
                 let [f1, f2] = check_intrinsic_arg_count(args)?;
                 let f1 = this.read_scalar(f1)?.to_f64()?;
                 let f2 = this.read_scalar(f2)?.to_f64()?;
@@ -335,6 +369,7 @@ pub trait EvalContextExt<'tcx>: crate::MiriInterpCxExt<'tcx> {
             }
 
             "powif32" => {
+                // FIXME: apply random relative error but without altering behaviour of powi
                 let [f, i] = check_intrinsic_arg_count(args)?;
                 let f = this.read_scalar(f)?.to_f32()?;
                 let i = this.read_scalar(i)?.to_i32()?;
@@ -344,6 +379,7 @@ pub trait EvalContextExt<'tcx>: crate::MiriInterpCxExt<'tcx> {
                 this.write_scalar(res, dest)?;
             }
             "powif64" => {
+                // FIXME: apply random relative error but without altering behaviour of powi
                 let [f, i] = check_intrinsic_arg_count(args)?;
                 let f = this.read_scalar(f)?.to_f64()?;
                 let i = this.read_scalar(i)?.to_i32()?;
@@ -372,7 +408,10 @@ pub trait EvalContextExt<'tcx>: crate::MiriInterpCxExt<'tcx> {
                     _ => bug!(),
                 };
                 let res = this.binary_op(op, &a, &b)?;
-                // `binary_op` already called `generate_nan` if necessary.
+                // `binary_op` already called `generate_nan` if needed.
+                // Apply a relative error of 16ULP to simulate non-deterministic precision loss
+                // due to optimizations.
+                let res = apply_random_float_error_to_imm(this, res)?;
                 this.write_immediate(*res, dest)?;
             }
 
@@ -418,11 +457,14 @@ pub trait EvalContextExt<'tcx>: crate::MiriInterpCxExt<'tcx> {
                     _ => {}
                 }
                 let res = this.binary_op(op, &a, &b)?;
+                // This cannot be a NaN so we also don't have to apply any non-determinism.
+                // (Also, `binary_op` already called `generate_nan` if needed.)
                 if !float_finite(&res)? {
                     throw_ub_format!("`{intrinsic_name}` intrinsic produced non-finite value as result");
                 }
-                // This cannot be a NaN so we also don't have to apply any non-determinism.
-                // (Also, `binary_op` already called `generate_nan` if needed.)
+                // Apply a relative error of 16ULP to simulate non-deterministic precision loss
+                // due to optimizations.
+                let res = apply_random_float_error_to_imm(this, res)?;
                 this.write_immediate(*res, dest)?;
             }
 
@@ -455,3 +497,21 @@ pub trait EvalContextExt<'tcx>: crate::MiriInterpCxExt<'tcx> {
         interp_ok(EmulateItemResult::NeedsReturn)
     }
 }
+
+/// Applies a random 16ULP floating point error to `val` and returns the new value.
+/// 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>,
+) -> 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(), 4).into(),
+        ty::Float(FloatTy::F32) => apply_random_float_error_ulp(ecx, scalar.to_f32(), 4).into(),
+        ty::Float(FloatTy::F64) => apply_random_float_error_ulp(ecx, scalar.to_f64(), 4).into(),
+        ty::Float(FloatTy::F128) => apply_random_float_error_ulp(ecx, scalar.to_f128(), 4).into(),
+        _ => bug!("intrinsic called with non-float input type"),
+    };
+
+    interp_ok(ImmTy::from_scalar_int(res, val.layout))
+}
diff --git a/src/tools/miri/src/math.rs b/src/tools/miri/src/math.rs
index 7117f722fee..fdd021f8539 100644
--- a/src/tools/miri/src/math.rs
+++ b/src/tools/miri/src/math.rs
@@ -27,6 +27,22 @@ pub(crate) fn apply_random_float_error<F: rustc_apfloat::Float>(
     (val * (F::from_u128(1).value + 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.
+pub(crate) fn apply_random_float_error_ulp<F: rustc_apfloat::Float>(
+    ecx: &mut crate::MiriInterpCx<'_>,
+    val: F,
+    ulp_exponent: 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)
+}
+
 pub(crate) fn sqrt<S: rustc_apfloat::ieee::Semantics>(x: IeeeFloat<S>) -> IeeeFloat<S> {
     match x.category() {
         // preserve zero sign
diff --git a/src/tools/miri/src/shims/foreign_items.rs b/src/tools/miri/src/shims/foreign_items.rs
index bedc1ebdc95..f7746ca81f2 100644
--- a/src/tools/miri/src/shims/foreign_items.rs
+++ b/src/tools/miri/src/shims/foreign_items.rs
@@ -765,7 +765,13 @@ trait EvalContextExtPriv<'tcx>: crate::MiriInterpCxExt<'tcx> {
                     "erfcf" => f_host.erfc(),
                     _ => bug!(),
                 };
-                let res = res.to_soft();
+                // Apply a relative error of 16ULP to introduce some non-determinism
+                // simulating imprecise implementations and optimizations.
+                let res = math::apply_random_float_error_ulp(
+                    this,
+                    res.to_soft(),
+                    4
+                );
                 let res = this.adjust_nan(res, &[f]);
                 this.write_scalar(res, dest)?;
             }
@@ -788,6 +794,13 @@ trait EvalContextExtPriv<'tcx>: crate::MiriInterpCxExt<'tcx> {
                     "fdimf" => f1.to_host().abs_sub(f2.to_host()).to_soft(),
                     _ => bug!(),
                 };
+                // Apply a relative error of 16ULP to introduce some non-determinism
+                // simulating imprecise implementations and optimizations.
+                let res = math::apply_random_float_error_ulp(
+                    this,
+                    res,
+                    4
+                );
                 let res = this.adjust_nan(res, &[f1, f2]);
                 this.write_scalar(res, dest)?;
             }
@@ -826,7 +839,13 @@ trait EvalContextExtPriv<'tcx>: crate::MiriInterpCxExt<'tcx> {
                     "erfc" => f_host.erfc(),
                     _ => bug!(),
                 };
-                let res = res.to_soft();
+                // Apply a relative error of 16ULP to introduce some non-determinism
+                // simulating imprecise implementations and optimizations.
+                let res = math::apply_random_float_error_ulp(
+                    this,
+                    res.to_soft(),
+                    4
+                );
                 let res = this.adjust_nan(res, &[f]);
                 this.write_scalar(res, dest)?;
             }
@@ -849,6 +868,13 @@ trait EvalContextExtPriv<'tcx>: crate::MiriInterpCxExt<'tcx> {
                     "fdim" => f1.to_host().abs_sub(f2.to_host()).to_soft(),
                     _ => bug!(),
                 };
+                // Apply a relative error of 16ULP to introduce some non-determinism
+                // simulating imprecise implementations and optimizations.
+                let res = math::apply_random_float_error_ulp(
+                    this,
+                    res,
+                    4
+                );
                 let res = this.adjust_nan(res, &[f1, f2]);
                 this.write_scalar(res, dest)?;
             }
@@ -874,7 +900,10 @@ trait EvalContextExtPriv<'tcx>: crate::MiriInterpCxExt<'tcx> {
                 // Using host floats (but it's fine, these operations do not have guaranteed precision).
                 let (res, sign) = x.to_host().ln_gamma();
                 this.write_int(sign, &signp)?;
-                let res = this.adjust_nan(res.to_soft(), &[x]);
+                // Apply a relative error of 16ULP to introduce some non-determinism
+                // simulating imprecise implementations and optimizations.
+                let res = math::apply_random_float_error_ulp(this, res.to_soft(), 4);
+                let res = this.adjust_nan(res, &[x]);
                 this.write_scalar(res, dest)?;
             }
             "lgamma_r" => {
@@ -885,7 +914,10 @@ trait EvalContextExtPriv<'tcx>: crate::MiriInterpCxExt<'tcx> {
                 // Using host floats (but it's fine, these operations do not have guaranteed precision).
                 let (res, sign) = x.to_host().ln_gamma();
                 this.write_int(sign, &signp)?;
-                let res = this.adjust_nan(res.to_soft(), &[x]);
+                // Apply a relative error of 16ULP to introduce some non-determinism
+                // simulating imprecise implementations and optimizations.
+                let res = math::apply_random_float_error_ulp(this, res.to_soft(), 4);
+                let res = this.adjust_nan(res, &[x]);
                 this.write_scalar(res, dest)?;
             }
 
diff --git a/src/tools/miri/tests/pass/float.rs b/src/tools/miri/tests/pass/float.rs
index 51cafbb3f43..808495c73d3 100644
--- a/src/tools/miri/tests/pass/float.rs
+++ b/src/tools/miri/tests/pass/float.rs
@@ -13,11 +13,85 @@ use std::fmt::{Debug, Display, LowerHex};
 use std::hint::black_box;
 use std::{f32, f64};
 
+/// Another way of checking if 2 floating-point numbers are almost equal to eachother.
+/// Using `a` and `b` as floating-point numbers:
+///
+/// Instead of performing a simple EPSILON check (which we used at first):
+/// The absolute difference between 'a' and 'b' must not be greater than some E (10^-6, ...)
+///
+/// We will now use ULP: `Units in the Last Place` or `Units of Least Precision`,
+/// more specific, the difference in ULP of `a` and `b`.
+/// First: The ULP of a float 'a' is the smallest possible change at 'a', so the ULP difference represents how
+/// many discrete floating-point steps are needed to reach 'b' from 'a'.
+///
+/// ULP(a) is the distance between the 2 closest floating-point numbers `x` and `y` around `a`, satisfying x < a < y, x != y.
+/// To use this to calculate the ULP difference we have to halve it (we need it at `a`, but we just went "up" and "down", halving it gives us this ULP).
+/// Then take the difference of `b` and `a` and divide it by that ULP and finally round it.
+/// We know now how many floating-point changes we have to apply to `a` to get to `b`.
+///
+/// So if this ULP difference is less than or equal to our chosen upper bound
+/// we can say that `a` and `b` are approximately equal, because they lie "close" enough to each other to be considered equal.
+///
+/// Note: We can see that checking `a` and `b` with different signs has no meaning, but we should not forget
+/// -0.0 and +0.0.
+///
+/// Essentially ULP can be seen as a distance metric of floating-point numbers, but with
+/// the same amount of "spacing" between all consecutive representable values. So even though 2 very large floating point numbers
+/// have a large value difference, their ULP can still be 1, so they are still "approximatly equal",
+/// but the EPSILON check would have failed.
+///
+fn approx_eq_check<F: Float>(
+    actual: F,
+    expected: F,
+    allowed_ulp: F::Int,
+) -> Result<(), NotApproxEq<F>>
+where
+    F::Int: PartialOrd,
+{
+    let actual_signum = actual.signum();
+    let expected_signum = expected.signum();
+
+    if actual_signum != expected_signum {
+        // Floats with different signs must both be 0.
+        if actual != expected {
+            return Err(NotApproxEq::SignsDiffer);
+        }
+    } else {
+        let ulp = (expected.next_up() - expected.next_down()).halve();
+        let ulp_diff = ((actual - expected) / ulp).round().as_int();
+
+        if ulp_diff > allowed_ulp {
+            return Err(NotApproxEq::UlpFail(ulp_diff));
+        }
+    }
+    Ok(())
+}
+
+/// Give more context to execution and result of [`approx_eq_check`].
+enum NotApproxEq<F: Float> {
+    SignsDiffer,
+
+    /// Contains the actual ulp value calculated.
+    UlpFail(F::Int),
+}
+
 macro_rules! assert_approx_eq {
-    ($a:expr, $b:expr) => {{
-        let (a, b) = (&$a, &$b);
-        assert!((*a - *b).abs() < 1.0e-6, "{} is not approximately equal to {}", *a, *b);
+    ($a:expr, $b:expr, $ulp:expr) => {{
+        let (a, b) = ($a, $b);
+        let allowed_ulp = $ulp;
+        match approx_eq_check(a, b, allowed_ulp) {
+            Err(NotApproxEq::SignsDiffer) =>
+                panic!("{a:?} is not approximately equal to {b:?}: signs differ"),
+            Err(NotApproxEq::UlpFail(actual_ulp)) =>
+                panic!("{a:?} is not approximately equal to {b:?}\nulp diff: {actual_ulp} > {allowed_ulp}"),
+            Ok(_) => {}
+        };
     }};
+
+    ($a:expr, $b: expr) => {
+        // accept up to 64ULP (16ULP for host floats and 16ULP for miri artificial error and 32 for any rounding errors)
+        assert_approx_eq!($a, $b, 64);
+    };
 }
 
 fn main() {
@@ -27,15 +101,23 @@ fn main() {
     ops();
     nan_casts();
     rounding();
-    mul_add();
     libm();
+    mul_add();
     test_fast();
     test_algebraic();
     test_fmuladd();
     test_min_max_nondet();
+    test_non_determinism();
 }
 
-trait Float: Copy + PartialEq + Debug {
+trait Float:
+    Copy
+    + PartialEq
+    + Debug
+    + std::ops::Sub<Output = Self>
+    + std::cmp::PartialOrd
+    + std::ops::Div<Output = Self>
+{
     /// The unsigned integer with the same bit width as this float
     type Int: Copy + PartialEq + LowerHex + Debug;
     const BITS: u32 = size_of::<Self>() as u32 * 8;
@@ -49,6 +131,15 @@ trait Float: Copy + PartialEq + Debug {
     const EXPONENT_BIAS: u32 = Self::EXPONENT_SAT >> 1;
 
     fn to_bits(self) -> Self::Int;
+
+    // to make "approx_eq_check" generic
+    fn signum(self) -> Self;
+    fn next_up(self) -> Self;
+    fn next_down(self) -> Self;
+    fn round(self) -> Self;
+    // self / 2
+    fn halve(self) -> Self;
+    fn as_int(self) -> Self::Int;
 }
 
 macro_rules! impl_float {
@@ -61,6 +152,27 @@ macro_rules! impl_float {
             fn to_bits(self) -> Self::Int {
                 self.to_bits()
             }
+
+            fn signum(self) -> Self {
+                self.signum()
+            }
+            fn next_up(self) -> Self {
+                self.next_up()
+            }
+            fn next_down(self) -> Self {
+                self.next_down()
+            }
+            fn round(self) -> Self {
+                self.round()
+            }
+
+            fn halve(self) -> Self {
+                self / 2.0
+            }
+
+            fn as_int(self) -> Self::Int {
+                self as Self::Int
+            }
         }
     };
 }
@@ -1005,8 +1117,8 @@ pub fn libm() {
 
     #[allow(deprecated)]
     {
-        assert_approx_eq!(5.0f32.abs_sub(3.0), 2.0);
-        assert_approx_eq!(3.0f64.abs_sub(5.0), 0.0);
+        assert_approx_eq!(5.0f32.abs_sub(3.0), 2.0f32);
+        assert_approx_eq!(3.0f64.abs_sub(5.0), 0.0f64);
     }
 
     assert_approx_eq!(27.0f32.cbrt(), 3.0f32);
@@ -1023,30 +1135,30 @@ pub fn libm() {
 
     assert_approx_eq!(0f32.sin(), 0f32);
     assert_approx_eq!((f64::consts::PI / 2f64).sin(), 1f64);
-    assert_approx_eq!(f32::consts::FRAC_PI_6.sin(), 0.5);
-    assert_approx_eq!(f64::consts::FRAC_PI_6.sin(), 0.5);
+    assert_approx_eq!(f32::consts::FRAC_PI_6.sin(), 0.5f32);
+    assert_approx_eq!(f64::consts::FRAC_PI_6.sin(), 0.5f64);
     assert_approx_eq!(f32::consts::FRAC_PI_4.sin().asin(), f32::consts::FRAC_PI_4);
     assert_approx_eq!(f64::consts::FRAC_PI_4.sin().asin(), f64::consts::FRAC_PI_4);
 
     assert_approx_eq!(1.0f32.sinh(), 1.1752012f32);
-    assert_approx_eq!(1.0f64.sinh(), 1.1752012f64);
+    assert_approx_eq!(1.0f64.sinh(), 1.1752011936438014f64);
     assert_approx_eq!(2.0f32.asinh(), 1.443635475178810342493276740273105f32);
     assert_approx_eq!((-2.0f64).asinh(), -1.443635475178810342493276740273105f64);
 
     assert_approx_eq!(0f32.cos(), 1f32);
     assert_approx_eq!((f64::consts::PI * 2f64).cos(), 1f64);
-    assert_approx_eq!(f32::consts::FRAC_PI_3.cos(), 0.5);
-    assert_approx_eq!(f64::consts::FRAC_PI_3.cos(), 0.5);
+    assert_approx_eq!(f32::consts::FRAC_PI_3.cos(), 0.5f32);
+    assert_approx_eq!(f64::consts::FRAC_PI_3.cos(), 0.5f64);
     assert_approx_eq!(f32::consts::FRAC_PI_4.cos().acos(), f32::consts::FRAC_PI_4);
     assert_approx_eq!(f64::consts::FRAC_PI_4.cos().acos(), f64::consts::FRAC_PI_4);
 
     assert_approx_eq!(1.0f32.cosh(), 1.54308f32);
-    assert_approx_eq!(1.0f64.cosh(), 1.54308f64);
+    assert_approx_eq!(1.0f64.cosh(), 1.5430806348152437f64);
     assert_approx_eq!(2.0f32.acosh(), 1.31695789692481670862504634730796844f32);
     assert_approx_eq!(3.0f64.acosh(), 1.76274717403908605046521864995958461f64);
 
     assert_approx_eq!(1.0f32.tan(), 1.557408f32);
-    assert_approx_eq!(1.0f64.tan(), 1.557408f64);
+    assert_approx_eq!(1.0f64.tan(), 1.5574077246549023f64);
     assert_approx_eq!(1.0_f32, 1.0_f32.tan().atan());
     assert_approx_eq!(1.0_f64, 1.0_f64.tan().atan());
     assert_approx_eq!(1.0f32.atan2(2.0f32), 0.46364761f32);
@@ -1063,8 +1175,8 @@ pub fn libm() {
     assert_approx_eq!(0.5f32.atanh(), 0.54930614433405484569762261846126285f32);
     assert_approx_eq!(0.5f64.atanh(), 0.54930614433405484569762261846126285f64);
 
-    assert_approx_eq!(5.0f32.gamma(), 24.0);
-    assert_approx_eq!(5.0f64.gamma(), 24.0);
+    assert_approx_eq!(5.0f32.gamma(), 24.0f32);
+    assert_approx_eq!(5.0f64.gamma(), 24.0f64);
     assert_approx_eq!((-0.5f32).gamma(), (-2.0) * f32::consts::PI.sqrt());
     assert_approx_eq!((-0.5f64).gamma(), (-2.0) * f64::consts::PI.sqrt());
 
@@ -1091,11 +1203,11 @@ fn test_fast() {
     pub fn test_operations_f16(a: f16, b: f16) {
         // make sure they all map to the correct operation
         unsafe {
-            assert_eq!(fadd_fast(a, b), a + b);
-            assert_eq!(fsub_fast(a, b), a - b);
-            assert_eq!(fmul_fast(a, b), a * b);
-            assert_eq!(fdiv_fast(a, b), a / b);
-            assert_eq!(frem_fast(a, b), a % b);
+            assert_approx_eq!(fadd_fast(a, b), a + b);
+            assert_approx_eq!(fsub_fast(a, b), a - b);
+            assert_approx_eq!(fmul_fast(a, b), a * b);
+            assert_approx_eq!(fdiv_fast(a, b), a / b);
+            assert_approx_eq!(frem_fast(a, b), a % b);
         }
     }
 
@@ -1103,11 +1215,11 @@ fn test_fast() {
     pub fn test_operations_f32(a: f32, b: f32) {
         // make sure they all map to the correct operation
         unsafe {
-            assert_eq!(fadd_fast(a, b), a + b);
-            assert_eq!(fsub_fast(a, b), a - b);
-            assert_eq!(fmul_fast(a, b), a * b);
-            assert_eq!(fdiv_fast(a, b), a / b);
-            assert_eq!(frem_fast(a, b), a % b);
+            assert_approx_eq!(fadd_fast(a, b), a + b);
+            assert_approx_eq!(fsub_fast(a, b), a - b);
+            assert_approx_eq!(fmul_fast(a, b), a * b);
+            assert_approx_eq!(fdiv_fast(a, b), a / b);
+            assert_approx_eq!(frem_fast(a, b), a % b);
         }
     }
 
@@ -1115,11 +1227,11 @@ fn test_fast() {
     pub fn test_operations_f64(a: f64, b: f64) {
         // make sure they all map to the correct operation
         unsafe {
-            assert_eq!(fadd_fast(a, b), a + b);
-            assert_eq!(fsub_fast(a, b), a - b);
-            assert_eq!(fmul_fast(a, b), a * b);
-            assert_eq!(fdiv_fast(a, b), a / b);
-            assert_eq!(frem_fast(a, b), a % b);
+            assert_approx_eq!(fadd_fast(a, b), a + b);
+            assert_approx_eq!(fsub_fast(a, b), a - b);
+            assert_approx_eq!(fmul_fast(a, b), a * b);
+            assert_approx_eq!(fdiv_fast(a, b), a / b);
+            assert_approx_eq!(frem_fast(a, b), a % b);
         }
     }
 
@@ -1127,11 +1239,11 @@ fn test_fast() {
     pub fn test_operations_f128(a: f128, b: f128) {
         // make sure they all map to the correct operation
         unsafe {
-            assert_eq!(fadd_fast(a, b), a + b);
-            assert_eq!(fsub_fast(a, b), a - b);
-            assert_eq!(fmul_fast(a, b), a * b);
-            assert_eq!(fdiv_fast(a, b), a / b);
-            assert_eq!(frem_fast(a, b), a % b);
+            assert_approx_eq!(fadd_fast(a, b), a + b);
+            assert_approx_eq!(fsub_fast(a, b), a - b);
+            assert_approx_eq!(fmul_fast(a, b), a * b);
+            assert_approx_eq!(fdiv_fast(a, b), a / b);
+            assert_approx_eq!(frem_fast(a, b), a % b);
         }
     }
 
@@ -1153,41 +1265,41 @@ fn test_algebraic() {
     #[inline(never)]
     pub fn test_operations_f16(a: f16, b: f16) {
         // make sure they all map to the correct operation
-        assert_eq!(fadd_algebraic(a, b), a + b);
-        assert_eq!(fsub_algebraic(a, b), a - b);
-        assert_eq!(fmul_algebraic(a, b), a * b);
-        assert_eq!(fdiv_algebraic(a, b), a / b);
-        assert_eq!(frem_algebraic(a, b), a % b);
+        assert_approx_eq!(fadd_algebraic(a, b), a + b);
+        assert_approx_eq!(fsub_algebraic(a, b), a - b);
+        assert_approx_eq!(fmul_algebraic(a, b), a * b);
+        assert_approx_eq!(fdiv_algebraic(a, b), a / b);
+        assert_approx_eq!(frem_algebraic(a, b), a % b);
     }
 
     #[inline(never)]
     pub fn test_operations_f32(a: f32, b: f32) {
         // make sure they all map to the correct operation
-        assert_eq!(fadd_algebraic(a, b), a + b);
-        assert_eq!(fsub_algebraic(a, b), a - b);
-        assert_eq!(fmul_algebraic(a, b), a * b);
-        assert_eq!(fdiv_algebraic(a, b), a / b);
-        assert_eq!(frem_algebraic(a, b), a % b);
+        assert_approx_eq!(fadd_algebraic(a, b), a + b);
+        assert_approx_eq!(fsub_algebraic(a, b), a - b);
+        assert_approx_eq!(fmul_algebraic(a, b), a * b);
+        assert_approx_eq!(fdiv_algebraic(a, b), a / b);
+        assert_approx_eq!(frem_algebraic(a, b), a % b);
     }
 
     #[inline(never)]
     pub fn test_operations_f64(a: f64, b: f64) {
         // make sure they all map to the correct operation
-        assert_eq!(fadd_algebraic(a, b), a + b);
-        assert_eq!(fsub_algebraic(a, b), a - b);
-        assert_eq!(fmul_algebraic(a, b), a * b);
-        assert_eq!(fdiv_algebraic(a, b), a / b);
-        assert_eq!(frem_algebraic(a, b), a % b);
+        assert_approx_eq!(fadd_algebraic(a, b), a + b);
+        assert_approx_eq!(fsub_algebraic(a, b), a - b);
+        assert_approx_eq!(fmul_algebraic(a, b), a * b);
+        assert_approx_eq!(fdiv_algebraic(a, b), a / b);
+        assert_approx_eq!(frem_algebraic(a, b), a % b);
     }
 
     #[inline(never)]
     pub fn test_operations_f128(a: f128, b: f128) {
         // make sure they all map to the correct operation
-        assert_eq!(fadd_algebraic(a, b), a + b);
-        assert_eq!(fsub_algebraic(a, b), a - b);
-        assert_eq!(fmul_algebraic(a, b), a * b);
-        assert_eq!(fdiv_algebraic(a, b), a / b);
-        assert_eq!(frem_algebraic(a, b), a % b);
+        assert_approx_eq!(fadd_algebraic(a, b), a + b);
+        assert_approx_eq!(fsub_algebraic(a, b), a - b);
+        assert_approx_eq!(fmul_algebraic(a, b), a * b);
+        assert_approx_eq!(fdiv_algebraic(a, b), a / b);
+        assert_approx_eq!(frem_algebraic(a, b), a % b);
     }
 
     test_operations_f16(11., 2.);
@@ -1245,3 +1357,102 @@ fn test_min_max_nondet() {
     ensure_both(|| f128::min(0.0, -0.0).is_sign_positive());
     ensure_both(|| f128::max(0.0, -0.0).is_sign_positive());
 }
+
+fn test_non_determinism() {
+    use std::intrinsics::{
+        fadd_algebraic, fadd_fast, fdiv_algebraic, fdiv_fast, fmul_algebraic, fmul_fast,
+        frem_algebraic, frem_fast, fsub_algebraic, fsub_fast,
+    };
+    use std::{f32, f64};
+    // TODO: Also test powi and powf when the non-determinism is implemented for them
+
+    /// Ensure that the operation is non-deterministic
+    #[track_caller]
+    fn ensure_nondet<T: PartialEq + std::fmt::Debug>(f: impl Fn() -> T) {
+        let rounds = 16;
+        let first = f();
+        for _ in 1..rounds {
+            if f() != first {
+                // We saw two different values!
+                return;
+            }
+        }
+        // We saw the same thing N times.
+        panic!("expected non-determinism, got {rounds} times the same result: {first:?}");
+    }
+
+    macro_rules! test_operations_f {
+        ($a:expr, $b:expr) => {
+            ensure_nondet(|| fadd_algebraic($a, $b));
+            ensure_nondet(|| fsub_algebraic($a, $b));
+            ensure_nondet(|| fmul_algebraic($a, $b));
+            ensure_nondet(|| fdiv_algebraic($a, $b));
+            ensure_nondet(|| frem_algebraic($a, $b));
+
+            unsafe {
+                ensure_nondet(|| fadd_fast($a, $b));
+                ensure_nondet(|| fsub_fast($a, $b));
+                ensure_nondet(|| fmul_fast($a, $b));
+                ensure_nondet(|| fdiv_fast($a, $b));
+                ensure_nondet(|| frem_fast($a, $b));
+            }
+        };
+    }
+
+    pub fn test_operations_f16(a: f16, b: f16) {
+        test_operations_f!(a, b);
+    }
+    pub fn test_operations_f32(a: f32, b: f32) {
+        test_operations_f!(a, b);
+        ensure_nondet(|| a.log(b));
+        ensure_nondet(|| a.exp());
+        ensure_nondet(|| 10f32.exp2());
+        ensure_nondet(|| f32::consts::E.ln());
+        ensure_nondet(|| 1f32.ln_1p());
+        ensure_nondet(|| 10f32.log10());
+        ensure_nondet(|| 8f32.log2());
+        ensure_nondet(|| 27.0f32.cbrt());
+        ensure_nondet(|| 3.0f32.hypot(4.0f32));
+        ensure_nondet(|| 1f32.sin());
+        ensure_nondet(|| 0f32.cos());
+        ensure_nondet(|| 1.0f32.sinh());
+        ensure_nondet(|| 1.0f32.asinh());
+        ensure_nondet(|| 1.0f32.cosh());
+        ensure_nondet(|| 2.0f32.acosh());
+        ensure_nondet(|| 1.0f32.tan());
+        ensure_nondet(|| 1.0f32.tanh());
+        ensure_nondet(|| 1.0f32.atan2(2.0f32));
+        ensure_nondet(|| 0.5f32.atanh());
+        ensure_nondet(|| 5.0f32.gamma());
+    }
+    pub fn test_operations_f64(a: f64, b: f64) {
+        test_operations_f!(a, b);
+        ensure_nondet(|| a.log(b));
+        ensure_nondet(|| a.exp());
+        ensure_nondet(|| 50f64.exp2());
+        ensure_nondet(|| 3f64.ln());
+        ensure_nondet(|| 1f64.ln_1p());
+        ensure_nondet(|| f64::consts::E.log10());
+        ensure_nondet(|| f64::consts::E.log2());
+        ensure_nondet(|| 1f64.sin());
+        ensure_nondet(|| 0f64.cos());
+        ensure_nondet(|| 27.0f64.cbrt());
+        ensure_nondet(|| 3.0f64.hypot(4.0f64));
+        ensure_nondet(|| 1.0f64.sinh());
+        ensure_nondet(|| 1.0f64.asinh());
+        ensure_nondet(|| 1.0f64.cosh());
+        ensure_nondet(|| 3.0f64.acosh());
+        ensure_nondet(|| 1.0f64.tan());
+        ensure_nondet(|| 1.0f64.tanh());
+        ensure_nondet(|| 0.5f64.atanh());
+        ensure_nondet(|| 5.0f64.gamma());
+    }
+    pub fn test_operations_f128(a: f128, b: f128) {
+        test_operations_f!(a, b);
+    }
+
+    test_operations_f16(5., 7.);
+    test_operations_f32(12., 5.);
+    test_operations_f64(19., 11.);
+    test_operations_f128(25., 18.);
+}