mesh_loader/utils/float/
slow.rs

Help
1//! Slow, fallback algorithm for cases the Eisel-Lemire algorithm cannot round.
2
3use super::{
4    common::BiasedFp,
5    decimal::{parse_decimal, Decimal},
6    float::RawFloat,
7};
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    static 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
35    let get_shift = |n| {
36        if n < NUM_POWERS {
37            POWERS[n] as usize
38        } else {
39            MAX_SHIFT
40        }
41    };
42
43    let fp_zero = BiasedFp::zero_pow2(0);
44    let fp_inf = BiasedFp::zero_pow2(F::INFINITE_POWER);
45
46    let mut d = parse_decimal(s);
47
48    // Short-circuit if the value can only be a literal 0 or infinity.
49    if d.num_digits == 0 || d.decimal_point < -324 {
50        return fp_zero;
51    } else if d.decimal_point >= 310 {
52        return fp_inf;
53    }
54    let mut exp2 = 0_i32;
55    // Shift right toward (1/2 ... 1].
56    while d.decimal_point > 0 {
57        let n = d.decimal_point as usize;
58        let shift = get_shift(n);
59        d.right_shift(shift);
60        if d.decimal_point < -Decimal::DECIMAL_POINT_RANGE {
61            return fp_zero;
62        }
63        exp2 += shift as i32;
64    }
65    // Shift left toward (1/2 ... 1].
66    while d.decimal_point <= 0 {
67        let shift = if d.decimal_point == 0 {
68            match d.digits[0] {
69                digit if digit >= 5 => break,
70                0 | 1 => 2,
71                _ => 1,
72            }
73        } else {
74            get_shift((-d.decimal_point) as usize)
75        };
76        d.left_shift(shift);
77        if d.decimal_point > Decimal::DECIMAL_POINT_RANGE {
78            return fp_inf;
79        }
80        exp2 -= shift as i32;
81    }
82    // We are now in the range [1/2 ... 1] but the binary format uses [1 ... 2].
83    exp2 -= 1;
84    while (F::MINIMUM_EXPONENT + 1) > exp2 {
85        let mut n = ((F::MINIMUM_EXPONENT + 1) - exp2) as usize;
86        if n > MAX_SHIFT {
87            n = MAX_SHIFT;
88        }
89        d.right_shift(n);
90        exp2 += n as i32;
91    }
92    if (exp2 - F::MINIMUM_EXPONENT) >= F::INFINITE_POWER {
93        return fp_inf;
94    }
95    // Shift the decimal to the hidden bit, and then round the value
96    // to get the high mantissa+1 bits.
97    d.left_shift(F::MANTISSA_EXPLICIT_BITS + 1);
98    let mut mantissa = d.round();
99    if mantissa >= (1_u64 << (F::MANTISSA_EXPLICIT_BITS + 1)) {
100        // Rounding up overflowed to the carry bit, need to
101        // shift back to the hidden bit.
102        d.right_shift(1);
103        exp2 += 1;
104        mantissa = d.round();
105        if (exp2 - F::MINIMUM_EXPONENT) >= F::INFINITE_POWER {
106            return fp_inf;
107        }
108    }
109    let mut power2 = exp2 - F::MINIMUM_EXPONENT;
110    if mantissa < (1_u64 << F::MANTISSA_EXPLICIT_BITS) {
111        power2 -= 1;
112    }
113    // Zero out all the bits above the explicit mantissa bits.
114    mantissa &= (1_u64 << F::MANTISSA_EXPLICIT_BITS) - 1;
115    BiasedFp {
116        f: mantissa,
117        e: power2,
118    }
119}
mesh_loader/utils/float/slow.rs

mesh_loader/utils/float/
slow.rs
