226 lines
8.9 KiB
Rust
226 lines
8.9 KiB
Rust
|
//! Implementation of the Eisel-Lemire algorithm.
|
||
|
|
//!
|
||
|
|
//! This is adapted from [fast-float-rust](https://github.com/aldanor/fast-float-rust),
|
||
|
|
//! a port of [fast_float](https://github.com/fastfloat/fast_float) to Rust.
|
||
|
|
|
||
|
|
#![cfg(not(feature = "compact"))]
|
||
|
|
#![doc(hidden)]
|
||
|
|
|
||
|
|
use crate::extended_float::ExtendedFloat;
|
||
|
|
use crate::num::Float;
|
||
|
|
use crate::number::Number;
|
||
|
|
use crate::table::{LARGEST_POWER_OF_FIVE, POWER_OF_FIVE_128, SMALLEST_POWER_OF_FIVE};
|
||
|
|
|
||
|
|
/// Ensure truncation of digits doesn't affect our computation, by doing 2 passes.
|
||
|
|
#[inline]
|
||
|
|
pub fn lemire<F: Float>(num: &Number) -> ExtendedFloat {
|
||
|
|
// If significant digits were truncated, then we can have rounding error
|
||
|
|
// only if `mantissa + 1` produces a different result. We also avoid
|
||
|
|
// redundantly using the Eisel-Lemire algorithm if it was unable to
|
||
|
|
// correctly round on the first pass.
|
||
|
|
let mut fp = compute_float::<F>(num.exponent, num.mantissa);
|
||
|
|
if num.many_digits && fp.exp >= 0 && fp != compute_float::<F>(num.exponent, num.mantissa + 1) {
|
||
|
|
// Need to re-calculate, since the previous values are rounded
|
||
|
|
// when the slow path algorithm expects a normalized extended float.
|
||
|
|
fp = compute_error::<F>(num.exponent, num.mantissa);
|
||
|
|
}
|
||
|
|
fp
|
||
|
|
}
|
||
|
|
|
||
|
|
/// Compute a float using an extended-precision representation.
|
||
|
|
///
|
||
|
|
/// Fast conversion of a the significant digits and decimal exponent
|
||
|
|
/// a float to a extended representation with a binary float. This
|
||
|
|
/// algorithm will accurately parse the vast majority of cases,
|
||
|
|
/// and uses a 128-bit representation (with a fallback 192-bit
|
||
|
|
/// representation).
|
||
|
|
///
|
||
|
|
/// This algorithm scales the exponent by the decimal exponent
|
||
|
|
/// using pre-computed powers-of-5, and calculates if the
|
||
|
|
/// representation can be unambiguously rounded to the nearest
|
||
|
|
/// machine float. Near-halfway cases are not handled here,
|
||
|
|
/// and are represented by a negative, biased binary exponent.
|
||
|
|
///
|
||
|
|
/// The algorithm is described in detail in "Daniel Lemire, Number Parsing
|
||
|
|
/// at a Gigabyte per Second" in section 5, "Fast Algorithm", and
|
||
|
|
/// section 6, "Exact Numbers And Ties", available online:
|
||
|
|
/// <https://arxiv.org/abs/2101.11408.pdf>.
|
||
|
|
pub fn compute_float<F: Float>(q: i32, mut w: u64) -> ExtendedFloat {
|
||
|
|
let fp_zero = ExtendedFloat {
|
||
|
|
mant: 0,
|
||
|
|
exp: 0,
|
||
|
|
};
|
||
|
|
let fp_inf = ExtendedFloat {
|
||
|
|
mant: 0,
|
||
|
|
exp: F::INFINITE_POWER,
|
||
|
|
};
|
||
|
|
|
||
|
|
// Short-circuit if the value can only be a literal 0 or infinity.
|
||
|
|
if w == 0 || q < F::SMALLEST_POWER_OF_TEN {
|
||
|
|
return fp_zero;
|
||
|
|
} else if q > F::LARGEST_POWER_OF_TEN {
|
||
|
|
return fp_inf;
|
||
|
|
}
|
||
|
|
// Normalize our significant digits, so the most-significant bit is set.
|
||
|
|
let lz = w.leading_zeros() as i32;
|
||
|
|
w <<= lz;
|
||
|
|
let (lo, hi) = compute_product_approx(q, w, F::MANTISSA_SIZE as usize + 3);
|
||
|
|
if lo == 0xFFFF_FFFF_FFFF_FFFF {
|
||
|
|
// If we have failed to approximate w x 5^-q with our 128-bit value.
|
||
|
|
// Since the addition of 1 could lead to an overflow which could then
|
||
|
|
// round up over the half-way point, this can lead to improper rounding
|
||
|
|
// of a float.
|
||
|
|
//
|
||
|
|
// However, this can only occur if q ∈ [-27, 55]. The upper bound of q
|
||
|
|
// is 55 because 5^55 < 2^128, however, this can only happen if 5^q > 2^64,
|
||
|
|
// since otherwise the product can be represented in 64-bits, producing
|
||
|
|
// an exact result. For negative exponents, rounding-to-even can
|
||
|
|
// only occur if 5^-q < 2^64.
|
||
|
|
//
|
||
|
|
// For detailed explanations of rounding for negative exponents, see
|
||
|
|
// <https://arxiv.org/pdf/2101.11408.pdf#section.9.1>. For detailed
|
||
|
|
// explanations of rounding for positive exponents, see
|
||
|
|
// <https://arxiv.org/pdf/2101.11408.pdf#section.8>.
|
||
|
|
let inside_safe_exponent = (q >= -27) && (q <= 55);
|
||
|
|
if !inside_safe_exponent {
|
||
|
|
return compute_error_scaled::<F>(q, hi, lz);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
let upperbit = (hi >> 63) as i32;
|
||
|
|
let mut mantissa = hi >> (upperbit + 64 - F::MANTISSA_SIZE - 3);
|
||
|
|
let mut power2 = power(q) + upperbit - lz - F::MINIMUM_EXPONENT;
|
||
|
|
if power2 <= 0 {
|
||
|
|
if -power2 + 1 >= 64 {
|
||
|
|
// Have more than 64 bits below the minimum exponent, must be 0.
|
||
|
|
return fp_zero;
|
||
|
|
}
|
||
|
|
// Have a subnormal value.
|
||
|
|
mantissa >>= -power2 + 1;
|
||
|
|
mantissa += mantissa & 1;
|
||
|
|
mantissa >>= 1;
|
||
|
|
power2 = (mantissa >= (1_u64 << F::MANTISSA_SIZE)) as i32;
|
||
|
|
return ExtendedFloat {
|
||
|
|
mant: mantissa,
|
||
|
|
exp: power2,
|
||
|
|
};
|
||
|
|
}
|
||
|
|
// Need to handle rounding ties. Normally, we need to round up,
|
||
|
|
// but if we fall right in between and and we have an even basis, we
|
||
|
|
// need to round down.
|
||
|
|
//
|
||
|
|
// This will only occur if:
|
||
|
|
// 1. The lower 64 bits of the 128-bit representation is 0.
|
||
|
|
// IE, 5^q fits in single 64-bit word.
|
||
|
|
// 2. The least-significant bit prior to truncated mantissa is odd.
|
||
|
|
// 3. All the bits truncated when shifting to mantissa bits + 1 are 0.
|
||
|
|
//
|
||
|
|
// Or, we may fall between two floats: we are exactly halfway.
|
||
|
|
if lo <= 1
|
||
|
|
&& q >= F::MIN_EXPONENT_ROUND_TO_EVEN
|
||
|
|
&& q <= F::MAX_EXPONENT_ROUND_TO_EVEN
|
||
|
|
&& mantissa & 3 == 1
|
||
|
|
&& (mantissa << (upperbit + 64 - F::MANTISSA_SIZE - 3)) == hi
|
||
|
|
{
|
||
|
|
// Zero the lowest bit, so we don't round up.
|
||
|
|
mantissa &= !1_u64;
|
||
|
|
}
|
||
|
|
// Round-to-even, then shift the significant digits into place.
|
||
|
|
mantissa += mantissa & 1;
|
||
|
|
mantissa >>= 1;
|
||
|
|
if mantissa >= (2_u64 << F::MANTISSA_SIZE) {
|
||
|
|
// Rounding up overflowed, so the carry bit is set. Set the
|
||
|
|
// mantissa to 1 (only the implicit, hidden bit is set) and
|
||
|
|
// increase the exponent.
|
||
|
|
mantissa = 1_u64 << F::MANTISSA_SIZE;
|
||
|
|
power2 += 1;
|
||
|
|
}
|
||
|
|
// Zero out the hidden bit.
|
||
|
|
mantissa &= !(1_u64 << F::MANTISSA_SIZE);
|
||
|
|
if power2 >= F::INFINITE_POWER {
|
||
|
|
// Exponent is above largest normal value, must be infinite.
|
||
|
|
return fp_inf;
|
||
|
|
}
|
||
|
|
ExtendedFloat {
|
||
|
|
mant: mantissa,
|
||
|
|
exp: power2,
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
/// Fallback algorithm to calculate the non-rounded representation.
|
||
|
|
/// This calculates the extended representation, and then normalizes
|
||
|
|
/// the resulting representation, so the high bit is set.
|
||
|
|
#[inline]
|
||
|
|
pub fn compute_error<F: Float>(q: i32, mut w: u64) -> ExtendedFloat {
|
||
|
|
let lz = w.leading_zeros() as i32;
|
||
|
|
w <<= lz;
|
||
|
|
let hi = compute_product_approx(q, w, F::MANTISSA_SIZE as usize + 3).1;
|
||
|
|
compute_error_scaled::<F>(q, hi, lz)
|
||
|
|
}
|
||
|
|
|
||
|
|
/// Compute the error from a mantissa scaled to the exponent.
|
||
|
|
#[inline]
|
||
|
|
pub fn compute_error_scaled<F: Float>(q: i32, mut w: u64, lz: i32) -> ExtendedFloat {
|
||
|
|
// Want to normalize the float, but this is faster than ctlz on most architectures.
|
||
|
|
let hilz = (w >> 63) as i32 ^ 1;
|
||
|
|
w <<= hilz;
|
||
|
|
let power2 = power(q as i32) + F::EXPONENT_BIAS - hilz - lz - 62;
|
||
|
|
|
||
|
|
ExtendedFloat {
|
||
|
|
mant: w,
|
||
|
|
exp: power2 + F::INVALID_FP,
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
/// Calculate a base 2 exponent from a decimal exponent.
|
||
|
|
/// This uses a pre-computed integer approximation for
|
||
|
|
/// log2(10), where 217706 / 2^16 is accurate for the
|
||
|
|
/// entire range of non-finite decimal exponents.
|
||
|
|
#[inline]
|
||
|
|
fn power(q: i32) -> i32 {
|
||
|
|
(q.wrapping_mul(152_170 + 65536) >> 16) + 63
|
||
|
|
}
|
||
|
|
|
||
|
|
#[inline]
|
||
|
|
fn full_multiplication(a: u64, b: u64) -> (u64, u64) {
|
||
|
|
let r = (a as u128) * (b as u128);
|
||
|
|
(r as u64, (r >> 64) as u64)
|
||
|
|
}
|
||
|
|
|
||
|
|
// This will compute or rather approximate w * 5**q and return a pair of 64-bit words
|
||
|
|
// approximating the result, with the "high" part corresponding to the most significant
|
||
|
|
// bits and the low part corresponding to the least significant bits.
|
||
|
|
fn compute_product_approx(q: i32, w: u64, precision: usize) -> (u64, u64) {
|
||
|
|
debug_assert!(q >= SMALLEST_POWER_OF_FIVE);
|
||
|
|
debug_assert!(q <= LARGEST_POWER_OF_FIVE);
|
||
|
|
debug_assert!(precision <= 64);
|
||
|
|
|
||
|
|
let mask = if precision < 64 {
|
||
|
|
0xFFFF_FFFF_FFFF_FFFF_u64 >> precision
|
||
|
|
} else {
|
||
|
|
0xFFFF_FFFF_FFFF_FFFF_u64
|
||
|
|
};
|
||
|
|
|
||
|
|
// 5^q < 2^64, then the multiplication always provides an exact value.
|
||
|
|
// That means whenever we need to round ties to even, we always have
|
||
|
|
// an exact value.
|
||
|
|
let index = (q - SMALLEST_POWER_OF_FIVE) as usize;
|
||
|
|
let (lo5, hi5) = POWER_OF_FIVE_128[index];
|
||
|
|
// Only need one multiplication as long as there is 1 zero but
|
||
|
|
// in the explicit mantissa bits, +1 for the hidden bit, +1 to
|
||
|
|
// determine the rounding direction, +1 for if the computed
|
||
|
|
// product has a leading zero.
|
||
|
|
let (mut first_lo, mut first_hi) = full_multiplication(w, lo5);
|
||
|
|
if first_hi & mask == mask {
|
||
|
|
// Need to do a second multiplication to get better precision
|
||
|
|
// for the lower product. This will always be exact
|
||
|
|
// where q is < 55, since 5^55 < 2^128. If this wraps,
|
||
|
|
// then we need to need to round up the hi product.
|
||
|
|
let (_, second_hi) = full_multiplication(w, hi5);
|
||
|
|
first_lo = first_lo.wrapping_add(second_hi);
|
||
|
|
if second_hi > first_lo {
|
||
|
|
first_hi += 1;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
(first_lo, first_hi)
|
||
|
|
}
|