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}