Adding upstream version 115.7.0esr.upstream/115.7.0esrupstream

Signed-off-by: Daniel Baumann <daniel.baumann@progress-linux.org>

1 files changed, 391 insertions, 0 deletions

diff --git a/third_party/rust/minimal-lexical/src/bellerophon.rs b/third_party/rust/minimal-lexical/src/bellerophon.rs
new file mode 100644
index 0000000000..86a2023d09
--- /dev/null
+++ b/third_party/rust/minimal-lexical/src/bellerophon.rs
@@ -0,0 +1,391 @@
+//! An implementation of Clinger's Bellerophon algorithm.
+//!
+//! This is a moderate path algorithm that uses an extended-precision
+//! float, represented in 80 bits, by calculating the bits of slop
+//! and determining if those bits could prevent unambiguous rounding.
+//!
+//! This algorithm requires less static storage than the Lemire algorithm,
+//! and has decent performance, and is therefore used when non-decimal,
+//! non-power-of-two strings need to be parsed. Clinger's algorithm
+//! is described in depth in "How to Read Floating Point Numbers Accurately.",
+//! available online [here](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.45.4152&rep=rep1&type=pdf).
+//!
+//! This implementation is loosely based off the Golang implementation,
+//! found [here](https://github.com/golang/go/blob/b10849fbb97a2244c086991b4623ae9f32c212d0/src/strconv/extfloat.go).
+//! This code is therefore subject to a 3-clause BSD license.
+
+#![cfg(feature = "compact")]
+#![doc(hidden)]
+
+use crate::extended_float::ExtendedFloat;
+use crate::mask::{lower_n_halfway, lower_n_mask};
+use crate::num::Float;
+use crate::number::Number;
+use crate::rounding::{round, round_nearest_tie_even};
+use crate::table::BASE10_POWERS;
+
+// ALGORITHM
+// ---------
+
+/// Core implementation of the Bellerophon algorithm.
+///
+/// Create an extended-precision float, scale it to the proper radix power,
+/// calculate the bits of slop, and return the representation. The value
+/// will always be guaranteed to be within 1 bit, rounded-down, of the real
+/// value. If a negative exponent is returned, this represents we were
+/// unable to unambiguously round the significant digits.
+///
+/// This has been modified to return a biased, rather than unbiased exponent.
+pub fn bellerophon<F: Float>(num: &Number) -> ExtendedFloat {
+    let fp_zero = ExtendedFloat {
+        mant: 0,
+        exp: 0,
+    };
+    let fp_inf = ExtendedFloat {
+        mant: 0,
+        exp: F::INFINITE_POWER,
+    };
+
+    // Early short-circuit, in case of literal 0 or infinity.
+    // This allows us to avoid narrow casts causing numeric overflow,
+    // and is a quick check for any radix.
+    if num.mantissa == 0 || num.exponent <= -0x1000 {
+        return fp_zero;
+    } else if num.exponent >= 0x1000 {
+        return fp_inf;
+    }
+
+    // Calculate our indexes for our extended-precision multiplication.
+    // This narrowing cast is safe, since exponent must be in a valid range.
+    let exponent = num.exponent as i32 + BASE10_POWERS.bias;
+    let small_index = exponent % BASE10_POWERS.step;
+    let large_index = exponent / BASE10_POWERS.step;
+
+    if exponent < 0 {
+        // Guaranteed underflow (assign 0).
+        return fp_zero;
+    }
+    if large_index as usize >= BASE10_POWERS.large.len() {
+        // Overflow (assign infinity)
+        return fp_inf;
+    }
+
+    // Within the valid exponent range, multiply by the large and small
+    // exponents and return the resulting value.
+
+    // Track errors to as a factor of unit in last-precision.
+    let mut errors: u32 = 0;
+    if num.many_digits {
+        errors += error_halfscale();
+    }
+
+    // Multiply by the small power.
+    // Check if we can directly multiply by an integer, if not,
+    // use extended-precision multiplication.
+    let mut fp = ExtendedFloat {
+        mant: num.mantissa,
+        exp: 0,
+    };
+    match fp.mant.overflowing_mul(BASE10_POWERS.get_small_int(small_index as usize)) {
+        // Overflow, multiplication unsuccessful, go slow path.
+        (_, true) => {
+            normalize(&mut fp);
+            fp = mul(&fp, &BASE10_POWERS.get_small(small_index as usize));
+            errors += error_halfscale();
+        },
+        // No overflow, multiplication successful.
+        (mant, false) => {
+            fp.mant = mant;
+            normalize(&mut fp);
+        },
+    }
+
+    // Multiply by the large power.
+    fp = mul(&fp, &BASE10_POWERS.get_large(large_index as usize));
+    if errors > 0 {
+        errors += 1;
+    }
+    errors += error_halfscale();
+
+    // Normalize the floating point (and the errors).
+    let shift = normalize(&mut fp);
+    errors <<= shift;
+    fp.exp += F::EXPONENT_BIAS;
+
+    // Check for literal overflow, even with halfway cases.
+    if -fp.exp + 1 > 65 {
+        return fp_zero;
+    }
+
+    // Too many errors accumulated, return an error.
+    if !error_is_accurate::<F>(errors, &fp) {
+        // Bias the exponent so we know it's invalid.
+        fp.exp += F::INVALID_FP;
+        return fp;
+    }
+
+    // Check if we have a literal 0 or overflow here.
+    // If we have an exponent of -63, we can still have a valid shift,
+    // giving a case where we have too many errors and need to round-up.
+    if -fp.exp + 1 == 65 {
+        // Have more than 64 bits below the minimum exponent, must be 0.
+        return fp_zero;
+    }
+
+    round::<F, _>(&mut fp, |f, s| {
+        round_nearest_tie_even(f, s, |is_odd, is_halfway, is_above| {
+            is_above || (is_odd && is_halfway)
+        });
+    });
+    fp
+}
+
+// ERRORS
+// ------
+
+// Calculate if the errors in calculating the extended-precision float.
+//
+// Specifically, we want to know if we are close to a halfway representation,
+// or halfway between `b` and `b+1`, or `b+h`. The halfway representation
+// has the form:
+//     SEEEEEEEHMMMMMMMMMMMMMMMMMMMMMMM100...
+// where:
+//     S = Sign Bit
+//     E = Exponent Bits
+//     H = Hidden Bit
+//     M = Mantissa Bits
+//
+// The halfway representation has a bit set 1-after the mantissa digits,
+// and no bits set immediately afterward, making it impossible to
+// round between `b` and `b+1` with this representation.
+
+/// Get the full error scale.
+#[inline(always)]
+const fn error_scale() -> u32 {
+    8
+}
+
+/// Get the half error scale.
+#[inline(always)]
+const fn error_halfscale() -> u32 {
+    error_scale() / 2
+}
+
+/// Determine if the number of errors is tolerable for float precision.
+fn error_is_accurate<F: Float>(errors: u32, fp: &ExtendedFloat) -> bool {
+    // Check we can't have a literal 0 denormal float.
+    debug_assert!(fp.exp >= -64);
+
+    // Determine if extended-precision float is a good approximation.
+    // If the error has affected too many units, the float will be
+    // inaccurate, or if the representation is too close to halfway
+    // that any operations could affect this halfway representation.
+    // See the documentation for dtoa for more information.
+
+    // This is always a valid u32, since `fp.exp >= -64`
+    // will always be positive and the significand size is {23, 52}.
+    let mantissa_shift = 64 - F::MANTISSA_SIZE - 1;
+
+    // The unbiased exponent checks is `unbiased_exp <= F::MANTISSA_SIZE
+    // - F::EXPONENT_BIAS -64 + 1`, or `biased_exp <= F::MANTISSA_SIZE - 63`,
+    // or `biased_exp <= mantissa_shift`.
+    let extrabits = match fp.exp <= -mantissa_shift {
+        // Denormal, since shifting to the hidden bit still has a negative exponent.
+        // The unbiased check calculation for bits is `1 - F::EXPONENT_BIAS - unbiased_exp`,
+        // or `1 - biased_exp`.
+        true => 1 - fp.exp,
+        false => 64 - F::MANTISSA_SIZE - 1,
+    };
+
+    // Our logic is as follows: we want to determine if the actual
+    // mantissa and the errors during calculation differ significantly
+    // from the rounding point. The rounding point for round-nearest
+    // is the halfway point, IE, this when the truncated bits start
+    // with b1000..., while the rounding point for the round-toward
+    // is when the truncated bits are equal to 0.
+    // To do so, we can check whether the rounding point +/- the error
+    // are >/< the actual lower n bits.
+    //
+    // For whether we need to use signed or unsigned types for this
+    // analysis, see this example, using u8 rather than u64 to simplify
+    // things.
+    //
+    // # Comparisons
+    //      cmp1 = (halfway - errors) < extra
+    //      cmp1 = extra < (halfway + errors)
+    //
+    // # Large Extrabits, Low Errors
+    //
+    //      extrabits = 8
+    //      halfway          =  0b10000000
+    //      extra            =  0b10000010
+    //      errors           =  0b00000100
+    //      halfway - errors =  0b01111100
+    //      halfway + errors =  0b10000100
+    //
+    //      Unsigned:
+    //          halfway - errors = 124
+    //          halfway + errors = 132
+    //          extra            = 130
+    //          cmp1             = true
+    //          cmp2             = true
+    //      Signed:
+    //          halfway - errors = 124
+    //          halfway + errors = -124
+    //          extra            = -126
+    //          cmp1             = false
+    //          cmp2             = true
+    //
+    // # Conclusion
+    //
+    // Since errors will always be small, and since we want to detect
+    // if the representation is accurate, we need to use an **unsigned**
+    // type for comparisons.
+    let maskbits = extrabits as u64;
+    let errors = errors as u64;
+
+    // Round-to-nearest, need to use the halfway point.
+    if extrabits > 64 {
+        // Underflow, we have a shift larger than the mantissa.
+        // Representation is valid **only** if the value is close enough
+        // overflow to the next bit within errors. If it overflows,
+        // the representation is **not** valid.
+        !fp.mant.overflowing_add(errors).1
+    } else {
+        let mask = lower_n_mask(maskbits);
+        let extra = fp.mant & mask;
+
+        // Round-to-nearest, need to check if we're close to halfway.
+        // IE, b10100 | 100000, where `|` signifies the truncation point.
+        let halfway = lower_n_halfway(maskbits);
+        let cmp1 = halfway.wrapping_sub(errors) < extra;
+        let cmp2 = extra < halfway.wrapping_add(errors);
+
+        // If both comparisons are true, we have significant rounding error,
+        // and the value cannot be exactly represented. Otherwise, the
+        // representation is valid.
+        !(cmp1 && cmp2)
+    }
+}
+
+// MATH
+// ----
+
+/// Normalize float-point number.
+///
+/// Shift the mantissa so the number of leading zeros is 0, or the value
+/// itself is 0.
+///
+/// Get the number of bytes shifted.
+pub fn normalize(fp: &mut ExtendedFloat) -> i32 {
+    // Note:
+    // Using the ctlz intrinsic via leading_zeros is way faster (~10x)
+    // than shifting 1-bit at a time, via while loop, and also way
+    // faster (~2x) than an unrolled loop that checks at 32, 16, 4,
+    // 2, and 1 bit.
+    //
+    // Using a modulus of pow2 (which will get optimized to a bitwise
+    // and with 0x3F or faster) is slightly slower than an if/then,
+    // however, removing the if/then will likely optimize more branched
+    // code as it removes conditional logic.
+
+    // Calculate the number of leading zeros, and then zero-out
+    // any overflowing bits, to avoid shl overflow when self.mant == 0.
+    if fp.mant != 0 {
+        let shift = fp.mant.leading_zeros() as i32;
+        fp.mant <<= shift;
+        fp.exp -= shift;
+        shift
+    } else {
+        0
+    }
+}
+
+/// Multiply two normalized extended-precision floats, as if by `a*b`.
+///
+/// The precision is maximal when the numbers are normalized, however,
+/// decent precision will occur as long as both values have high bits
+/// set. The result is not normalized.
+///
+/// Algorithm:
+///     1. Non-signed multiplication of mantissas (requires 2x as many bits as input).
+///     2. Normalization of the result (not done here).
+///     3. Addition of exponents.
+pub fn mul(x: &ExtendedFloat, y: &ExtendedFloat) -> ExtendedFloat {
+    // Logic check, values must be decently normalized prior to multiplication.
+    debug_assert!(x.mant >> 32 != 0);
+    debug_assert!(y.mant >> 32 != 0);
+
+    // Extract high-and-low masks.
+    // Mask is u32::MAX for older Rustc versions.
+    const LOMASK: u64 = 0xffff_ffff;
+    let x1 = x.mant >> 32;
+    let x0 = x.mant & LOMASK;
+    let y1 = y.mant >> 32;
+    let y0 = y.mant & LOMASK;
+
+    // Get our products
+    let x1_y0 = x1 * y0;
+    let x0_y1 = x0 * y1;
+    let x0_y0 = x0 * y0;
+    let x1_y1 = x1 * y1;
+
+    let mut tmp = (x1_y0 & LOMASK) + (x0_y1 & LOMASK) + (x0_y0 >> 32);
+    // round up
+    tmp += 1 << (32 - 1);
+
+    ExtendedFloat {
+        mant: x1_y1 + (x1_y0 >> 32) + (x0_y1 >> 32) + (tmp >> 32),
+        exp: x.exp + y.exp + 64,
+    }
+}
+
+// POWERS
+// ------
+
+/// Precalculated powers of base N for the Bellerophon algorithm.
+pub struct BellerophonPowers {
+    // Pre-calculated small powers.
+    pub small: &'static [u64],
+    // Pre-calculated large powers.
+    pub large: &'static [u64],
+    /// Pre-calculated small powers as 64-bit integers
+    pub small_int: &'static [u64],
+    // Step between large powers and number of small powers.
+    pub step: i32,
+    // Exponent bias for the large powers.
+    pub bias: i32,
+    /// ceil(log2(radix)) scaled as a multiplier.
+    pub log2: i64,
+    /// Bitshift for the log2 multiplier.
+    pub log2_shift: i32,
+}
+
+/// Allow indexing of values without bounds checking
+impl BellerophonPowers {
+    #[inline]
+    pub fn get_small(&self, index: usize) -> ExtendedFloat {
+        let mant = self.small[index];
+        let exp = (1 - 64) + ((self.log2 * index as i64) >> self.log2_shift);
+        ExtendedFloat {
+            mant,
+            exp: exp as i32,
+        }
+    }
+
+    #[inline]
+    pub fn get_large(&self, index: usize) -> ExtendedFloat {
+        let mant = self.large[index];
+        let biased_e = index as i64 * self.step as i64 - self.bias as i64;
+        let exp = (1 - 64) + ((self.log2 * biased_e) >> self.log2_shift);
+        ExtendedFloat {
+            mant,
+            exp: exp as i32,
+        }
+    }
+
+    #[inline]
+    pub fn get_small_int(&self, index: usize) -> u64 {
+        self.small_int[index]
+    }
+}