Auto merge of #138062 - LorrensP-2158466:miri-enable-float-nondet, r=RalfJung

Enable Non-determinism of float operations in Miri and change std tests

Links to [#4208](https://github.com/rust-lang/miri/issues/4208) and [#3555](https://github.com/rust-lang/miri/issues/3555) in Miri.

Non-determinism of floating point operations was disabled in rust-lang/rust#137594 because it breaks the tests and doc-tests in core/coretests and std. This PR enables some of them.

This pr includes the following changes:

- Enables the float non-determinism but with a lower relative error of 4ULP instead of 16ULP
- These operations now have a fixed output based on the C23 standard, except the pow operations, this is tracked in [#4286](https://github.com/rust-lang/miri/issues/4286#issue-3010677983)
- Changes tests that made incorrect assumptions about the operations, not to make that assumption anymore (from `assert_eq!` to `assert_approx_eq!`.
- Changed the doctests of the stdlib of these operations to compare against fixed constants instead of `f*::EPSILON`, which now succeed with Miri and `-Zmiri-many-seeds`
- Added a constant `APPROX_DELTA` in `std/tests/floats/f32.rs` which is used for approximation tests, but with a different value when run in Miri. This is to make these tests succeed.
- Added tests in the float tests of Miri to test the C23 behaviour.

Fixes https://github.com/rust-lang/miri/issues/4208

3 files changed, 329 insertions, 62 deletions

diff --git a/src/tools/miri/src/intrinsics/mod.rs b/src/tools/miri/src/intrinsics/mod.rs
index a4882a20148..9957e351ff1 100644
--- a/src/tools/miri/src/intrinsics/mod.rs
+++ b/src/tools/miri/src/intrinsics/mod.rs
@@ -3,17 +3,20 @@
 mod atomic;
 mod simd;
 
+use std::ops::Neg;
+
 use rand::Rng;
 use rustc_abi::Size;
-use rustc_apfloat::{Float, Round};
+use rustc_apfloat::ieee::{IeeeFloat, Semantics};
+use rustc_apfloat::{self, 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_to_imm;
+use crate::math::{IeeeExt, apply_random_float_error_ulp};
 use crate::*;
 
 impl<'tcx> EvalContextExt<'tcx> for crate::MiriInterpCx<'tcx> {}
@@ -187,31 +190,39 @@ 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 (but it's fine, these operations do not have
-                // guaranteed precision).
-                let host = f.to_host();
-                let res = match intrinsic_name {
-                    "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.
-                // FIXME: temporarily disabled as it breaks std tests.
-                // let res = apply_random_float_error_ulp(
-                //     this,
-                //     res,
-                //     4, // log2(16)
-                // );
+
+                let res = fixed_float_value(intrinsic_name, &[f]).unwrap_or_else(||{
+                    // 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" => 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 4ULP to introduce some non-determinism
+                    // simulating imprecise implementations and optimizations.
+                    let res = apply_random_float_error_ulp(
+                        this,
+                        res,
+                        2, // log2(4)
+                    );
+
+                    // Clamp the result to the guaranteed range of this function according to the C standard,
+                    // if any.
+                    clamp_float_value(intrinsic_name, res)
+                });
                 let res = this.adjust_nan(res, &[f]);
                 this.write_scalar(res, dest)?;
             }
+
             #[rustfmt::skip]
             | "sinf64"
             | "cosf64"
@@ -223,28 +234,35 @@ 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 (but it's fine, these operations do not have
-                // guaranteed precision).
-                let host = f.to_host();
-                let res = match intrinsic_name {
-                    "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.
-                // FIXME: temporarily disabled as it breaks std tests.
-                // let res = apply_random_float_error_ulp(
-                //     this,
-                //     res,
-                //     4, // log2(16)
-                // );
+
+                let res = fixed_float_value(intrinsic_name, &[f]).unwrap_or_else(||{
+                    // 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" => 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 4ULP to introduce some non-determinism
+                    // simulating imprecise implementations and optimizations.
+                    let res = apply_random_float_error_ulp(
+                        this,
+                        res,
+                        2, // log2(4)
+                    );
+
+                    // Clamp the result to the guaranteed range of this function according to the C standard,
+                    // if any.
+                    clamp_float_value(intrinsic_name, res)
+                });
                 let res = this.adjust_nan(res, &[f]);
                 this.write_scalar(res, dest)?;
             }
@@ -302,43 +320,75 @@ 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()?;
-                // Using host floats (but it's fine, this operation does not have guaranteed precision).
-                let res = f1.to_host().powf(f2.to_host()).to_soft();
+
+                let res = fixed_float_value(intrinsic_name, &[f1, f2]).unwrap_or_else(|| {
+                    // Using host floats (but it's fine, this operation does not have guaranteed precision).
+                    let res = f1.to_host().powf(f2.to_host()).to_soft();
+
+                    // 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)
+                    )
+                });
                 let res = this.adjust_nan(res, &[f1, f2]);
                 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()?;
-                // Using host floats (but it's fine, this operation does not have guaranteed precision).
-                let res = f1.to_host().powf(f2.to_host()).to_soft();
+
+                let res = fixed_float_value(intrinsic_name, &[f1, f2]).unwrap_or_else(|| {
+                    // Using host floats (but it's fine, this operation does not have guaranteed precision).
+                    let res = f1.to_host().powf(f2.to_host()).to_soft();
+
+                    // 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)
+                    )
+                });
                 let res = this.adjust_nan(res, &[f1, f2]);
                 this.write_scalar(res, dest)?;
             }
 
             "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()?;
-                // Using host floats (but it's fine, this operation does not have guaranteed precision).
-                let res = f.to_host().powi(i).to_soft();
+
+                let res = fixed_powi_float_value(f, i).unwrap_or_else(|| {
+                    // Using host floats (but it's fine, this operation does not have guaranteed precision).
+                    let res = f.to_host().powi(i).to_soft();
+
+                    // 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)
+                    )
+                });
                 let res = this.adjust_nan(res, &[f]);
                 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()?;
-                // Using host floats (but it's fine, this operation does not have guaranteed precision).
-                let res = f.to_host().powi(i).to_soft();
+
+                let res = fixed_powi_float_value(f, i).unwrap_or_else(|| {
+                    // Using host floats (but it's fine, this operation does not have guaranteed precision).
+                    let res = f.to_host().powi(i).to_soft();
+
+                    // 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)
+                    )
+                });
                 let res = this.adjust_nan(res, &[f]);
                 this.write_scalar(res, dest)?;
             }
