miri/
math.rs

1use rand::Rng as _;
2use rustc_apfloat::Float as _;
3use rustc_apfloat::ieee::IeeeFloat;
4use rustc_middle::ty::{self, FloatTy, ScalarInt};
5
6use crate::*;
7
8/// Disturbes a floating-point result by a relative error in the range (-2^scale, 2^scale).
9///
10/// For a 2^N ULP error, you can use an `err_scale` of `-(F::PRECISION - 1 - N)`.
11/// In other words, a 1 ULP (absolute) error is the same as a `2^-(F::PRECISION-1)` relative error.
12/// (Subtracting 1 compensates for the integer bit.)
13pub(crate) fn apply_random_float_error<F: rustc_apfloat::Float>(
14    ecx: &mut crate::MiriInterpCx<'_>,
15    val: F,
16    err_scale: i32,
17) -> F {
18    let rng = ecx.machine.rng.get_mut();
19    // Generate a random integer in the range [0, 2^PREC).
20    // (When read as binary, the position of the first `1` determines the exponent,
21    // and the remaining bits fill the mantissa. `PREC` is one plus the size of the mantissa,
22    // so this all works out.)
23    let r = F::from_u128(rng.random_range(0..(1 << F::PRECISION))).value;
24    // Multiply this with 2^(scale - PREC). The result is between 0 and
25    // 2^PREC * 2^(scale - PREC) = 2^scale.
26    let err = r.scalbn(err_scale.strict_sub(F::PRECISION.try_into().unwrap()));
27    // give it a random sign
28    let err = if rng.random() { -err } else { err };
29    // multiple the value with (1+err)
30    (val * (F::from_u128(1).value + err).value).value
31}
32
33/// [`apply_random_float_error`] gives instructions to apply a 2^N ULP error.
34/// This function implements these instructions such that applying a 2^N ULP error is less error prone.
35/// So for a 2^N ULP error, you would pass N as the `ulp_exponent` argument.
36pub(crate) fn apply_random_float_error_ulp<F: rustc_apfloat::Float>(
37    ecx: &mut crate::MiriInterpCx<'_>,
38    val: F,
39    ulp_exponent: u32,
40) -> F {
41    let n = i32::try_from(ulp_exponent)
42        .expect("`err_scale_for_ulp`: exponent is too large to create an error scale");
43    // we know this fits
44    let prec = i32::try_from(F::PRECISION).unwrap();
45    let err_scale = -(prec - n - 1);
46    apply_random_float_error(ecx, val, err_scale)
47}
48
49/// Applies a random 16ULP floating point error to `val` and returns the new value.
50/// Will fail if `val` is not a floating point number.
51pub(crate) fn apply_random_float_error_to_imm<'tcx>(
52    ecx: &mut MiriInterpCx<'tcx>,
53    val: ImmTy<'tcx>,
54    ulp_exponent: u32,
55) -> InterpResult<'tcx, ImmTy<'tcx>> {
56    let scalar = val.to_scalar_int()?;
57    let res: ScalarInt = match val.layout.ty.kind() {
58        ty::Float(FloatTy::F16) =>
59            apply_random_float_error_ulp(ecx, scalar.to_f16(), ulp_exponent).into(),
60        ty::Float(FloatTy::F32) =>
61            apply_random_float_error_ulp(ecx, scalar.to_f32(), ulp_exponent).into(),
62        ty::Float(FloatTy::F64) =>
63            apply_random_float_error_ulp(ecx, scalar.to_f64(), ulp_exponent).into(),
64        ty::Float(FloatTy::F128) =>
65            apply_random_float_error_ulp(ecx, scalar.to_f128(), ulp_exponent).into(),
66        _ => bug!("intrinsic called with non-float input type"),
67    };
68
69    interp_ok(ImmTy::from_scalar_int(res, val.layout))
70}
71
72pub(crate) fn sqrt<S: rustc_apfloat::ieee::Semantics>(x: IeeeFloat<S>) -> IeeeFloat<S> {
73    match x.category() {
74        // preserve zero sign
75        rustc_apfloat::Category::Zero => x,
76        // propagate NaN
77        rustc_apfloat::Category::NaN => x,
78        // sqrt of negative number is NaN
79        _ if x.is_negative() => IeeeFloat::NAN,
80        // sqrt(∞) = ∞
81        rustc_apfloat::Category::Infinity => IeeeFloat::INFINITY,
82        rustc_apfloat::Category::Normal => {
83            // Floating point precision, excluding the integer bit
84            let prec = i32::try_from(S::PRECISION).unwrap() - 1;
85
86            // x = 2^(exp - prec) * mant
87            // where mant is an integer with prec+1 bits
88            // mant is a u128, which should be large enough for the largest prec (112 for f128)
89            let mut exp = x.ilogb();
90            let mut mant = x.scalbn(prec - exp).to_u128(128).value;
91
92            if exp % 2 != 0 {
93                // Make exponent even, so it can be divided by 2
94                exp -= 1;
95                mant <<= 1;
96            }
97
98            // Bit-by-bit (base-2 digit-by-digit) sqrt of mant.
99            // mant is treated here as a fixed point number with prec fractional bits.
100            // mant will be shifted left by one bit to have an extra fractional bit, which
101            // will be used to determine the rounding direction.
102
103            // res is the truncated sqrt of mant, where one bit is added at each iteration.
104            let mut res = 0u128;
105            // rem is the remainder with the current res
106            // rem_i = 2^i * ((mant<<1) - res_i^2)
107            // starting with res = 0, rem = mant<<1
108            let mut rem = mant << 1;
109            // s_i = 2*res_i
110            let mut s = 0u128;
111            // d is used to iterate over bits, from high to low (d_i = 2^(-i))
112            let mut d = 1u128 << (prec + 1);
113
114            // For iteration j=i+1, we need to find largest b_j = 0 or 1 such that
115            //  (res_i + b_j * 2^(-j))^2 <= mant<<1
116            // Expanding (a + b)^2 = a^2 + b^2 + 2*a*b:
117            //  res_i^2 + (b_j * 2^(-j))^2 + 2 * res_i * b_j * 2^(-j) <= mant<<1
118            // And rearranging the terms:
119            //  b_j^2 * 2^(-j) + 2 * res_i * b_j <= 2^j * (mant<<1 - res_i^2)
120            //  b_j^2 * 2^(-j) + 2 * res_i * b_j <= rem_i
121
122            while d != 0 {
123                // Probe b_j^2 * 2^(-j) + 2 * res_i * b_j <= rem_i with b_j = 1:
124                // t = 2*res_i + 2^(-j)
125                let t = s + d;
126                if rem >= t {
127                    // b_j should be 1, so make res_j = res_i + 2^(-j) and adjust rem
128                    res += d;
129                    s += d + d;
130                    rem -= t;
131                }
132                // Adjust rem for next iteration
133                rem <<= 1;
134                // Shift iterator
135                d >>= 1;
136            }
137
138            // Remove extra fractional bit from result, rounding to nearest.
139            // If the last bit is 0, then the nearest neighbor is definitely the lower one.
140            // If the last bit is 1, it sounds like this may either be a tie (if there's
141            // infinitely many 0s after this 1), or the nearest neighbor is the upper one.
142            // However, since square roots are either exact or irrational, and an exact root
143            // would lead to the last "extra" bit being 0, we can exclude a tie in this case.
144            // We therefore always round up if the last bit is 1. When the last bit is 0,
145            // adding 1 will not do anything since the shift will discard it.
146            res = (res + 1) >> 1;
147
148            // Build resulting value with res as mantissa and exp/2 as exponent
149            IeeeFloat::from_u128(res).value.scalbn(exp / 2 - prec)
150        }
151    }
152}
153
154#[cfg(test)]
155mod tests {
156    use rustc_apfloat::ieee::{DoubleS, HalfS, IeeeFloat, QuadS, SingleS};
157
158    use super::sqrt;
159
160    #[test]
161    fn test_sqrt() {
162        #[track_caller]
163        fn test<S: rustc_apfloat::ieee::Semantics>(x: &str, expected: &str) {
164            let x: IeeeFloat<S> = x.parse().unwrap();
165            let expected: IeeeFloat<S> = expected.parse().unwrap();
166            let result = sqrt(x);
167            assert_eq!(result, expected);
168        }
169
170        fn exact_tests<S: rustc_apfloat::ieee::Semantics>() {
171            test::<S>("0", "0");
172            test::<S>("1", "1");
173            test::<S>("1.5625", "1.25");
174            test::<S>("2.25", "1.5");
175            test::<S>("4", "2");
176            test::<S>("5.0625", "2.25");
177            test::<S>("9", "3");
178            test::<S>("16", "4");
179            test::<S>("25", "5");
180            test::<S>("36", "6");
181            test::<S>("49", "7");
182            test::<S>("64", "8");
183            test::<S>("81", "9");
184            test::<S>("100", "10");
185
186            test::<S>("0.5625", "0.75");
187            test::<S>("0.25", "0.5");
188            test::<S>("0.0625", "0.25");
189            test::<S>("0.00390625", "0.0625");
190        }
191
192        exact_tests::<HalfS>();
193        exact_tests::<SingleS>();
194        exact_tests::<DoubleS>();
195        exact_tests::<QuadS>();
196
197        test::<SingleS>("2", "1.4142135");
198        test::<DoubleS>("2", "1.4142135623730951");
199
200        test::<SingleS>("1.1", "1.0488088");
201        test::<DoubleS>("1.1", "1.0488088481701516");
202
203        test::<SingleS>("2.2", "1.4832398");
204        test::<DoubleS>("2.2", "1.4832396974191326");
205
206        test::<SingleS>("1.22101e-40", "1.10499205e-20");
207        test::<DoubleS>("1.22101e-310", "1.1049932126488395e-155");
208
209        test::<SingleS>("3.4028235e38", "1.8446743e19");
210        test::<DoubleS>("1.7976931348623157e308", "1.3407807929942596e154");
211    }
212}