|
| 1 | +use rand::Rng as _; |
| 2 | +use rand::distributions::Distribution as _; |
| 3 | +use rustc_apfloat::Float as _; |
| 4 | +use rustc_apfloat::ieee::IeeeFloat; |
| 5 | + |
| 6 | +/// Disturbes a floating-point result by a relative error on the order of (-2^scale, 2^scale). |
| 7 | +pub(crate) fn apply_random_float_error<F: rustc_apfloat::Float>( |
| 8 | + ecx: &mut crate::MiriInterpCx<'_>, |
| 9 | + val: F, |
| 10 | + err_scale: i32, |
| 11 | +) -> F { |
| 12 | + let rng = ecx.machine.rng.get_mut(); |
| 13 | + // Generate a random integer in the range [0, 2^PREC). |
| 14 | + let dist = rand::distributions::Uniform::new(0, 1 << F::PRECISION); |
| 15 | + let err = F::from_u128(dist.sample(rng)) |
| 16 | + .value |
| 17 | + .scalbn(err_scale.strict_sub(F::PRECISION.try_into().unwrap())); |
| 18 | + // give it a random sign |
| 19 | + let err = if rng.gen::<bool>() { -err } else { err }; |
| 20 | + // multiple the value with (1+err) |
| 21 | + (val * (F::from_u128(1).value + err).value).value |
| 22 | +} |
| 23 | + |
| 24 | +pub(crate) fn sqrt<S: rustc_apfloat::ieee::Semantics>(x: IeeeFloat<S>) -> IeeeFloat<S> { |
| 25 | + match x.category() { |
| 26 | + // preserve zero sign |
| 27 | + rustc_apfloat::Category::Zero => x, |
| 28 | + // propagate NaN |
| 29 | + rustc_apfloat::Category::NaN => x, |
| 30 | + // sqrt of negative number is NaN |
| 31 | + _ if x.is_negative() => IeeeFloat::NAN, |
| 32 | + // sqrt(∞) = ∞ |
| 33 | + rustc_apfloat::Category::Infinity => IeeeFloat::INFINITY, |
| 34 | + rustc_apfloat::Category::Normal => { |
| 35 | + // Floating point precision, excluding the integer bit |
| 36 | + let prec = i32::try_from(S::PRECISION).unwrap() - 1; |
| 37 | + |
| 38 | + // x = 2^(exp - prec) * mant |
| 39 | + // where mant is an integer with prec+1 bits |
| 40 | + // mant is a u128, which should be large enough for the largest prec (112 for f128) |
| 41 | + let mut exp = x.ilogb(); |
| 42 | + let mut mant = x.scalbn(prec - exp).to_u128(128).value; |
| 43 | + |
| 44 | + if exp % 2 != 0 { |
| 45 | + // Make exponent even, so it can be divided by 2 |
| 46 | + exp -= 1; |
| 47 | + mant <<= 1; |
| 48 | + } |
| 49 | + |
| 50 | + // Bit-by-bit (base-2 digit-by-digit) sqrt of mant. |
| 51 | + // mant is treated here as a fixed point number with prec fractional bits. |
| 52 | + // mant will be shifted left by one bit to have an extra fractional bit, which |
| 53 | + // will be used to determine the rounding direction. |
| 54 | + |
| 55 | + // res is the truncated sqrt of mant, where one bit is added at each iteration. |
| 56 | + let mut res = 0u128; |
| 57 | + // rem is the remainder with the current res |
| 58 | + // rem_i = 2^i * ((mant<<1) - res_i^2) |
| 59 | + // starting with res = 0, rem = mant<<1 |
| 60 | + let mut rem = mant << 1; |
| 61 | + // s_i = 2*res_i |
| 62 | + let mut s = 0u128; |
| 63 | + // d is used to iterate over bits, from high to low (d_i = 2^(-i)) |
| 64 | + let mut d = 1u128 << (prec + 1); |
| 65 | + |
| 66 | + // For iteration j=i+1, we need to find largest b_j = 0 or 1 such that |
| 67 | + // (res_i + b_j * 2^(-j))^2 <= mant<<1 |
| 68 | + // Expanding (a + b)^2 = a^2 + b^2 + 2*a*b: |
| 69 | + // res_i^2 + (b_j * 2^(-j))^2 + 2 * res_i * b_j * 2^(-j) <= mant<<1 |
| 70 | + // And rearranging the terms: |
| 71 | + // b_j^2 * 2^(-j) + 2 * res_i * b_j <= 2^j * (mant<<1 - res_i^2) |
| 72 | + // b_j^2 * 2^(-j) + 2 * res_i * b_j <= rem_i |
| 73 | + |
| 74 | + while d != 0 { |
| 75 | + // Probe b_j^2 * 2^(-j) + 2 * res_i * b_j <= rem_i with b_j = 1: |
| 76 | + // t = 2*res_i + 2^(-j) |
| 77 | + let t = s + d; |
| 78 | + if rem >= t { |
| 79 | + // b_j should be 1, so make res_j = res_i + 2^(-j) and adjust rem |
| 80 | + res += d; |
| 81 | + s += d + d; |
| 82 | + rem -= t; |
| 83 | + } |
| 84 | + // Adjust rem for next iteration |
| 85 | + rem <<= 1; |
| 86 | + // Shift iterator |
| 87 | + d >>= 1; |
| 88 | + } |
| 89 | + |
| 90 | + // Remove extra fractional bit from result, rounding to nearest. |
| 91 | + // If the last bit is 0, then the nearest neighbor is definitely the lower one. |
| 92 | + // If the last bit is 1, it sounds like this may either be a tie (if there's |
| 93 | + // infinitely many 0s after this 1), or the nearest neighbor is the upper one. |
| 94 | + // However, since square roots are either exact or irrational, and an exact root |
| 95 | + // would lead to the last "extra" bit being 0, we can exclude a tie in this case. |
| 96 | + // We therefore always round up if the last bit is 1. When the last bit is 0, |
| 97 | + // adding 1 will not do anything since the shift will discard it. |
| 98 | + res = (res + 1) >> 1; |
| 99 | + |
| 100 | + // Build resulting value with res as mantissa and exp/2 as exponent |
| 101 | + IeeeFloat::from_u128(res).value.scalbn(exp / 2 - prec) |
| 102 | + } |
| 103 | + } |
| 104 | +} |
| 105 | + |
| 106 | +#[cfg(test)] |
| 107 | +mod tests { |
| 108 | + use rustc_apfloat::ieee::{DoubleS, HalfS, IeeeFloat, QuadS, SingleS}; |
| 109 | + |
| 110 | + use super::sqrt; |
| 111 | + |
| 112 | + #[test] |
| 113 | + fn test_sqrt() { |
| 114 | + #[track_caller] |
| 115 | + fn test<S: rustc_apfloat::ieee::Semantics>(x: &str, expected: &str) { |
| 116 | + let x: IeeeFloat<S> = x.parse().unwrap(); |
| 117 | + let expected: IeeeFloat<S> = expected.parse().unwrap(); |
| 118 | + let result = sqrt(x); |
| 119 | + assert_eq!(result, expected); |
| 120 | + } |
| 121 | + |
| 122 | + fn exact_tests<S: rustc_apfloat::ieee::Semantics>() { |
| 123 | + test::<S>("0", "0"); |
| 124 | + test::<S>("1", "1"); |
| 125 | + test::<S>("1.5625", "1.25"); |
| 126 | + test::<S>("2.25", "1.5"); |
| 127 | + test::<S>("4", "2"); |
| 128 | + test::<S>("5.0625", "2.25"); |
| 129 | + test::<S>("9", "3"); |
| 130 | + test::<S>("16", "4"); |
| 131 | + test::<S>("25", "5"); |
| 132 | + test::<S>("36", "6"); |
| 133 | + test::<S>("49", "7"); |
| 134 | + test::<S>("64", "8"); |
| 135 | + test::<S>("81", "9"); |
| 136 | + test::<S>("100", "10"); |
| 137 | + |
| 138 | + test::<S>("0.5625", "0.75"); |
| 139 | + test::<S>("0.25", "0.5"); |
| 140 | + test::<S>("0.0625", "0.25"); |
| 141 | + test::<S>("0.00390625", "0.0625"); |
| 142 | + } |
| 143 | + |
| 144 | + exact_tests::<HalfS>(); |
| 145 | + exact_tests::<SingleS>(); |
| 146 | + exact_tests::<DoubleS>(); |
| 147 | + exact_tests::<QuadS>(); |
| 148 | + |
| 149 | + test::<SingleS>("2", "1.4142135"); |
| 150 | + test::<DoubleS>("2", "1.4142135623730951"); |
| 151 | + |
| 152 | + test::<SingleS>("1.1", "1.0488088"); |
| 153 | + test::<DoubleS>("1.1", "1.0488088481701516"); |
| 154 | + |
| 155 | + test::<SingleS>("2.2", "1.4832398"); |
| 156 | + test::<DoubleS>("2.2", "1.4832396974191326"); |
| 157 | + |
| 158 | + test::<SingleS>("1.22101e-40", "1.10499205e-20"); |
| 159 | + test::<DoubleS>("1.22101e-310", "1.1049932126488395e-155"); |
| 160 | + |
| 161 | + test::<SingleS>("3.4028235e38", "1.8446743e19"); |
| 162 | + test::<DoubleS>("1.7976931348623157e308", "1.3407807929942596e154"); |
| 163 | + } |
| 164 | +} |
0 commit comments