@@ -425,3 +475,97 @@ pub trait EvalContextExt<'tcx>: crate::MiriInterpCxExt<'tcx> {
         interp_ok(EmulateItemResult::NeedsReturn)
     }
 }
+
+/// 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
+/// - expf32, expf64, exp2f32, exp2f64
+/// - logf32, logf64, log2f32, log2f64, log10f32, log10f64
+/// - powf32, powf64
+///
+/// Returns `Some(output)` if the `intrinsic` results in a defined fixed `output` specified in the C standard
+/// (specifically, C23 annex F.10)  when given `args` as arguments. Outputs that are unaffected by a relative error
+/// (such as INF and zero) are not handled here, they are assumed to be handled by the underlying
+/// implementation. Returns `None` if no specific value is guaranteed.
+fn fixed_float_value<S: Semantics>(
+    intrinsic_name: &str,
+    args: &[IeeeFloat<S>],
+) -> Option<IeeeFloat<S>> {
+    let one = IeeeFloat::<S>::one();
+    match (intrinsic_name, args) {
+        // cos(+- 0) = 1
+        ("cosf32" | "cosf64", [input]) if input.is_zero() => Some(one),
+
+        // e^0 = 1
+        ("expf32" | "expf64" | "exp2f32" | "exp2f64", [input]) if input.is_zero() => Some(one),
+
+        // 1^y = 1 for any y, even a NaN.
+        ("powf32" | "powf64", [base, _]) if *base == one => Some(one),
+
+        // (-1)^(±INF) = 1
+        ("powf32" | "powf64", [base, exp]) if *base == -one && exp.is_infinite() => Some(one),
+
+        // FIXME(#4286): The C ecosystem is inconsistent with handling sNaN's, some return 1 others propogate
+        // the NaN. We should return either 1 or the NaN non-deterministically here.
+        // But for now, just handle them all the same.
+        // x^(±0) = 1 for any x, even a NaN
+        ("powf32" | "powf64", [_, exp]) if exp.is_zero() => Some(one),
+
+        // There are a lot of cases for fixed outputs according to the C Standard, but these are mainly INF or zero
+        // which are not affected by the applied error.
+        _ => None,
+    }
+}
+
+/// Returns `Some(output)` if `powi` (called `pown` in C) results in a fixed value specified in the C standard
+/// (specifically, C23 annex F.10.4.6) when doing `base^exp`. Otherwise, returns `None`.
+fn fixed_powi_float_value<S: Semantics>(base: IeeeFloat<S>, exp: i32) -> Option<IeeeFloat<S>> {
+    match (base.category(), exp) {
+        // x^0 = 1, if x is not a Signaling NaN
+        // FIXME(#4286): The C ecosystem is inconsistent with handling sNaN's, some return 1 others propogate
+        // the NaN. We should return either 1 or the NaN non-deterministically here.
+        // But for now, just handle them all the same.
+        (_, 0) => Some(IeeeFloat::<S>::one()),
+
+        _ => None,
+    }
+}
+
+/// Given an floating-point operation and a floating-point value, clamps the result to the output
+/// range of the given operation.
+fn clamp_float_value<S: Semantics>(intrinsic_name: &str, val: IeeeFloat<S>) -> IeeeFloat<S> {
+    match intrinsic_name {
+        // sin and cos: [-1, 1]
+        "sinf32" | "cosf32" | "sinf64" | "cosf64" =>
+            val.clamp(IeeeFloat::<S>::one().neg(), IeeeFloat::<S>::one()),
+        // exp: [0, +INF]
+        "expf32" | "exp2f32" | "expf64" | "exp2f64" =>
+            IeeeFloat::<S>::maximum(val, IeeeFloat::<S>::ZERO),
+        _ => val,
+    }
+}
diff --git a/src/tools/miri/src/math.rs b/src/tools/miri/src/math.rs
index 2ff29c7ac1a..d1355a21684 100644
--- a/src/tools/miri/src/math.rs
+++ b/src/tools/miri/src/math.rs
@@ -151,6 +151,20 @@ pub(crate) fn sqrt<S: rustc_apfloat::ieee::Semantics>(x: IeeeFloat<S>) -> IeeeFl
     }
 }
 
