core/num/imp/dec2flt/
slow.rs

1//! Slow, fallback algorithm for cases the Eisel-Lemire algorithm cannot round.
2
3use dec2flt::common::BiasedFp;
4use dec2flt::decimal_seq::{DecimalSeq, parse_decimal_seq};
5use dec2flt::float::RawFloat;
6
7use crate::num::imp::dec2flt;
8
9/// Parse the significant digits and biased, binary exponent of a float.
10///
11/// This is a fallback algorithm that uses a big-integer representation
12/// of the float, and therefore is considerably slower than faster
13/// approximations. However, it will always determine how to round
14/// the significant digits to the nearest machine float, allowing
15/// use to handle near half-way cases.
16///
17/// Near half-way cases are halfway between two consecutive machine floats.
18/// For example, the float `16777217.0` has a bitwise representation of
19/// `100000000000000000000000 1`. Rounding to a single-precision float,
20/// the trailing `1` is truncated. Using round-nearest, tie-even, any
21/// value above `16777217.0` must be rounded up to `16777218.0`, while
22/// any value before or equal to `16777217.0` must be rounded down
23/// to `16777216.0`. These near-halfway conversions therefore may require
24/// a large number of digits to unambiguously determine how to round.
25///
26/// The algorithms described here are based on "Processing Long Numbers Quickly",
27/// available here: <https://arxiv.org/pdf/2101.11408.pdf#section.11>.
28pub(crate) fn parse_long_mantissa<F: RawFloat>(s: &[u8]) -> BiasedFp {
29    const MAX_SHIFT: usize = 60;
30    const NUM_POWERS: usize = 19;
31    const POWERS: [u8; 19] =
32        [0, 3, 6, 9, 13, 16, 19, 23, 26, 29, 33, 36, 39, 43, 46, 49, 53, 56, 59];
33
34    let get_shift = |n| {
35        if n < NUM_POWERS { POWERS[n] as usize } else { MAX_SHIFT }
36    };
37
38    let fp_zero = BiasedFp::zero_pow2(0);
39    let fp_inf = BiasedFp::zero_pow2(F::INFINITE_POWER);
40
41    let mut d = parse_decimal_seq(s);
42
43    // Short-circuit if the value can only be a literal 0 or infinity.
44    if d.num_digits == 0 || d.decimal_point < -324 {
45        return fp_zero;
46    } else if d.decimal_point >= 310 {
47        return fp_inf;
48    }
49    let mut exp2 = 0_i32;
50    // Shift right toward (1/2 ... 1].
51    while d.decimal_point > 0 {
52        let n = d.decimal_point as usize;
53        let shift = get_shift(n);
54        d.right_shift(shift);
55        if d.decimal_point < -DecimalSeq::DECIMAL_POINT_RANGE {
56            return fp_zero;
57        }
58        exp2 += shift as i32;
59    }
60    // Shift left toward (1/2 ... 1].
61    while d.decimal_point <= 0 {
62        let shift = if d.decimal_point == 0 {
63            match d.digits[0] {
64                digit if digit >= 5 => break,
65                0 | 1 => 2,
66                _ => 1,
67            }
68        } else {
69            get_shift((-d.decimal_point) as _)
70        };
71        d.left_shift(shift);
72        if d.decimal_point > DecimalSeq::DECIMAL_POINT_RANGE {
73            return fp_inf;
74        }
75        exp2 -= shift as i32;
76    }
77    // We are now in the range [1/2 ... 1] but the binary format uses [1 ... 2].
78    exp2 -= 1;
79    while F::EXP_MIN > exp2 {
80        let mut n = (F::EXP_MIN - exp2) as usize;
81        if n > MAX_SHIFT {
82            n = MAX_SHIFT;
83        }
84        d.right_shift(n);
85        exp2 += n as i32;
86    }
87    if (exp2 - F::EXP_MIN + 1) >= F::INFINITE_POWER {
88        return fp_inf;
89    }
90    // Shift the decimal to the hidden bit, and then round the value
91    // to get the high mantissa+1 bits.
92    d.left_shift(F::SIG_BITS as usize + 1);
93    let mut mantissa = d.round();
94    if mantissa >= (1_u64 << (F::SIG_BITS + 1)) {
95        // Rounding up overflowed to the carry bit, need to
96        // shift back to the hidden bit.
97        d.right_shift(1);
98        exp2 += 1;
99        mantissa = d.round();
100        if (exp2 - F::EXP_MIN + 1) >= F::INFINITE_POWER {
101            return fp_inf;
102        }
103    }
104    let mut power2 = exp2 - F::EXP_MIN + 1;
105    if mantissa < (1_u64 << F::SIG_BITS) {
106        power2 -= 1;
107    }
108    // Zero out all the bits above the explicit mantissa bits.
109    mantissa &= (1_u64 << F::SIG_BITS) - 1;
110    BiasedFp { m: mantissa, p_biased: power2 }
111}
core/num/imp/dec2flt/slow.rs

core/num/imp/dec2flt/
slow.rs