+/// Extend functionality of rustc_apfloat softfloats
+pub trait IeeeExt: rustc_apfloat::Float {
+    #[inline]
+    fn one() -> Self {
+        Self::from_u128(1).value
+    }
+
+    #[inline]
+    fn clamp(self, min: Self, max: Self) -> Self {
+        self.maximum(min).minimum(max)
+    }
+}
+impl<S: rustc_apfloat::ieee::Semantics> IeeeExt for IeeeFloat<S> {}
+
 #[cfg(test)]
 mod tests {
     use rustc_apfloat::ieee::{DoubleS, HalfS, IeeeFloat, QuadS, SingleS};
diff --git a/src/tools/miri/tests/pass/float.rs b/src/tools/miri/tests/pass/float.rs
index 98a88cfd62d..7ce0bc88517 100644
--- a/src/tools/miri/tests/pass/float.rs
+++ b/src/tools/miri/tests/pass/float.rs
@@ -45,6 +45,31 @@ macro_rules! assert_approx_eq {
     };
 }
 
+
+/// From IEEE 754 a Signaling NaN for single precision has the following representation:
+/// ```
+/// s | 1111 1111 | 0x..x
+/// ````
+/// Were at least one `x` is a 1.
+///
+/// This sNaN has the following representation and is used for testing purposes.:
+/// ```
+/// 0 | 1111111 | 01..0
+/// ```
+const SNAN_F32: f32 = f32::from_bits(0x7fa00000);
+
+/// From IEEE 754 a Signaling NaN for double precision has the following representation:
+/// ```
+/// s | 1111 1111 111 | 0x..x
+/// ````
+/// Were at least one `x` is a 1.
+///
+/// This sNaN has the following representation and is used for testing purposes.:
+/// ```
+/// 0 | 1111 1111 111 | 01..0
+/// ```
+const SNAN_F64: f64 = f64::from_bits(0x7ff4000000000000);
+
 fn main() {
     basic();
     casts();
@@ -1008,17 +1033,84 @@ pub fn libm() {
     assert_approx_eq!(25f32.powf(-2f32), 0.0016f32);
     assert_approx_eq!(400f64.powf(0.5f64), 20f64);
 
+    // Some inputs to powf and powi result in fixed outputs
+    // and thus must be exactly equal to that value.
+    // C standard says:
+    // 1^y = 1 for any y, even a NaN.
+    assert_eq!(1f32.powf(10.0), 1.0);
+    assert_eq!(1f64.powf(100.0), 1.0);
+    assert_eq!(1f32.powf(f32::INFINITY), 1.0);
+    assert_eq!(1f64.powf(f64::INFINITY), 1.0);
+    assert_eq!(1f32.powf(f32::NAN), 1.0);
+    assert_eq!(1f64.powf(f64::NAN), 1.0);
+
+    // f*::NAN is a quiet NAN and should return 1 as well.
+    assert_eq!(f32::NAN.powf(0.0), 1.0);
+    assert_eq!(f64::NAN.powf(0.0), 1.0);
+
+    assert_eq!(42f32.powf(0.0), 1.0);
+    assert_eq!(42f64.powf(0.0), 1.0);
+    assert_eq!(f32::INFINITY.powf(0.0), 1.0);
+    assert_eq!(f64::INFINITY.powf(0.0), 1.0);
+
+    // f*::NAN is a quiet NAN and should return 1 as well.
+    assert_eq!(f32::NAN.powi(0), 1.0);
+    assert_eq!(f64::NAN.powi(0), 1.0);
+
+    assert_eq!(10.0f32.powi(0), 1.0);
+    assert_eq!(10.0f64.powi(0), 1.0);
+    assert_eq!(f32::INFINITY.powi(0), 1.0);
+    assert_eq!(f64::INFINITY.powi(0), 1.0);
+
+    assert_eq!((-1f32).powf(f32::INFINITY), 1.0);
+    assert_eq!((-1f64).powf(f64::INFINITY), 1.0);
+    assert_eq!((-1f32).powf(f32::NEG_INFINITY), 1.0);
+    assert_eq!((-1f64).powf(f64::NEG_INFINITY), 1.0);
+
+    // For pow (powf in rust) the C standard says:
+    // x^0 = 1 for all x even a sNaN
+    // FIXME(#4286): this does not match the behavior of all implementations.
+    assert_eq!(SNAN_F32.powf(0.0), 1.0);
+    assert_eq!(SNAN_F64.powf(0.0), 1.0);
+
+    // For pown (powi in rust) the C standard says:
+    // x^0 = 1 for all x even a sNaN
+    // FIXME(#4286): this does not match the behavior of all implementations.
+    assert_eq!(SNAN_F32.powi(0), 1.0);
+    assert_eq!(SNAN_F64.powi(0), 1.0);
+
+    assert_eq!(0f32.powi(10), 0.0);
+    assert_eq!(0f64.powi(100), 0.0);
+    assert_eq!(0f32.powi(9), 0.0);
+    assert_eq!(0f64.powi(99), 0.0);
+
+    assert_biteq((-0f32).powf(10.0), 0.0, "-0^x = +0 where x is positive");
+    assert_biteq((-0f64).powf(100.0), 0.0, "-0^x = +0 where x is positive");
+    assert_biteq((-0f32).powf(9.0), -0.0, "-0^x = -0 where x is negative");
+    assert_biteq((-0f64).powf(99.0), -0.0, "-0^x = -0 where x is negative");
+
+    assert_biteq((-0f32).powi(10), 0.0, "-0^x = +0 where x is positive");
+    assert_biteq((-0f64).powi(100), 0.0, "-0^x = +0 where x is positive");
+    assert_biteq((-0f32).powi(9), -0.0, "-0^x = -0 where x is negative");
+    assert_biteq((-0f64).powi(99), -0.0, "-0^x = -0 where x is negative");
+
     assert_approx_eq!(1f32.exp(), f32::consts::E);
     assert_approx_eq!(1f64.exp(), f64::consts::E);
+    assert_eq!(0f32.exp(), 1.0);
+    assert_eq!(0f64.exp(), 1.0);
 
     assert_approx_eq!(1f32.exp_m1(), f32::consts::E - 1.0);
     assert_approx_eq!(1f64.exp_m1(), f64::consts::E - 1.0);
 
     assert_approx_eq!(10f32.exp2(), 1024f32);
     assert_approx_eq!(50f64.exp2(), 1125899906842624f64);
+    assert_eq!(0f32.exp2(), 1.0);
+    assert_eq!(0f64.exp2(), 1.0);
 
     assert_approx_eq!(f32::consts::E.ln(), 1f32);
-    assert_approx_eq!(1f64.ln(), 0f64);
+    assert_approx_eq!(f64::consts::E.ln(), 1f64);
+    assert_eq!(1f32.ln(), 0.0);
+    assert_eq!(1f64.ln(), 0.0);
 
     assert_approx_eq!(0f32.ln_1p(), 0f32);
     assert_approx_eq!(0f64.ln_1p(), 0f64);
@@ -1047,7 +1139,8 @@ pub fn libm() {
 
     // Trigonometric functions.
 
-    assert_approx_eq!(0f32.sin(), 0f32);
+    assert_eq!(0f32.sin(), 0f32);
+    assert_eq!(0f64.sin(), 0f64);
     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);
@@ -1059,7 +1152,23 @@ pub fn libm() {
     assert_approx_eq!(2.0f32.asinh(), 1.443635475178810342493276740273105f32);
     assert_approx_eq!((-2.0f64).asinh(), -1.443635475178810342493276740273105f64);
 
-    assert_approx_eq!(0f32.cos(), 1f32);
+    // Ensure `sin` always returns something that is a valid input for `asin`, and same for
+    // `cos` and `acos`.
+    let halve_pi_f32 = std::f32::consts::FRAC_PI_2;
+    let halve_pi_f64 = std::f64::consts::FRAC_PI_2;
+    let pi_f32 = std::f32::consts::PI;
+    let pi_f64 = std::f64::consts::PI;
+    for _ in 0..64 {
+        // sin() should be clamped to [-1, 1] so asin() can never return NaN
+        assert!(!halve_pi_f32.sin().asin().is_nan());
+        assert!(!halve_pi_f64.sin().asin().is_nan());
+        // cos() should be clamped to [-1, 1] so acos() can never return NaN
+        assert!(!pi_f32.cos().acos().is_nan());
+        assert!(!pi_f64.cos().acos().is_nan());
+    }
+
+    assert_eq!(0f32.cos(), 1f32);
+    assert_eq!(0f64.cos(), 1f64);
     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);