Rollup merge of #134063 - tgross35:dec2flt-refactoring, r=Noratrieb

dec2flt: Clean up float parsing modules

This is the first portion of my work adding support for parsing and printing `f16`. Changes in `float.rs` replace the magic constants with expressions and add some use of generics to better support the new float types. Everything else is related to documentation or naming; there are no functional changes in this PR.

This can be reviewed by commit.
This commit is contained in:
许杰友 Jieyou Xu (Joe) 2025-03-05 21:46:31 +08:00 committed by GitHub
commit 9b8accbeb6
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
18 changed files with 851 additions and 600 deletions

View File

@ -8,12 +8,12 @@ pub(crate) trait ByteSlice {
/// Writes a 64-bit integer as 8 bytes in little-endian order.
fn write_u64(&mut self, value: u64);
/// Calculate the offset of a slice from another.
/// Calculate the difference in length between two slices.
fn offset_from(&self, other: &Self) -> isize;
/// Iteratively parse and consume digits from bytes.
/// Returns the same bytes with consumed digits being
/// elided.
///
/// Returns the same bytes with consumed digits being elided. Breaks on invalid digits.
fn parse_digits(&self, func: impl FnMut(u8)) -> &Self;
}
@ -39,11 +39,11 @@ impl ByteSlice for [u8] {
fn parse_digits(&self, mut func: impl FnMut(u8)) -> &Self {
let mut s = self;
while let Some((c, s_next)) = s.split_first() {
while let Some((c, rest)) = s.split_first() {
let c = c.wrapping_sub(b'0');
if c < 10 {
func(c);
s = s_next;
s = rest;
} else {
break;
}
@ -53,7 +53,9 @@ impl ByteSlice for [u8] {
}
}
/// Determine if 8 bytes are all decimal digits.
/// Determine if all characters in an 8-byte byte string (represented as a `u64`) are all decimal
/// digits.
///
/// This does not care about the order in which the bytes were loaded.
pub(crate) fn is_8digits(v: u64) -> bool {
let a = v.wrapping_add(0x4646_4646_4646_4646);
@ -61,19 +63,20 @@ pub(crate) fn is_8digits(v: u64) -> bool {
(a | b) & 0x8080_8080_8080_8080 == 0
}
/// A custom 64-bit floating point type, representing `f * 2^e`.
/// e is biased, so it be directly shifted into the exponent bits.
/// A custom 64-bit floating point type, representing `m * 2^p`.
/// p is biased, so it be directly shifted into the exponent bits.
#[derive(Debug, Copy, Clone, PartialEq, Eq, Default)]
pub struct BiasedFp {
/// The significant digits.
pub f: u64,
pub m: u64,
/// The biased, binary exponent.
pub e: i32,
pub p_biased: i32,
}
impl BiasedFp {
/// Represent `0 ^ p`
#[inline]
pub const fn zero_pow2(e: i32) -> Self {
Self { f: 0, e }
pub const fn zero_pow2(p_biased: i32) -> Self {
Self { m: 0, p_biased }
}
}

View File

@ -1,358 +1,87 @@
//! Arbitrary-precision decimal class for fallback algorithms.
//!
//! This is only used if the fast-path (native floats) and
//! the Eisel-Lemire algorithm are unable to unambiguously
//! determine the float.
//!
//! The technique used is "Simple Decimal Conversion", developed
//! by Nigel Tao and Ken Thompson. A detailed description of the
//! algorithm can be found in "ParseNumberF64 by Simple Decimal Conversion",
//! available online: <https://nigeltao.github.io/blog/2020/parse-number-f64-simple.html>.
//! Representation of a float as the significant digits and exponent.
use crate::num::dec2flt::common::{ByteSlice, is_8digits};
use crate::num::dec2flt::float::RawFloat;
use crate::num::dec2flt::fpu::set_precision;
#[derive(Clone)]
pub(super) struct Decimal {
/// The number of significant digits in the decimal.
pub num_digits: usize,
/// The offset of the decimal point in the significant digits.
pub decimal_point: i32,
/// If the number of significant digits stored in the decimal is truncated.
pub truncated: bool,
/// Buffer of the raw digits, in the range [0, 9].
pub digits: [u8; Self::MAX_DIGITS],
}
const INT_POW10: [u64; 16] = [
1,
10,
100,
1000,
10000,
100000,
1000000,
10000000,
100000000,
1000000000,
10000000000,
100000000000,
1000000000000,
10000000000000,
100000000000000,
1000000000000000,
];
impl Default for Decimal {
fn default() -> Self {
Self { num_digits: 0, decimal_point: 0, truncated: false, digits: [0; Self::MAX_DIGITS] }
}
/// A floating point number with up to 64 bits of mantissa and an `i64` exponent.
#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
pub struct Decimal {
pub exponent: i64,
pub mantissa: u64,
pub negative: bool,
pub many_digits: bool,
}
impl Decimal {
/// The maximum number of digits required to unambiguously round a float.
///
/// For a double-precision IEEE 754 float, this required 767 digits,
/// so we store the max digits + 1.
///
/// We can exactly represent a float in radix `b` from radix 2 if
/// `b` is divisible by 2. This function calculates the exact number of
/// digits required to exactly represent that float.
///
/// According to the "Handbook of Floating Point Arithmetic",
/// for IEEE754, with emin being the min exponent, p2 being the
/// precision, and b being the radix, the number of digits follows as:
///
/// `emin + p2 + ⌊(emin + 1) log(2, b) log(1 2^(p2), b)⌋`
///
/// For f32, this follows as:
/// emin = -126
/// p2 = 24
///
/// For f64, this follows as:
/// emin = -1022
/// p2 = 53
///
/// In Python:
/// `-emin + p2 + math.floor((emin+ 1)*math.log(2, b)-math.log(1-2**(-p2), b))`
pub(super) const MAX_DIGITS: usize = 768;
/// The max digits that can be exactly represented in a 64-bit integer.
pub(super) const MAX_DIGITS_WITHOUT_OVERFLOW: usize = 19;
pub(super) const DECIMAL_POINT_RANGE: i32 = 2047;
/// Append a digit to the buffer.
pub(super) fn try_add_digit(&mut self, digit: u8) {
if self.num_digits < Self::MAX_DIGITS {
self.digits[self.num_digits] = digit;
}
self.num_digits += 1;
/// Detect if the float can be accurately reconstructed from native floats.
#[inline]
fn can_use_fast_path<F: RawFloat>(&self) -> bool {
F::MIN_EXPONENT_FAST_PATH <= self.exponent
&& self.exponent <= F::MAX_EXPONENT_DISGUISED_FAST_PATH
&& self.mantissa <= F::MAX_MANTISSA_FAST_PATH
&& !self.many_digits
}
/// Trim trailing zeros from the buffer.
pub(super) fn trim(&mut self) {
// All of the following calls to `Decimal::trim` can't panic because:
//
// 1. `parse_decimal` sets `num_digits` to a max of `Decimal::MAX_DIGITS`.
// 2. `right_shift` sets `num_digits` to `write_index`, which is bounded by `num_digits`.
// 3. `left_shift` `num_digits` to a max of `Decimal::MAX_DIGITS`.
//
// Trim is only called in `right_shift` and `left_shift`.
debug_assert!(self.num_digits <= Self::MAX_DIGITS);
while self.num_digits != 0 && self.digits[self.num_digits - 1] == 0 {
self.num_digits -= 1;
}
}
/// Try turning the decimal into an exact float representation, using machine-sized integers
/// and floats.
///
/// This is extracted into a separate function so that it can be attempted before constructing
/// a Decimal. This only works if both the mantissa and the exponent
/// can be exactly represented as a machine float, since IEE-754 guarantees
/// no rounding will occur.
///
/// There is an exception: disguised fast-path cases, where we can shift
/// powers-of-10 from the exponent to the significant digits.
pub fn try_fast_path<F: RawFloat>(&self) -> Option<F> {
// Here we need to work around <https://github.com/rust-lang/rust/issues/114479>.
// The fast path crucially depends on arithmetic being rounded to the correct number of bits
// without any intermediate rounding. On x86 (without SSE or SSE2) this requires the precision
// of the x87 FPU stack to be changed so that it directly rounds to 64/32 bit.
// The `set_precision` function takes care of setting the precision on architectures which
// require setting it by changing the global state (like the control word of the x87 FPU).
let _cw = set_precision::<F>();
pub(super) fn round(&self) -> u64 {
if self.num_digits == 0 || self.decimal_point < 0 {
return 0;
} else if self.decimal_point > 18 {
return 0xFFFF_FFFF_FFFF_FFFF_u64;
if !self.can_use_fast_path::<F>() {
return None;
}
let dp = self.decimal_point as usize;
let mut n = 0_u64;
for i in 0..dp {
n *= 10;
if i < self.num_digits {
n += self.digits[i] as u64;
}
}
let mut round_up = false;
if dp < self.num_digits {
round_up = self.digits[dp] >= 5;
if self.digits[dp] == 5 && dp + 1 == self.num_digits {
round_up = self.truncated || ((dp != 0) && (1 & self.digits[dp - 1] != 0))
}
}
if round_up {
n += 1;
}
n
}
/// Computes decimal * 2^shift.
pub(super) fn left_shift(&mut self, shift: usize) {
if self.num_digits == 0 {
return;
}
let num_new_digits = number_of_digits_decimal_left_shift(self, shift);
let mut read_index = self.num_digits;
let mut write_index = self.num_digits + num_new_digits;
let mut n = 0_u64;
while read_index != 0 {
read_index -= 1;
write_index -= 1;
n += (self.digits[read_index] as u64) << shift;
let quotient = n / 10;
let remainder = n - (10 * quotient);
if write_index < Self::MAX_DIGITS {
self.digits[write_index] = remainder as u8;
} else if remainder > 0 {
self.truncated = true;
}
n = quotient;
}
while n > 0 {
write_index -= 1;
let quotient = n / 10;
let remainder = n - (10 * quotient);
if write_index < Self::MAX_DIGITS {
self.digits[write_index] = remainder as u8;
} else if remainder > 0 {
self.truncated = true;
}
n = quotient;
}
self.num_digits += num_new_digits;
if self.num_digits > Self::MAX_DIGITS {
self.num_digits = Self::MAX_DIGITS;
}
self.decimal_point += num_new_digits as i32;
self.trim();
}
/// Computes decimal * 2^-shift.
pub(super) fn right_shift(&mut self, shift: usize) {
let mut read_index = 0;
let mut write_index = 0;
let mut n = 0_u64;
while (n >> shift) == 0 {
if read_index < self.num_digits {
n = (10 * n) + self.digits[read_index] as u64;
read_index += 1;
} else if n == 0 {
return;
let value = if self.exponent <= F::MAX_EXPONENT_FAST_PATH {
// normal fast path
let value = F::from_u64(self.mantissa);
if self.exponent < 0 {
value / F::pow10_fast_path((-self.exponent) as _)
} else {
while (n >> shift) == 0 {
n *= 10;
read_index += 1;
}
break;
value * F::pow10_fast_path(self.exponent as _)
}
}
self.decimal_point -= read_index as i32 - 1;
if self.decimal_point < -Self::DECIMAL_POINT_RANGE {
// `self = Self::Default()`, but without the overhead of clearing `digits`.
self.num_digits = 0;
self.decimal_point = 0;
self.truncated = false;
return;
}
let mask = (1_u64 << shift) - 1;
while read_index < self.num_digits {
let new_digit = (n >> shift) as u8;
n = (10 * (n & mask)) + self.digits[read_index] as u64;
read_index += 1;
self.digits[write_index] = new_digit;
write_index += 1;
}
while n > 0 {
let new_digit = (n >> shift) as u8;
n = 10 * (n & mask);
if write_index < Self::MAX_DIGITS {
self.digits[write_index] = new_digit;
write_index += 1;
} else if new_digit > 0 {
self.truncated = true;
}
}
self.num_digits = write_index;
self.trim();
}
}
/// Parse a big integer representation of the float as a decimal.
pub(super) fn parse_decimal(mut s: &[u8]) -> Decimal {
let mut d = Decimal::default();
let start = s;
while let Some((&b'0', s_next)) = s.split_first() {
s = s_next;
}
s = s.parse_digits(|digit| d.try_add_digit(digit));
if let Some((b'.', s_next)) = s.split_first() {
s = s_next;
let first = s;
// Skip leading zeros.
if d.num_digits == 0 {
while let Some((&b'0', s_next)) = s.split_first() {
s = s_next;
}
}
while s.len() >= 8 && d.num_digits + 8 < Decimal::MAX_DIGITS {
let v = s.read_u64();
if !is_8digits(v) {
break;
}
d.digits[d.num_digits..].write_u64(v - 0x3030_3030_3030_3030);
d.num_digits += 8;
s = &s[8..];
}
s = s.parse_digits(|digit| d.try_add_digit(digit));
d.decimal_point = s.len() as i32 - first.len() as i32;
}
if d.num_digits != 0 {
// Ignore the trailing zeros if there are any
let mut n_trailing_zeros = 0;
for &c in start[..(start.len() - s.len())].iter().rev() {
if c == b'0' {
n_trailing_zeros += 1;
} else if c != b'.' {
break;
}
}
d.decimal_point += n_trailing_zeros as i32;
d.num_digits -= n_trailing_zeros;
d.decimal_point += d.num_digits as i32;
if d.num_digits > Decimal::MAX_DIGITS {
d.truncated = true;
d.num_digits = Decimal::MAX_DIGITS;
}
}
if let Some((&ch, s_next)) = s.split_first() {
if ch == b'e' || ch == b'E' {
s = s_next;
let mut neg_exp = false;
if let Some((&ch, s_next)) = s.split_first() {
neg_exp = ch == b'-';
if ch == b'-' || ch == b'+' {
s = s_next;
}
}
let mut exp_num = 0_i32;
s.parse_digits(|digit| {
if exp_num < 0x10000 {
exp_num = 10 * exp_num + digit as i32;
}
});
d.decimal_point += if neg_exp { -exp_num } else { exp_num };
}
}
for i in d.num_digits..Decimal::MAX_DIGITS_WITHOUT_OVERFLOW {
d.digits[i] = 0;
}
d
}
fn number_of_digits_decimal_left_shift(d: &Decimal, mut shift: usize) -> usize {
#[rustfmt::skip]
const TABLE: [u16; 65] = [
0x0000, 0x0800, 0x0801, 0x0803, 0x1006, 0x1009, 0x100D, 0x1812, 0x1817, 0x181D, 0x2024,
0x202B, 0x2033, 0x203C, 0x2846, 0x2850, 0x285B, 0x3067, 0x3073, 0x3080, 0x388E, 0x389C,
0x38AB, 0x38BB, 0x40CC, 0x40DD, 0x40EF, 0x4902, 0x4915, 0x4929, 0x513E, 0x5153, 0x5169,
0x5180, 0x5998, 0x59B0, 0x59C9, 0x61E3, 0x61FD, 0x6218, 0x6A34, 0x6A50, 0x6A6D, 0x6A8B,
0x72AA, 0x72C9, 0x72E9, 0x7B0A, 0x7B2B, 0x7B4D, 0x8370, 0x8393, 0x83B7, 0x83DC, 0x8C02,
0x8C28, 0x8C4F, 0x9477, 0x949F, 0x94C8, 0x9CF2, 0x051C, 0x051C, 0x051C, 0x051C,
];
#[rustfmt::skip]
const TABLE_POW5: [u8; 0x051C] = [
5, 2, 5, 1, 2, 5, 6, 2, 5, 3, 1, 2, 5, 1, 5, 6, 2, 5, 7, 8, 1, 2, 5, 3, 9, 0, 6, 2, 5, 1,
9, 5, 3, 1, 2, 5, 9, 7, 6, 5, 6, 2, 5, 4, 8, 8, 2, 8, 1, 2, 5, 2, 4, 4, 1, 4, 0, 6, 2, 5,
1, 2, 2, 0, 7, 0, 3, 1, 2, 5, 6, 1, 0, 3, 5, 1, 5, 6, 2, 5, 3, 0, 5, 1, 7, 5, 7, 8, 1, 2,
5, 1, 5, 2, 5, 8, 7, 8, 9, 0, 6, 2, 5, 7, 6, 2, 9, 3, 9, 4, 5, 3, 1, 2, 5, 3, 8, 1, 4, 6,
9, 7, 2, 6, 5, 6, 2, 5, 1, 9, 0, 7, 3, 4, 8, 6, 3, 2, 8, 1, 2, 5, 9, 5, 3, 6, 7, 4, 3, 1,
6, 4, 0, 6, 2, 5, 4, 7, 6, 8, 3, 7, 1, 5, 8, 2, 0, 3, 1, 2, 5, 2, 3, 8, 4, 1, 8, 5, 7, 9,
1, 0, 1, 5, 6, 2, 5, 1, 1, 9, 2, 0, 9, 2, 8, 9, 5, 5, 0, 7, 8, 1, 2, 5, 5, 9, 6, 0, 4, 6,
4, 4, 7, 7, 5, 3, 9, 0, 6, 2, 5, 2, 9, 8, 0, 2, 3, 2, 2, 3, 8, 7, 6, 9, 5, 3, 1, 2, 5, 1,
4, 9, 0, 1, 1, 6, 1, 1, 9, 3, 8, 4, 7, 6, 5, 6, 2, 5, 7, 4, 5, 0, 5, 8, 0, 5, 9, 6, 9, 2,
3, 8, 2, 8, 1, 2, 5, 3, 7, 2, 5, 2, 9, 0, 2, 9, 8, 4, 6, 1, 9, 1, 4, 0, 6, 2, 5, 1, 8, 6,
2, 6, 4, 5, 1, 4, 9, 2, 3, 0, 9, 5, 7, 0, 3, 1, 2, 5, 9, 3, 1, 3, 2, 2, 5, 7, 4, 6, 1, 5,
4, 7, 8, 5, 1, 5, 6, 2, 5, 4, 6, 5, 6, 6, 1, 2, 8, 7, 3, 0, 7, 7, 3, 9, 2, 5, 7, 8, 1, 2,
5, 2, 3, 2, 8, 3, 0, 6, 4, 3, 6, 5, 3, 8, 6, 9, 6, 2, 8, 9, 0, 6, 2, 5, 1, 1, 6, 4, 1, 5,
3, 2, 1, 8, 2, 6, 9, 3, 4, 8, 1, 4, 4, 5, 3, 1, 2, 5, 5, 8, 2, 0, 7, 6, 6, 0, 9, 1, 3, 4,
6, 7, 4, 0, 7, 2, 2, 6, 5, 6, 2, 5, 2, 9, 1, 0, 3, 8, 3, 0, 4, 5, 6, 7, 3, 3, 7, 0, 3, 6,
1, 3, 2, 8, 1, 2, 5, 1, 4, 5, 5, 1, 9, 1, 5, 2, 2, 8, 3, 6, 6, 8, 5, 1, 8, 0, 6, 6, 4, 0,
6, 2, 5, 7, 2, 7, 5, 9, 5, 7, 6, 1, 4, 1, 8, 3, 4, 2, 5, 9, 0, 3, 3, 2, 0, 3, 1, 2, 5, 3,
6, 3, 7, 9, 7, 8, 8, 0, 7, 0, 9, 1, 7, 1, 2, 9, 5, 1, 6, 6, 0, 1, 5, 6, 2, 5, 1, 8, 1, 8,
9, 8, 9, 4, 0, 3, 5, 4, 5, 8, 5, 6, 4, 7, 5, 8, 3, 0, 0, 7, 8, 1, 2, 5, 9, 0, 9, 4, 9, 4,
7, 0, 1, 7, 7, 2, 9, 2, 8, 2, 3, 7, 9, 1, 5, 0, 3, 9, 0, 6, 2, 5, 4, 5, 4, 7, 4, 7, 3, 5,
0, 8, 8, 6, 4, 6, 4, 1, 1, 8, 9, 5, 7, 5, 1, 9, 5, 3, 1, 2, 5, 2, 2, 7, 3, 7, 3, 6, 7, 5,
4, 4, 3, 2, 3, 2, 0, 5, 9, 4, 7, 8, 7, 5, 9, 7, 6, 5, 6, 2, 5, 1, 1, 3, 6, 8, 6, 8, 3, 7,
7, 2, 1, 6, 1, 6, 0, 2, 9, 7, 3, 9, 3, 7, 9, 8, 8, 2, 8, 1, 2, 5, 5, 6, 8, 4, 3, 4, 1, 8,
8, 6, 0, 8, 0, 8, 0, 1, 4, 8, 6, 9, 6, 8, 9, 9, 4, 1, 4, 0, 6, 2, 5, 2, 8, 4, 2, 1, 7, 0,
9, 4, 3, 0, 4, 0, 4, 0, 0, 7, 4, 3, 4, 8, 4, 4, 9, 7, 0, 7, 0, 3, 1, 2, 5, 1, 4, 2, 1, 0,
8, 5, 4, 7, 1, 5, 2, 0, 2, 0, 0, 3, 7, 1, 7, 4, 2, 2, 4, 8, 5, 3, 5, 1, 5, 6, 2, 5, 7, 1,
0, 5, 4, 2, 7, 3, 5, 7, 6, 0, 1, 0, 0, 1, 8, 5, 8, 7, 1, 1, 2, 4, 2, 6, 7, 5, 7, 8, 1, 2,
5, 3, 5, 5, 2, 7, 1, 3, 6, 7, 8, 8, 0, 0, 5, 0, 0, 9, 2, 9, 3, 5, 5, 6, 2, 1, 3, 3, 7, 8,
9, 0, 6, 2, 5, 1, 7, 7, 6, 3, 5, 6, 8, 3, 9, 4, 0, 0, 2, 5, 0, 4, 6, 4, 6, 7, 7, 8, 1, 0,
6, 6, 8, 9, 4, 5, 3, 1, 2, 5, 8, 8, 8, 1, 7, 8, 4, 1, 9, 7, 0, 0, 1, 2, 5, 2, 3, 2, 3, 3,
8, 9, 0, 5, 3, 3, 4, 4, 7, 2, 6, 5, 6, 2, 5, 4, 4, 4, 0, 8, 9, 2, 0, 9, 8, 5, 0, 0, 6, 2,
6, 1, 6, 1, 6, 9, 4, 5, 2, 6, 6, 7, 2, 3, 6, 3, 2, 8, 1, 2, 5, 2, 2, 2, 0, 4, 4, 6, 0, 4,
9, 2, 5, 0, 3, 1, 3, 0, 8, 0, 8, 4, 7, 2, 6, 3, 3, 3, 6, 1, 8, 1, 6, 4, 0, 6, 2, 5, 1, 1,
1, 0, 2, 2, 3, 0, 2, 4, 6, 2, 5, 1, 5, 6, 5, 4, 0, 4, 2, 3, 6, 3, 1, 6, 6, 8, 0, 9, 0, 8,
2, 0, 3, 1, 2, 5, 5, 5, 5, 1, 1, 1, 5, 1, 2, 3, 1, 2, 5, 7, 8, 2, 7, 0, 2, 1, 1, 8, 1, 5,
8, 3, 4, 0, 4, 5, 4, 1, 0, 1, 5, 6, 2, 5, 2, 7, 7, 5, 5, 5, 7, 5, 6, 1, 5, 6, 2, 8, 9, 1,
3, 5, 1, 0, 5, 9, 0, 7, 9, 1, 7, 0, 2, 2, 7, 0, 5, 0, 7, 8, 1, 2, 5, 1, 3, 8, 7, 7, 7, 8,
7, 8, 0, 7, 8, 1, 4, 4, 5, 6, 7, 5, 5, 2, 9, 5, 3, 9, 5, 8, 5, 1, 1, 3, 5, 2, 5, 3, 9, 0,
6, 2, 5, 6, 9, 3, 8, 8, 9, 3, 9, 0, 3, 9, 0, 7, 2, 2, 8, 3, 7, 7, 6, 4, 7, 6, 9, 7, 9, 2,
5, 5, 6, 7, 6, 2, 6, 9, 5, 3, 1, 2, 5, 3, 4, 6, 9, 4, 4, 6, 9, 5, 1, 9, 5, 3, 6, 1, 4, 1,
8, 8, 8, 2, 3, 8, 4, 8, 9, 6, 2, 7, 8, 3, 8, 1, 3, 4, 7, 6, 5, 6, 2, 5, 1, 7, 3, 4, 7, 2,
3, 4, 7, 5, 9, 7, 6, 8, 0, 7, 0, 9, 4, 4, 1, 1, 9, 2, 4, 4, 8, 1, 3, 9, 1, 9, 0, 6, 7, 3,
8, 2, 8, 1, 2, 5, 8, 6, 7, 3, 6, 1, 7, 3, 7, 9, 8, 8, 4, 0, 3, 5, 4, 7, 2, 0, 5, 9, 6, 2,
2, 4, 0, 6, 9, 5, 9, 5, 3, 3, 6, 9, 1, 4, 0, 6, 2, 5,
];
shift &= 63;
let x_a = TABLE[shift];
let x_b = TABLE[shift + 1];
let num_new_digits = (x_a >> 11) as _;
let pow5_a = (0x7FF & x_a) as usize;
let pow5_b = (0x7FF & x_b) as usize;
let pow5 = &TABLE_POW5[pow5_a..];
for (i, &p5) in pow5.iter().enumerate().take(pow5_b - pow5_a) {
if i >= d.num_digits {
return num_new_digits - 1;
} else if d.digits[i] == p5 {
continue;
} else if d.digits[i] < p5 {
return num_new_digits - 1;
} else {
return num_new_digits;
}
// disguised fast path
let shift = self.exponent - F::MAX_EXPONENT_FAST_PATH;
let mantissa = self.mantissa.checked_mul(INT_POW10[shift as usize])?;
if mantissa > F::MAX_MANTISSA_FAST_PATH {
return None;
}
F::from_u64(mantissa) * F::pow10_fast_path(F::MAX_EXPONENT_FAST_PATH as _)
};
if self.negative { Some(-value) } else { Some(value) }
}
num_new_digits
}

View File

@ -0,0 +1,379 @@
//! Arbitrary-precision decimal type used by fallback algorithms.
//!
//! This is only used if the fast-path (native floats) and
//! the Eisel-Lemire algorithm are unable to unambiguously
//! determine the float.
//!
//! The technique used is "Simple Decimal Conversion", developed
//! by Nigel Tao and Ken Thompson. A detailed description of the
//! algorithm can be found in "ParseNumberF64 by Simple Decimal Conversion",
//! available online: <https://nigeltao.github.io/blog/2020/parse-number-f64-simple.html>.
use crate::num::dec2flt::common::{ByteSlice, is_8digits};
/// A decimal floating-point number, represented as a sequence of decimal digits.
#[derive(Clone, Debug, PartialEq)]
pub struct DecimalSeq {
/// The number of significant digits in the decimal.
pub num_digits: usize,
/// The offset of the decimal point in the significant digits.
pub decimal_point: i32,
/// If the number of significant digits stored in the decimal is truncated.
pub truncated: bool,
/// Buffer of the raw digits, in the range [0, 9].
pub digits: [u8; Self::MAX_DIGITS],
}
impl Default for DecimalSeq {
fn default() -> Self {
Self { num_digits: 0, decimal_point: 0, truncated: false, digits: [0; Self::MAX_DIGITS] }
}
}
impl DecimalSeq {
/// The maximum number of digits required to unambiguously round up to a 64-bit float.
///
/// For an IEEE 754 binary64 float, this required 767 digits. So we store the max digits + 1.
///
/// We can exactly represent a float in radix `b` from radix 2 if
/// `b` is divisible by 2. This function calculates the exact number of
/// digits required to exactly represent that float.
///
/// According to the "Handbook of Floating Point Arithmetic",
/// for IEEE754, with `emin` being the min exponent, `p2` being the
/// precision, and `b` being the radix, the number of digits follows as:
///
/// `emin + p2 + ⌊(emin + 1) log(2, b) log(1 2^(p2), b)⌋`
///
/// For f32, this follows as:
/// emin = -126
/// p2 = 24
///
/// For f64, this follows as:
/// emin = -1022
/// p2 = 53
///
/// In Python:
/// `-emin + p2 + math.floor((emin+ 1)*math.log(2, b)-math.log(1-2**(-p2), b))`
pub const MAX_DIGITS: usize = 768;
/// The max decimal digits that can be exactly represented in a 64-bit integer.
pub(super) const MAX_DIGITS_WITHOUT_OVERFLOW: usize = 19;
pub(super) const DECIMAL_POINT_RANGE: i32 = 2047;
/// Append a digit to the buffer if it fits.
// FIXME(tgross35): it may be better for this to return an option
// FIXME(tgross35): incrementing the digit counter even if we don't push anything
// seems incorrect.
pub(super) fn try_add_digit(&mut self, digit: u8) {
if self.num_digits < Self::MAX_DIGITS {
self.digits[self.num_digits] = digit;
}
self.num_digits += 1;
}
/// Trim trailing zeros from the buffer.
// FIXME(tgross35): this could be `.rev().position()` if perf is okay
pub fn trim(&mut self) {
// All of the following calls to `DecimalSeq::trim` can't panic because:
//
// 1. `parse_decimal` sets `num_digits` to a max of `DecimalSeq::MAX_DIGITS`.
// 2. `right_shift` sets `num_digits` to `write_index`, which is bounded by `num_digits`.
// 3. `left_shift` `num_digits` to a max of `DecimalSeq::MAX_DIGITS`.
//
// Trim is only called in `right_shift` and `left_shift`.
debug_assert!(self.num_digits <= Self::MAX_DIGITS);
while self.num_digits != 0 && self.digits[self.num_digits - 1] == 0 {
self.num_digits -= 1;
}
}
pub(super) fn round(&self) -> u64 {
if self.num_digits == 0 || self.decimal_point < 0 {
return 0;
} else if self.decimal_point >= Self::MAX_DIGITS_WITHOUT_OVERFLOW as i32 {
return 0xFFFF_FFFF_FFFF_FFFF_u64;
}
let dp = self.decimal_point as usize;
let mut n = 0_u64;
for i in 0..dp {
n *= 10;
if i < self.num_digits {
n += self.digits[i] as u64;
}
}
let mut round_up = false;
if dp < self.num_digits {
round_up = self.digits[dp] >= 5;
if self.digits[dp] == 5 && dp + 1 == self.num_digits {
round_up = self.truncated || ((dp != 0) && (1 & self.digits[dp - 1] != 0))
}
}
if round_up {
n += 1;
}
n
}
/// Computes decimal * 2^shift.
pub(super) fn left_shift(&mut self, shift: usize) {
if self.num_digits == 0 {
return;
}
let num_new_digits = number_of_digits_decimal_left_shift(self, shift);
let mut read_index = self.num_digits;
let mut write_index = self.num_digits + num_new_digits;
let mut n = 0_u64;
while read_index != 0 {
read_index -= 1;
write_index -= 1;
n += (self.digits[read_index] as u64) << shift;
let quotient = n / 10;
let remainder = n - (10 * quotient);
if write_index < Self::MAX_DIGITS {
self.digits[write_index] = remainder as u8;
} else if remainder > 0 {
self.truncated = true;
}
n = quotient;
}
while n > 0 {
write_index -= 1;
let quotient = n / 10;
let remainder = n - (10 * quotient);
if write_index < Self::MAX_DIGITS {
self.digits[write_index] = remainder as u8;
} else if remainder > 0 {
self.truncated = true;
}
n = quotient;
}
self.num_digits += num_new_digits;
if self.num_digits > Self::MAX_DIGITS {
self.num_digits = Self::MAX_DIGITS;
}
self.decimal_point += num_new_digits as i32;
self.trim();
}
/// Computes decimal * 2^-shift.
pub(super) fn right_shift(&mut self, shift: usize) {
let mut read_index = 0;
let mut write_index = 0;
let mut n = 0_u64;
while (n >> shift) == 0 {
if read_index < self.num_digits {
n = (10 * n) + self.digits[read_index] as u64;
read_index += 1;
} else if n == 0 {
return;
} else {
while (n >> shift) == 0 {
n *= 10;
read_index += 1;
}
break;
}
}
self.decimal_point -= read_index as i32 - 1;
if self.decimal_point < -Self::DECIMAL_POINT_RANGE {
// `self = Self::Default()`, but without the overhead of clearing `digits`.
self.num_digits = 0;
self.decimal_point = 0;
self.truncated = false;
return;
}
let mask = (1_u64 << shift) - 1;
while read_index < self.num_digits {
let new_digit = (n >> shift) as u8;
n = (10 * (n & mask)) + self.digits[read_index] as u64;
read_index += 1;
self.digits[write_index] = new_digit;
write_index += 1;
}
while n > 0 {
let new_digit = (n >> shift) as u8;
n = 10 * (n & mask);
if write_index < Self::MAX_DIGITS {
self.digits[write_index] = new_digit;
write_index += 1;
} else if new_digit > 0 {
self.truncated = true;
}
}
self.num_digits = write_index;
self.trim();
}
}
/// Parse a big integer representation of the float as a decimal.
pub fn parse_decimal_seq(mut s: &[u8]) -> DecimalSeq {
let mut d = DecimalSeq::default();
let start = s;
while let Some((&b'0', s_next)) = s.split_first() {
s = s_next;
}
s = s.parse_digits(|digit| d.try_add_digit(digit));
if let Some((b'.', s_next)) = s.split_first() {
s = s_next;
let first = s;
// Skip leading zeros.
if d.num_digits == 0 {
while let Some((&b'0', s_next)) = s.split_first() {
s = s_next;
}
}
while s.len() >= 8 && d.num_digits + 8 < DecimalSeq::MAX_DIGITS {
let v = s.read_u64();
if !is_8digits(v) {
break;
}
d.digits[d.num_digits..].write_u64(v - 0x3030_3030_3030_3030);
d.num_digits += 8;
s = &s[8..];
}
s = s.parse_digits(|digit| d.try_add_digit(digit));
d.decimal_point = s.len() as i32 - first.len() as i32;
}
if d.num_digits != 0 {
// Ignore the trailing zeros if there are any
let mut n_trailing_zeros = 0;
for &c in start[..(start.len() - s.len())].iter().rev() {
if c == b'0' {
n_trailing_zeros += 1;
} else if c != b'.' {
break;
}
}
d.decimal_point += n_trailing_zeros as i32;
d.num_digits -= n_trailing_zeros;
d.decimal_point += d.num_digits as i32;
if d.num_digits > DecimalSeq::MAX_DIGITS {
d.truncated = true;
d.num_digits = DecimalSeq::MAX_DIGITS;
}
}
if let Some((&ch, s_next)) = s.split_first() {
if ch == b'e' || ch == b'E' {
s = s_next;
let mut neg_exp = false;
if let Some((&ch, s_next)) = s.split_first() {
neg_exp = ch == b'-';
if ch == b'-' || ch == b'+' {
s = s_next;
}
}
let mut exp_num = 0_i32;
s.parse_digits(|digit| {
if exp_num < 0x10000 {
exp_num = 10 * exp_num + digit as i32;
}
});
d.decimal_point += if neg_exp { -exp_num } else { exp_num };
}
}
for i in d.num_digits..DecimalSeq::MAX_DIGITS_WITHOUT_OVERFLOW {
d.digits[i] = 0;
}
d
}
fn number_of_digits_decimal_left_shift(d: &DecimalSeq, mut shift: usize) -> usize {
#[rustfmt::skip]
const TABLE: [u16; 65] = [
0x0000, 0x0800, 0x0801, 0x0803, 0x1006, 0x1009, 0x100D, 0x1812, 0x1817, 0x181D, 0x2024,
0x202B, 0x2033, 0x203C, 0x2846, 0x2850, 0x285B, 0x3067, 0x3073, 0x3080, 0x388E, 0x389C,
0x38AB, 0x38BB, 0x40CC, 0x40DD, 0x40EF, 0x4902, 0x4915, 0x4929, 0x513E, 0x5153, 0x5169,
0x5180, 0x5998, 0x59B0, 0x59C9, 0x61E3, 0x61FD, 0x6218, 0x6A34, 0x6A50, 0x6A6D, 0x6A8B,
0x72AA, 0x72C9, 0x72E9, 0x7B0A, 0x7B2B, 0x7B4D, 0x8370, 0x8393, 0x83B7, 0x83DC, 0x8C02,
0x8C28, 0x8C4F, 0x9477, 0x949F, 0x94C8, 0x9CF2, 0x051C, 0x051C, 0x051C, 0x051C,
];
#[rustfmt::skip]
const TABLE_POW5: [u8; 0x051C] = [
5, 2, 5, 1, 2, 5, 6, 2, 5, 3, 1, 2, 5, 1, 5, 6, 2, 5, 7, 8, 1, 2, 5, 3, 9, 0, 6, 2, 5, 1,
9, 5, 3, 1, 2, 5, 9, 7, 6, 5, 6, 2, 5, 4, 8, 8, 2, 8, 1, 2, 5, 2, 4, 4, 1, 4, 0, 6, 2, 5,
1, 2, 2, 0, 7, 0, 3, 1, 2, 5, 6, 1, 0, 3, 5, 1, 5, 6, 2, 5, 3, 0, 5, 1, 7, 5, 7, 8, 1, 2,
5, 1, 5, 2, 5, 8, 7, 8, 9, 0, 6, 2, 5, 7, 6, 2, 9, 3, 9, 4, 5, 3, 1, 2, 5, 3, 8, 1, 4, 6,
9, 7, 2, 6, 5, 6, 2, 5, 1, 9, 0, 7, 3, 4, 8, 6, 3, 2, 8, 1, 2, 5, 9, 5, 3, 6, 7, 4, 3, 1,
6, 4, 0, 6, 2, 5, 4, 7, 6, 8, 3, 7, 1, 5, 8, 2, 0, 3, 1, 2, 5, 2, 3, 8, 4, 1, 8, 5, 7, 9,
1, 0, 1, 5, 6, 2, 5, 1, 1, 9, 2, 0, 9, 2, 8, 9, 5, 5, 0, 7, 8, 1, 2, 5, 5, 9, 6, 0, 4, 6,
4, 4, 7, 7, 5, 3, 9, 0, 6, 2, 5, 2, 9, 8, 0, 2, 3, 2, 2, 3, 8, 7, 6, 9, 5, 3, 1, 2, 5, 1,
4, 9, 0, 1, 1, 6, 1, 1, 9, 3, 8, 4, 7, 6, 5, 6, 2, 5, 7, 4, 5, 0, 5, 8, 0, 5, 9, 6, 9, 2,
3, 8, 2, 8, 1, 2, 5, 3, 7, 2, 5, 2, 9, 0, 2, 9, 8, 4, 6, 1, 9, 1, 4, 0, 6, 2, 5, 1, 8, 6,
2, 6, 4, 5, 1, 4, 9, 2, 3, 0, 9, 5, 7, 0, 3, 1, 2, 5, 9, 3, 1, 3, 2, 2, 5, 7, 4, 6, 1, 5,
4, 7, 8, 5, 1, 5, 6, 2, 5, 4, 6, 5, 6, 6, 1, 2, 8, 7, 3, 0, 7, 7, 3, 9, 2, 5, 7, 8, 1, 2,
5, 2, 3, 2, 8, 3, 0, 6, 4, 3, 6, 5, 3, 8, 6, 9, 6, 2, 8, 9, 0, 6, 2, 5, 1, 1, 6, 4, 1, 5,
3, 2, 1, 8, 2, 6, 9, 3, 4, 8, 1, 4, 4, 5, 3, 1, 2, 5, 5, 8, 2, 0, 7, 6, 6, 0, 9, 1, 3, 4,
6, 7, 4, 0, 7, 2, 2, 6, 5, 6, 2, 5, 2, 9, 1, 0, 3, 8, 3, 0, 4, 5, 6, 7, 3, 3, 7, 0, 3, 6,
1, 3, 2, 8, 1, 2, 5, 1, 4, 5, 5, 1, 9, 1, 5, 2, 2, 8, 3, 6, 6, 8, 5, 1, 8, 0, 6, 6, 4, 0,
6, 2, 5, 7, 2, 7, 5, 9, 5, 7, 6, 1, 4, 1, 8, 3, 4, 2, 5, 9, 0, 3, 3, 2, 0, 3, 1, 2, 5, 3,
6, 3, 7, 9, 7, 8, 8, 0, 7, 0, 9, 1, 7, 1, 2, 9, 5, 1, 6, 6, 0, 1, 5, 6, 2, 5, 1, 8, 1, 8,
9, 8, 9, 4, 0, 3, 5, 4, 5, 8, 5, 6, 4, 7, 5, 8, 3, 0, 0, 7, 8, 1, 2, 5, 9, 0, 9, 4, 9, 4,
7, 0, 1, 7, 7, 2, 9, 2, 8, 2, 3, 7, 9, 1, 5, 0, 3, 9, 0, 6, 2, 5, 4, 5, 4, 7, 4, 7, 3, 5,
0, 8, 8, 6, 4, 6, 4, 1, 1, 8, 9, 5, 7, 5, 1, 9, 5, 3, 1, 2, 5, 2, 2, 7, 3, 7, 3, 6, 7, 5,
4, 4, 3, 2, 3, 2, 0, 5, 9, 4, 7, 8, 7, 5, 9, 7, 6, 5, 6, 2, 5, 1, 1, 3, 6, 8, 6, 8, 3, 7,
7, 2, 1, 6, 1, 6, 0, 2, 9, 7, 3, 9, 3, 7, 9, 8, 8, 2, 8, 1, 2, 5, 5, 6, 8, 4, 3, 4, 1, 8,
8, 6, 0, 8, 0, 8, 0, 1, 4, 8, 6, 9, 6, 8, 9, 9, 4, 1, 4, 0, 6, 2, 5, 2, 8, 4, 2, 1, 7, 0,
9, 4, 3, 0, 4, 0, 4, 0, 0, 7, 4, 3, 4, 8, 4, 4, 9, 7, 0, 7, 0, 3, 1, 2, 5, 1, 4, 2, 1, 0,
8, 5, 4, 7, 1, 5, 2, 0, 2, 0, 0, 3, 7, 1, 7, 4, 2, 2, 4, 8, 5, 3, 5, 1, 5, 6, 2, 5, 7, 1,
0, 5, 4, 2, 7, 3, 5, 7, 6, 0, 1, 0, 0, 1, 8, 5, 8, 7, 1, 1, 2, 4, 2, 6, 7, 5, 7, 8, 1, 2,
5, 3, 5, 5, 2, 7, 1, 3, 6, 7, 8, 8, 0, 0, 5, 0, 0, 9, 2, 9, 3, 5, 5, 6, 2, 1, 3, 3, 7, 8,
9, 0, 6, 2, 5, 1, 7, 7, 6, 3, 5, 6, 8, 3, 9, 4, 0, 0, 2, 5, 0, 4, 6, 4, 6, 7, 7, 8, 1, 0,
6, 6, 8, 9, 4, 5, 3, 1, 2, 5, 8, 8, 8, 1, 7, 8, 4, 1, 9, 7, 0, 0, 1, 2, 5, 2, 3, 2, 3, 3,
8, 9, 0, 5, 3, 3, 4, 4, 7, 2, 6, 5, 6, 2, 5, 4, 4, 4, 0, 8, 9, 2, 0, 9, 8, 5, 0, 0, 6, 2,
6, 1, 6, 1, 6, 9, 4, 5, 2, 6, 6, 7, 2, 3, 6, 3, 2, 8, 1, 2, 5, 2, 2, 2, 0, 4, 4, 6, 0, 4,
9, 2, 5, 0, 3, 1, 3, 0, 8, 0, 8, 4, 7, 2, 6, 3, 3, 3, 6, 1, 8, 1, 6, 4, 0, 6, 2, 5, 1, 1,
1, 0, 2, 2, 3, 0, 2, 4, 6, 2, 5, 1, 5, 6, 5, 4, 0, 4, 2, 3, 6, 3, 1, 6, 6, 8, 0, 9, 0, 8,
2, 0, 3, 1, 2, 5, 5, 5, 5, 1, 1, 1, 5, 1, 2, 3, 1, 2, 5, 7, 8, 2, 7, 0, 2, 1, 1, 8, 1, 5,
8, 3, 4, 0, 4, 5, 4, 1, 0, 1, 5, 6, 2, 5, 2, 7, 7, 5, 5, 5, 7, 5, 6, 1, 5, 6, 2, 8, 9, 1,
3, 5, 1, 0, 5, 9, 0, 7, 9, 1, 7, 0, 2, 2, 7, 0, 5, 0, 7, 8, 1, 2, 5, 1, 3, 8, 7, 7, 7, 8,
7, 8, 0, 7, 8, 1, 4, 4, 5, 6, 7, 5, 5, 2, 9, 5, 3, 9, 5, 8, 5, 1, 1, 3, 5, 2, 5, 3, 9, 0,
6, 2, 5, 6, 9, 3, 8, 8, 9, 3, 9, 0, 3, 9, 0, 7, 2, 2, 8, 3, 7, 7, 6, 4, 7, 6, 9, 7, 9, 2,
5, 5, 6, 7, 6, 2, 6, 9, 5, 3, 1, 2, 5, 3, 4, 6, 9, 4, 4, 6, 9, 5, 1, 9, 5, 3, 6, 1, 4, 1,
8, 8, 8, 2, 3, 8, 4, 8, 9, 6, 2, 7, 8, 3, 8, 1, 3, 4, 7, 6, 5, 6, 2, 5, 1, 7, 3, 4, 7, 2,
3, 4, 7, 5, 9, 7, 6, 8, 0, 7, 0, 9, 4, 4, 1, 1, 9, 2, 4, 4, 8, 1, 3, 9, 1, 9, 0, 6, 7, 3,
8, 2, 8, 1, 2, 5, 8, 6, 7, 3, 6, 1, 7, 3, 7, 9, 8, 8, 4, 0, 3, 5, 4, 7, 2, 0, 5, 9, 6, 2,
2, 4, 0, 6, 9, 5, 9, 5, 3, 3, 6, 9, 1, 4, 0, 6, 2, 5,
];
shift &= 63;
let x_a = TABLE[shift];
let x_b = TABLE[shift + 1];
let num_new_digits = (x_a >> 11) as _;
let pow5_a = (0x7FF & x_a) as usize;
let pow5_b = (0x7FF & x_b) as usize;
let pow5 = &TABLE_POW5[pow5_a..];
for (i, &p5) in pow5.iter().enumerate().take(pow5_b - pow5_a) {
if i >= d.num_digits {
return num_new_digits - 1;
} else if d.digits[i] == p5 {
continue;
} else if d.digits[i] < p5 {
return num_new_digits - 1;
} else {
return num_new_digits;
}
}
num_new_digits
}

View File

@ -1,14 +1,57 @@
//! Helper trait for generic float types.
use core::f64;
use crate::fmt::{Debug, LowerExp};
use crate::num::FpCategory;
use crate::ops::{Add, Div, Mul, Neg};
use crate::ops::{self, Add, Div, Mul, Neg};
/// A helper trait to avoid duplicating basically all the conversion code for `f32` and `f64`.
/// Lossy `as` casting between two types.
pub trait CastInto<T: Copy>: Copy {
fn cast(self) -> T;
}
/// Collection of traits that allow us to be generic over integer size.
pub trait Integer:
Sized
+ Clone
+ Copy
+ Debug
+ ops::Shr<u32, Output = Self>
+ ops::Shl<u32, Output = Self>
+ ops::BitAnd<Output = Self>
+ ops::BitOr<Output = Self>
+ PartialEq
+ CastInto<i16>
{
const ZERO: Self;
const ONE: Self;
}
macro_rules! int {
($($ty:ty),+) => {
$(
impl CastInto<i16> for $ty {
fn cast(self) -> i16 {
self as i16
}
}
impl Integer for $ty {
const ZERO: Self = 0;
const ONE: Self = 1;
}
)+
}
}
int!(u32, u64);
/// A helper trait to avoid duplicating basically all the conversion code for IEEE floats.
///
/// See the parent module's doc comment for why this is necessary.
///
/// Should **never ever** be implemented for other types or be used outside the dec2flt module.
/// Should **never ever** be implemented for other types or be used outside the `dec2flt` module.
#[doc(hidden)]
pub trait RawFloat:
Sized
@ -24,62 +67,107 @@ pub trait RawFloat:
+ Copy
+ Debug
{
/// The unsigned integer with the same size as the float
type Int: Integer + Into<u64>;
/* general constants */
const INFINITY: Self;
const NEG_INFINITY: Self;
const NAN: Self;
const NEG_NAN: Self;
/// The number of bits in the significand, *excluding* the hidden bit.
const MANTISSA_EXPLICIT_BITS: usize;
/// Bit width of the float
const BITS: u32;
// Round-to-even only happens for negative values of q
// when q ≥ 4 in the 64-bit case and when q ≥ 17 in
// the 32-bitcase.
//
// When q ≥ 0,we have that 5^q ≤ 2m+1. In the 64-bit case,we
// have 5^q ≤ 2m+1 ≤ 2^54 or q ≤ 23. In the 32-bit case,we have
// 5^q ≤ 2m+1 ≤ 2^25 or q ≤ 10.
//
// When q < 0, we have w ≥ (2m+1)×5^q. We must have that w < 2^64
// so (2m+1)×5^q < 2^64. We have that 2m+1 > 2^53 (64-bit case)
// or 2m+1 > 2^24 (32-bit case). Hence,we must have 2^53×5^q < 2^64
// (64-bit) and 2^24×5^q < 2^64 (32-bit). Hence we have 5^q < 2^11
// or q ≥ 4 (64-bit case) and 5^q < 2^40 or q ≥ 17 (32-bitcase).
//
// Thus we have that we only need to round ties to even when
// we have that q ∈ [4,23](in the 64-bit case) or q∈[17,10]
// (in the 32-bit case). In both cases,the power of five(5^|q|)
// fits in a 64-bit word.
/// The number of bits in the significand, *including* the hidden bit.
const SIG_TOTAL_BITS: u32;
const EXP_MASK: Self::Int;
const SIG_MASK: Self::Int;
/// The number of bits in the significand, *excluding* the hidden bit.
const SIG_BITS: u32 = Self::SIG_TOTAL_BITS - 1;
/// Number of bits in the exponent.
const EXP_BITS: u32 = Self::BITS - Self::SIG_BITS - 1;
/// The saturated (maximum bitpattern) value of the exponent, i.e. the infinite
/// representation.
///
/// This shifted fully right, use `EXP_MASK` for the shifted value.
const EXP_SAT: u32 = (1 << Self::EXP_BITS) - 1;
/// Signed version of `EXP_SAT` since we convert a lot.
const INFINITE_POWER: i32 = Self::EXP_SAT as i32;
/// The exponent bias value. This is also the maximum value of the exponent.
const EXP_BIAS: u32 = Self::EXP_SAT >> 1;
/// Minimum exponent value of normal values.
const EXP_MIN: i32 = -(Self::EXP_BIAS as i32 - 1);
/// Round-to-even only happens for negative values of q
/// when q ≥ 4 in the 64-bit case and when q ≥ 17 in
/// the 32-bitcase.
///
/// When q ≥ 0,we have that 5^q ≤ 2m+1. In the 64-bit case,we
/// have 5^q ≤ 2m+1 ≤ 2^54 or q ≤ 23. In the 32-bit case,we have
/// 5^q ≤ 2m+1 ≤ 2^25 or q ≤ 10.
///
/// When q < 0, we have w ≥ (2m+1)×5^q. We must have that w < 2^64
/// so (2m+1)×5^q < 2^64. We have that 2m+1 > 2^53 (64-bit case)
/// or 2m+1 > 2^24 (32-bit case). Hence,we must have 2^53×5^q < 2^64
/// (64-bit) and 2^24×5^q < 2^64 (32-bit). Hence we have 5^q < 2^11
/// or q ≥ 4 (64-bit case) and 5^q < 2^40 or q ≥ 17 (32-bitcase).
///
/// Thus we have that we only need to round ties to even when
/// we have that q ∈ [4,23](in the 64-bit case) or q∈[17,10]
/// (in the 32-bit case). In both cases,the power of five(5^|q|)
/// fits in a 64-bit word.
const MIN_EXPONENT_ROUND_TO_EVEN: i32;
const MAX_EXPONENT_ROUND_TO_EVEN: i32;
// Minimum exponent that for a fast path case, or `-⌊(MANTISSA_EXPLICIT_BITS+1)/log2(5)⌋`
const MIN_EXPONENT_FAST_PATH: i64;
/* limits related to Fast pathing */
// Maximum exponent that for a fast path case, or `⌊(MANTISSA_EXPLICIT_BITS+1)/log2(5)⌋`
const MAX_EXPONENT_FAST_PATH: i64;
/// Largest decimal exponent for a non-infinite value.
///
/// This is the max exponent in binary converted to the max exponent in decimal. Allows fast
/// pathing anything larger than `10^LARGEST_POWER_OF_TEN`, which will round to infinity.
const LARGEST_POWER_OF_TEN: i32 = {
let largest_pow2 = Self::EXP_BIAS + 1;
pow2_to_pow10(largest_pow2 as i64) as i32
};
// Maximum exponent that can be represented for a disguised-fast path case.
// This is `MAX_EXPONENT_FAST_PATH + ⌊(MANTISSA_EXPLICIT_BITS+1)/log2(10)⌋`
const MAX_EXPONENT_DISGUISED_FAST_PATH: i64;
// Minimum exponent value `-(1 << (EXP_BITS - 1)) + 1`.
const MINIMUM_EXPONENT: i32;
// Largest exponent value `(1 << EXP_BITS) - 1`.
const INFINITE_POWER: i32;
// Index (in bits) of the sign.
const SIGN_INDEX: usize;
// Smallest decimal exponent for a non-zero value.
/// Smallest decimal exponent for a non-zero value. This allows for fast pathing anything
/// smaller than `10^SMALLEST_POWER_OF_TEN`, which will round to zero.
///
/// The smallest power of ten is represented by `⌊log10(2^-n / (2^64 - 1))⌋`, where `n` is
/// the smallest power of two. The `2^64 - 1)` denomenator comes from the number of values
/// that are representable by the intermediate storage format. I don't actually know _why_
/// the storage format is relevant here.
///
/// The values may be calculated using the formula. Unfortunately we cannot calculate them at
/// compile time since intermediates exceed the range of an `f64`.
const SMALLEST_POWER_OF_TEN: i32;
// Largest decimal exponent for a non-infinite value.
const LARGEST_POWER_OF_TEN: i32;
/// Maximum exponent for a fast path case, or `⌊(SIG_BITS+1)/log2(5)⌋`
// assuming FLT_EVAL_METHOD = 0
const MAX_EXPONENT_FAST_PATH: i64 = {
let log2_5 = f64::consts::LOG2_10 - 1.0;
(Self::SIG_TOTAL_BITS as f64 / log2_5) as i64
};
// Maximum mantissa for the fast-path (`1 << 53` for f64).
const MAX_MANTISSA_FAST_PATH: u64 = 2_u64 << Self::MANTISSA_EXPLICIT_BITS;
/// Minimum exponent for a fast path case, or `-⌊(SIG_BITS+1)/log2(5)⌋`
const MIN_EXPONENT_FAST_PATH: i64 = -Self::MAX_EXPONENT_FAST_PATH;
/// Maximum exponent that can be represented for a disguised-fast path case.
/// This is `MAX_EXPONENT_FAST_PATH + ⌊(SIG_BITS+1)/log2(10)⌋`
const MAX_EXPONENT_DISGUISED_FAST_PATH: i64 =
Self::MAX_EXPONENT_FAST_PATH + (Self::SIG_TOTAL_BITS as f64 / f64::consts::LOG2_10) as i64;
/// Maximum mantissa for the fast-path (`1 << 53` for f64).
const MAX_MANTISSA_FAST_PATH: u64 = 1 << Self::SIG_TOTAL_BITS;
/// Converts integer into float through an as cast.
/// This is only called in the fast-path algorithm, and therefore
@ -96,27 +184,51 @@ pub trait RawFloat:
/// Returns the category that this number falls into.
fn classify(self) -> FpCategory;
/// Transmute to the integer representation
fn to_bits(self) -> Self::Int;
/// Returns the mantissa, exponent and sign as integers.
fn integer_decode(self) -> (u64, i16, i8);
///
/// That is, this returns `(m, p, s)` such that `s * m * 2^p` represents the original float.
/// For 0, the exponent will be `-(EXP_BIAS + SIG_BITS`, which is the
/// minimum subnormal power.
fn integer_decode(self) -> (u64, i16, i8) {
let bits = self.to_bits();
let sign: i8 = if bits >> (Self::BITS - 1) == Self::Int::ZERO { 1 } else { -1 };
let mut exponent: i16 = ((bits & Self::EXP_MASK) >> Self::SIG_BITS).cast();
let mantissa = if exponent == 0 {
(bits & Self::SIG_MASK) << 1
} else {
(bits & Self::SIG_MASK) | (Self::Int::ONE << Self::SIG_BITS)
};
// Exponent bias + mantissa shift
exponent -= (Self::EXP_BIAS + Self::SIG_BITS) as i16;
(mantissa.into(), exponent, sign)
}
}
/// Solve for `b` in `10^b = 2^a`
const fn pow2_to_pow10(a: i64) -> i64 {
let res = (a as f64) / f64::consts::LOG2_10;
res as i64
}
impl RawFloat for f32 {
type Int = u32;
const INFINITY: Self = f32::INFINITY;
const NEG_INFINITY: Self = f32::NEG_INFINITY;
const NAN: Self = f32::NAN;
const NEG_NAN: Self = -f32::NAN;
const MANTISSA_EXPLICIT_BITS: usize = 23;
const BITS: u32 = 32;
const SIG_TOTAL_BITS: u32 = Self::MANTISSA_DIGITS;
const EXP_MASK: Self::Int = Self::EXP_MASK;
const SIG_MASK: Self::Int = Self::MAN_MASK;
const MIN_EXPONENT_ROUND_TO_EVEN: i32 = -17;
const MAX_EXPONENT_ROUND_TO_EVEN: i32 = 10;
const MIN_EXPONENT_FAST_PATH: i64 = -10; // assuming FLT_EVAL_METHOD = 0
const MAX_EXPONENT_FAST_PATH: i64 = 10;
const MAX_EXPONENT_DISGUISED_FAST_PATH: i64 = 17;
const MINIMUM_EXPONENT: i32 = -127;
const INFINITE_POWER: i32 = 0xFF;
const SIGN_INDEX: usize = 31;
const SMALLEST_POWER_OF_TEN: i32 = -65;
const LARGEST_POWER_OF_TEN: i32 = 38;
#[inline]
fn from_u64(v: u64) -> Self {
@ -136,16 +248,8 @@ impl RawFloat for f32 {
TABLE[exponent & 15]
}
/// Returns the mantissa, exponent and sign as integers.
fn integer_decode(self) -> (u64, i16, i8) {
let bits = self.to_bits();
let sign: i8 = if bits >> 31 == 0 { 1 } else { -1 };
let mut exponent: i16 = ((bits >> 23) & 0xff) as i16;
let mantissa =
if exponent == 0 { (bits & 0x7fffff) << 1 } else { (bits & 0x7fffff) | 0x800000 };
// Exponent bias + mantissa shift
exponent -= 127 + 23;
(mantissa as u64, exponent, sign)
fn to_bits(self) -> Self::Int {
self.to_bits()
}
fn classify(self) -> FpCategory {
@ -154,22 +258,21 @@ impl RawFloat for f32 {
}
impl RawFloat for f64 {
const INFINITY: Self = f64::INFINITY;
const NEG_INFINITY: Self = f64::NEG_INFINITY;
const NAN: Self = f64::NAN;
const NEG_NAN: Self = -f64::NAN;
type Int = u64;
const INFINITY: Self = Self::INFINITY;
const NEG_INFINITY: Self = Self::NEG_INFINITY;
const NAN: Self = Self::NAN;
const NEG_NAN: Self = -Self::NAN;
const BITS: u32 = 64;
const SIG_TOTAL_BITS: u32 = Self::MANTISSA_DIGITS;
const EXP_MASK: Self::Int = Self::EXP_MASK;
const SIG_MASK: Self::Int = Self::MAN_MASK;
const MANTISSA_EXPLICIT_BITS: usize = 52;
const MIN_EXPONENT_ROUND_TO_EVEN: i32 = -4;
const MAX_EXPONENT_ROUND_TO_EVEN: i32 = 23;
const MIN_EXPONENT_FAST_PATH: i64 = -22; // assuming FLT_EVAL_METHOD = 0
const MAX_EXPONENT_FAST_PATH: i64 = 22;
const MAX_EXPONENT_DISGUISED_FAST_PATH: i64 = 37;
const MINIMUM_EXPONENT: i32 = -1023;
const INFINITE_POWER: i32 = 0x7FF;
const SIGN_INDEX: usize = 63;
const SMALLEST_POWER_OF_TEN: i32 = -342;
const LARGEST_POWER_OF_TEN: i32 = 308;
#[inline]
fn from_u64(v: u64) -> Self {
@ -190,19 +293,8 @@ impl RawFloat for f64 {
TABLE[exponent & 31]
}
/// Returns the mantissa, exponent and sign as integers.
fn integer_decode(self) -> (u64, i16, i8) {
let bits = self.to_bits();
let sign: i8 = if bits >> 63 == 0 { 1 } else { -1 };
let mut exponent: i16 = ((bits >> 52) & 0x7ff) as i16;
let mantissa = if exponent == 0 {
(bits & 0xfffffffffffff) << 1
} else {
(bits & 0xfffffffffffff) | 0x10000000000000
};
// Exponent bias + mantissa shift
exponent -= 1023 + 52;
(mantissa, exponent, sign)
fn to_bits(self) -> Self::Int {
self.to_bits()
}
fn classify(self) -> FpCategory {

View File

@ -38,7 +38,7 @@ pub fn compute_float<F: RawFloat>(q: i64, mut w: u64) -> BiasedFp {
// Normalize our significant digits, so the most-significant bit is set.
let lz = w.leading_zeros();
w <<= lz;
let (lo, hi) = compute_product_approx(q, w, F::MANTISSA_EXPLICIT_BITS + 3);
let (lo, hi) = compute_product_approx(q, w, F::SIG_BITS 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
@ -61,8 +61,8 @@ pub fn compute_float<F: RawFloat>(q: i64, mut w: u64) -> BiasedFp {
}
}
let upperbit = (hi >> 63) as i32;
let mut mantissa = hi >> (upperbit + 64 - F::MANTISSA_EXPLICIT_BITS as i32 - 3);
let mut power2 = power(q as i32) + upperbit - lz as i32 - F::MINIMUM_EXPONENT;
let mut mantissa = hi >> (upperbit + 64 - F::SIG_BITS as i32 - 3);
let mut power2 = power(q as i32) + upperbit - lz as i32 - F::EXP_MIN + 1;
if power2 <= 0 {
if -power2 + 1 >= 64 {
// Have more than 64 bits below the minimum exponent, must be 0.
@ -72,8 +72,8 @@ pub fn compute_float<F: RawFloat>(q: i64, mut w: u64) -> BiasedFp {
mantissa >>= -power2 + 1;
mantissa += mantissa & 1;
mantissa >>= 1;
power2 = (mantissa >= (1_u64 << F::MANTISSA_EXPLICIT_BITS)) as i32;
return BiasedFp { f: mantissa, e: power2 };
power2 = (mantissa >= (1_u64 << F::SIG_BITS)) as i32;
return BiasedFp { m: mantissa, p_biased: power2 };
}
// Need to handle rounding ties. Normally, we need to round up,
// but if we fall right in between and we have an even basis, we
@ -89,8 +89,8 @@ pub fn compute_float<F: RawFloat>(q: i64, mut w: u64) -> BiasedFp {
if lo <= 1
&& q >= F::MIN_EXPONENT_ROUND_TO_EVEN as i64
&& q <= F::MAX_EXPONENT_ROUND_TO_EVEN as i64
&& mantissa & 3 == 1
&& (mantissa << (upperbit + 64 - F::MANTISSA_EXPLICIT_BITS as i32 - 3)) == hi
&& mantissa & 0b11 == 0b01
&& (mantissa << (upperbit + 64 - F::SIG_BITS as i32 - 3)) == hi
{
// Zero the lowest bit, so we don't round up.
mantissa &= !1_u64;
@ -98,20 +98,20 @@ pub fn compute_float<F: RawFloat>(q: i64, mut w: u64) -> BiasedFp {
// Round-to-even, then shift the significant digits into place.
mantissa += mantissa & 1;
mantissa >>= 1;
if mantissa >= (2_u64 << F::MANTISSA_EXPLICIT_BITS) {
if mantissa >= (2_u64 << F::SIG_BITS) {
// 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_EXPLICIT_BITS;
mantissa = 1_u64 << F::SIG_BITS;
power2 += 1;
}
// Zero out the hidden bit.
mantissa &= !(1_u64 << F::MANTISSA_EXPLICIT_BITS);
mantissa &= !(1_u64 << F::SIG_BITS);
if power2 >= F::INFINITE_POWER {
// Exponent is above largest normal value, must be infinite.
return fp_inf;
}
BiasedFp { f: mantissa, e: power2 }
BiasedFp { m: mantissa, p_biased: power2 }
}
/// Calculate a base 2 exponent from a decimal exponent.

View File

@ -3,8 +3,8 @@
//! # Problem statement
//!
//! We are given a decimal string such as `12.34e56`. This string consists of integral (`12`),
//! fractional (`34`), and exponent (`56`) parts. All parts are optional and interpreted as zero
//! when missing.
//! fractional (`34`), and exponent (`56`) parts. All parts are optional and interpreted as a
//! default value (1 or 0) when missing.
//!
//! We seek the IEEE 754 floating point number that is closest to the exact value of the decimal
//! string. It is well-known that many decimal strings do not have terminating representations in
@ -67,6 +67,18 @@
//! "such that the exponent +/- the number of decimal digits fits into a 64 bit integer".
//! Larger exponents are accepted, but we don't do arithmetic with them, they are immediately
//! turned into {positive,negative} {zero,infinity}.
//!
//! # Notation
//!
//! This module uses the same notation as the Lemire paper:
//!
//! - `m`: binary mantissa; always nonnegative
//! - `p`: binary exponent; a signed integer
//! - `w`: decimal significand; always nonnegative
//! - `q`: decimal exponent; a signed integer
//!
//! This gives `m * 2^p` for the binary floating-point number, with `w * 10^q` as the decimal
//! equivalent.
#![doc(hidden)]
#![unstable(
@ -85,14 +97,14 @@ use crate::fmt;
use crate::str::FromStr;
mod common;
mod decimal;
pub mod decimal;
pub mod decimal_seq;
mod fpu;
mod slow;
mod table;
// float is used in flt2dec, and all are used in unit tests.
pub mod float;
pub mod lemire;
pub mod number;
pub mod parse;
macro_rules! from_str_float_impl {
@ -220,10 +232,10 @@ pub fn pfe_invalid() -> ParseFloatError {
}
/// Converts a `BiasedFp` to the closest machine float type.
fn biased_fp_to_float<T: RawFloat>(x: BiasedFp) -> T {
let mut word = x.f;
word |= (x.e as u64) << T::MANTISSA_EXPLICIT_BITS;
T::from_u64_bits(word)
fn biased_fp_to_float<F: RawFloat>(x: BiasedFp) -> F {
let mut word = x.m;
word |= (x.p_biased as u64) << F::SIG_BITS;
F::from_u64_bits(word)
}
/// Converts a decimal string into a floating point number.
@ -260,12 +272,15 @@ pub fn dec2flt<F: RawFloat>(s: &str) -> Result<F, ParseFloatError> {
// 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.e >= 0 && fp != compute_float::<F>(num.exponent, num.mantissa + 1) {
fp.e = -1;
if num.many_digits
&& fp.p_biased >= 0
&& fp != compute_float::<F>(num.exponent, num.mantissa + 1)
{
fp.p_biased = -1;
}
// Unable to correctly round the float using the Eisel-Lemire algorithm.
// Fallback to a slower, but always correct algorithm.
if fp.e < 0 {
if fp.p_biased < 0 {
fp = parse_long_mantissa::<F>(s);
}

View File

@ -1,88 +0,0 @@
//! Representation of a float as the significant digits and exponent.
use crate::num::dec2flt::float::RawFloat;
use crate::num::dec2flt::fpu::set_precision;
#[rustfmt::skip]
const INT_POW10: [u64; 16] = [
1,
10,
100,
1000,
10000,
100000,
1000000,
10000000,
100000000,
1000000000,
10000000000,
100000000000,
1000000000000,
10000000000000,
100000000000000,
1000000000000000,
];
#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
pub struct Number {
pub exponent: i64,
pub mantissa: u64,
pub negative: bool,
pub many_digits: bool,
}
impl Number {
/// Detect if the float can be accurately reconstructed from native floats.
#[inline]
fn is_fast_path<F: RawFloat>(&self) -> bool {
F::MIN_EXPONENT_FAST_PATH <= self.exponent
&& self.exponent <= F::MAX_EXPONENT_DISGUISED_FAST_PATH
&& self.mantissa <= F::MAX_MANTISSA_FAST_PATH
&& !self.many_digits
}
/// The fast path algorithm using machine-sized integers and floats.
///
/// This is extracted into a separate function so that it can be attempted before constructing
/// a Decimal. This only works if both the mantissa and the exponent
/// can be exactly represented as a machine float, since IEE-754 guarantees
/// no rounding will occur.
///
/// There is an exception: disguised fast-path cases, where we can shift
/// powers-of-10 from the exponent to the significant digits.
pub fn try_fast_path<F: RawFloat>(&self) -> Option<F> {
// Here we need to work around <https://github.com/rust-lang/rust/issues/114479>.
// The fast path crucially depends on arithmetic being rounded to the correct number of bits
// without any intermediate rounding. On x86 (without SSE or SSE2) this requires the precision
// of the x87 FPU stack to be changed so that it directly rounds to 64/32 bit.
// The `set_precision` function takes care of setting the precision on architectures which
// require setting it by changing the global state (like the control word of the x87 FPU).
let _cw = set_precision::<F>();
if self.is_fast_path::<F>() {
let mut value = if self.exponent <= F::MAX_EXPONENT_FAST_PATH {
// normal fast path
let value = F::from_u64(self.mantissa);
if self.exponent < 0 {
value / F::pow10_fast_path((-self.exponent) as _)
} else {
value * F::pow10_fast_path(self.exponent as _)
}
} else {
// disguised fast path
let shift = self.exponent - F::MAX_EXPONENT_FAST_PATH;
let mantissa = self.mantissa.checked_mul(INT_POW10[shift as usize])?;
if mantissa > F::MAX_MANTISSA_FAST_PATH {
return None;
}
F::from_u64(mantissa) * F::pow10_fast_path(F::MAX_EXPONENT_FAST_PATH as _)
};
if self.negative {
value = -value;
}
Some(value)
} else {
None
}
}
}

View File

@ -1,8 +1,8 @@
//! Functions to parse floating-point numbers.
use crate::num::dec2flt::common::{ByteSlice, is_8digits};
use crate::num::dec2flt::decimal::Decimal;
use crate::num::dec2flt::float::RawFloat;
use crate::num::dec2flt::number::Number;
const MIN_19DIGIT_INT: u64 = 100_0000_0000_0000_0000;
@ -100,7 +100,7 @@ fn parse_scientific(s_ref: &mut &[u8]) -> Option<i64> {
///
/// This creates a representation of the float as the
/// significant digits and the decimal exponent.
fn parse_partial_number(mut s: &[u8]) -> Option<(Number, usize)> {
fn parse_partial_number(mut s: &[u8]) -> Option<(Decimal, usize)> {
debug_assert!(!s.is_empty());
// parse initial digits before dot
@ -146,7 +146,7 @@ fn parse_partial_number(mut s: &[u8]) -> Option<(Number, usize)> {
// handle uncommon case with many digits
if n_digits <= 19 {
return Some((Number { exponent, mantissa, negative: false, many_digits: false }, len));
return Some((Decimal { exponent, mantissa, negative: false, many_digits: false }, len));
}
n_digits -= 19;
@ -179,13 +179,13 @@ fn parse_partial_number(mut s: &[u8]) -> Option<(Number, usize)> {
exponent += exp_number;
}
Some((Number { exponent, mantissa, negative: false, many_digits }, len))
Some((Decimal { exponent, mantissa, negative: false, many_digits }, len))
}
/// Try to parse a non-special floating point number,
/// as well as two slices with integer and fractional parts
/// and the parsed exponent.
pub fn parse_number(s: &[u8]) -> Option<Number> {
pub fn parse_number(s: &[u8]) -> Option<Decimal> {
if let Some((float, rest)) = parse_partial_number(s) {
if rest == s.len() {
return Some(float);

View File

@ -1,7 +1,7 @@
//! Slow, fallback algorithm for cases the Eisel-Lemire algorithm cannot round.
use crate::num::dec2flt::common::BiasedFp;
use crate::num::dec2flt::decimal::{Decimal, parse_decimal};
use crate::num::dec2flt::decimal_seq::{DecimalSeq, parse_decimal_seq};
use crate::num::dec2flt::float::RawFloat;
/// Parse the significant digits and biased, binary exponent of a float.
@ -36,7 +36,7 @@ pub(crate) fn parse_long_mantissa<F: RawFloat>(s: &[u8]) -> BiasedFp {
let fp_zero = BiasedFp::zero_pow2(0);
let fp_inf = BiasedFp::zero_pow2(F::INFINITE_POWER);
let mut d = parse_decimal(s);
let mut d = parse_decimal_seq(s);
// Short-circuit if the value can only be a literal 0 or infinity.
if d.num_digits == 0 || d.decimal_point < -324 {
@ -50,7 +50,7 @@ pub(crate) fn parse_long_mantissa<F: RawFloat>(s: &[u8]) -> BiasedFp {
let n = d.decimal_point as usize;
let shift = get_shift(n);
d.right_shift(shift);
if d.decimal_point < -Decimal::DECIMAL_POINT_RANGE {
if d.decimal_point < -DecimalSeq::DECIMAL_POINT_RANGE {
return fp_zero;
}
exp2 += shift as i32;
@ -67,43 +67,43 @@ pub(crate) fn parse_long_mantissa<F: RawFloat>(s: &[u8]) -> BiasedFp {
get_shift((-d.decimal_point) as _)
};
d.left_shift(shift);
if d.decimal_point > Decimal::DECIMAL_POINT_RANGE {
if d.decimal_point > DecimalSeq::DECIMAL_POINT_RANGE {
return fp_inf;
}
exp2 -= shift as i32;
}
// We are now in the range [1/2 ... 1] but the binary format uses [1 ... 2].
exp2 -= 1;
while (F::MINIMUM_EXPONENT + 1) > exp2 {
let mut n = ((F::MINIMUM_EXPONENT + 1) - exp2) as usize;
while F::EXP_MIN > exp2 {
let mut n = (F::EXP_MIN - exp2) as usize;
if n > MAX_SHIFT {
n = MAX_SHIFT;
}
d.right_shift(n);
exp2 += n as i32;
}
if (exp2 - F::MINIMUM_EXPONENT) >= F::INFINITE_POWER {
if (exp2 - F::EXP_MIN + 1) >= F::INFINITE_POWER {
return fp_inf;
}
// Shift the decimal to the hidden bit, and then round the value
// to get the high mantissa+1 bits.
d.left_shift(F::MANTISSA_EXPLICIT_BITS + 1);
d.left_shift(F::SIG_BITS as usize + 1);
let mut mantissa = d.round();
if mantissa >= (1_u64 << (F::MANTISSA_EXPLICIT_BITS + 1)) {
if mantissa >= (1_u64 << (F::SIG_BITS + 1)) {
// Rounding up overflowed to the carry bit, need to
// shift back to the hidden bit.
d.right_shift(1);
exp2 += 1;
mantissa = d.round();
if (exp2 - F::MINIMUM_EXPONENT) >= F::INFINITE_POWER {
if (exp2 - F::EXP_MIN + 1) >= F::INFINITE_POWER {
return fp_inf;
}
}
let mut power2 = exp2 - F::MINIMUM_EXPONENT;
if mantissa < (1_u64 << F::MANTISSA_EXPLICIT_BITS) {
let mut power2 = exp2 - F::EXP_MIN + 1;
if mantissa < (1_u64 << F::SIG_BITS) {
power2 -= 1;
}
// Zero out all the bits above the explicit mantissa bits.
mantissa &= (1_u64 << F::MANTISSA_EXPLICIT_BITS) - 1;
BiasedFp { f: mantissa, e: power2 }
mantissa &= (1_u64 << F::SIG_BITS) - 1;
BiasedFp { m: mantissa, p_biased: power2 }
}

View File

@ -493,13 +493,13 @@ impl f32 {
pub const NEG_INFINITY: f32 = -1.0_f32 / 0.0_f32;
/// Sign bit
const SIGN_MASK: u32 = 0x8000_0000;
pub(crate) const SIGN_MASK: u32 = 0x8000_0000;
/// Exponent mask
const EXP_MASK: u32 = 0x7f80_0000;
pub(crate) const EXP_MASK: u32 = 0x7f80_0000;
/// Mantissa mask
const MAN_MASK: u32 = 0x007f_ffff;
pub(crate) const MAN_MASK: u32 = 0x007f_ffff;
/// Minimum representable positive value (min subnormal)
const TINY_BITS: u32 = 0x1;

View File

@ -492,13 +492,13 @@ impl f64 {
pub const NEG_INFINITY: f64 = -1.0_f64 / 0.0_f64;
/// Sign bit
const SIGN_MASK: u64 = 0x8000_0000_0000_0000;
pub(crate) const SIGN_MASK: u64 = 0x8000_0000_0000_0000;
/// Exponent mask
const EXP_MASK: u64 = 0x7ff0_0000_0000_0000;
pub(crate) const EXP_MASK: u64 = 0x7ff0_0000_0000_0000;
/// Mantissa mask
const MAN_MASK: u64 = 0x000f_ffff_ffff_ffff;
pub(crate) const MAN_MASK: u64 = 0x000f_ffff_ffff_ffff;
/// Minimum representable positive value (min subnormal)
const TINY_BITS: u64 = 0x1;

View File

@ -0,0 +1,28 @@
use core::num::dec2flt::decimal::Decimal;
type FPath<F> = ((i64, u64, bool, bool), Option<F>);
const FPATHS_F32: &[FPath<f32>] =
&[((0, 0, false, false), Some(0.0)), ((0, 0, false, false), Some(0.0))];
const FPATHS_F64: &[FPath<f64>] =
&[((0, 0, false, false), Some(0.0)), ((0, 0, false, false), Some(0.0))];
#[test]
fn check_fast_path_f32() {
for ((exponent, mantissa, negative, many_digits), expected) in FPATHS_F32.iter().copied() {
let dec = Decimal { exponent, mantissa, negative, many_digits };
let actual = dec.try_fast_path::<f32>();
assert_eq!(actual, expected);
}
}
#[test]
fn check_fast_path_f64() {
for ((exponent, mantissa, negative, many_digits), expected) in FPATHS_F64.iter().copied() {
let dec = Decimal { exponent, mantissa, negative, many_digits };
let actual = dec.try_fast_path::<f64>();
assert_eq!(actual, expected);
}
}

View File

@ -0,0 +1,30 @@
use core::num::dec2flt::decimal_seq::{DecimalSeq, parse_decimal_seq};
#[test]
fn test_trim() {
let mut dec = DecimalSeq::default();
let digits = [1, 2, 3, 4];
dec.digits[0..4].copy_from_slice(&digits);
dec.num_digits = 8;
dec.trim();
assert_eq!(dec.digits[0..4], digits);
assert_eq!(dec.num_digits, 4);
}
#[test]
fn test_parse() {
let tests = [("1.234", [1, 2, 3, 4], 1)];
for (s, exp_digits, decimal_point) in tests {
let actual = parse_decimal_seq(s.as_bytes());
let mut digits = [0; DecimalSeq::MAX_DIGITS];
digits[..exp_digits.len()].copy_from_slice(&exp_digits);
let expected =
DecimalSeq { num_digits: exp_digits.len(), decimal_point, truncated: false, digits };
assert_eq!(actual, expected);
}
}

View File

@ -12,8 +12,8 @@ fn test_f32_integer_decode() {
// Ignore the "sign" (quiet / signalling flag) of NAN.
// It can vary between runtime operations and LLVM folding.
let (nan_m, nan_e, _nan_s) = f32::NAN.integer_decode();
assert_eq!((nan_m, nan_e), (12582912, 105));
let (nan_m, nan_p, _nan_s) = f32::NAN.integer_decode();
assert_eq!((nan_m, nan_p), (12582912, 105));
}
#[test]
@ -28,6 +28,46 @@ fn test_f64_integer_decode() {
// Ignore the "sign" (quiet / signalling flag) of NAN.
// It can vary between runtime operations and LLVM folding.
let (nan_m, nan_e, _nan_s) = f64::NAN.integer_decode();
assert_eq!((nan_m, nan_e), (6755399441055744, 972));
let (nan_m, nan_p, _nan_s) = f64::NAN.integer_decode();
assert_eq!((nan_m, nan_p), (6755399441055744, 972));
}
/* Sanity checks of computed magic numbers */
#[test]
fn test_f32_consts() {
assert_eq!(<f32 as RawFloat>::INFINITY, f32::INFINITY);
assert_eq!(<f32 as RawFloat>::NEG_INFINITY, -f32::INFINITY);
assert_eq!(<f32 as RawFloat>::NAN.to_bits(), f32::NAN.to_bits());
assert_eq!(<f32 as RawFloat>::NEG_NAN.to_bits(), (-f32::NAN).to_bits());
assert_eq!(<f32 as RawFloat>::SIG_BITS, 23);
assert_eq!(<f32 as RawFloat>::MIN_EXPONENT_ROUND_TO_EVEN, -17);
assert_eq!(<f32 as RawFloat>::MAX_EXPONENT_ROUND_TO_EVEN, 10);
assert_eq!(<f32 as RawFloat>::MIN_EXPONENT_FAST_PATH, -10);
assert_eq!(<f32 as RawFloat>::MAX_EXPONENT_FAST_PATH, 10);
assert_eq!(<f32 as RawFloat>::MAX_EXPONENT_DISGUISED_FAST_PATH, 17);
assert_eq!(<f32 as RawFloat>::EXP_MIN, -126);
assert_eq!(<f32 as RawFloat>::EXP_SAT, 0xff);
assert_eq!(<f32 as RawFloat>::SMALLEST_POWER_OF_TEN, -65);
assert_eq!(<f32 as RawFloat>::LARGEST_POWER_OF_TEN, 38);
assert_eq!(<f32 as RawFloat>::MAX_MANTISSA_FAST_PATH, 16777216);
}
#[test]
fn test_f64_consts() {
assert_eq!(<f64 as RawFloat>::INFINITY, f64::INFINITY);
assert_eq!(<f64 as RawFloat>::NEG_INFINITY, -f64::INFINITY);
assert_eq!(<f64 as RawFloat>::NAN.to_bits(), f64::NAN.to_bits());
assert_eq!(<f64 as RawFloat>::NEG_NAN.to_bits(), (-f64::NAN).to_bits());
assert_eq!(<f64 as RawFloat>::SIG_BITS, 52);
assert_eq!(<f64 as RawFloat>::MIN_EXPONENT_ROUND_TO_EVEN, -4);
assert_eq!(<f64 as RawFloat>::MAX_EXPONENT_ROUND_TO_EVEN, 23);
assert_eq!(<f64 as RawFloat>::MIN_EXPONENT_FAST_PATH, -22);
assert_eq!(<f64 as RawFloat>::MAX_EXPONENT_FAST_PATH, 22);
assert_eq!(<f64 as RawFloat>::MAX_EXPONENT_DISGUISED_FAST_PATH, 37);
assert_eq!(<f64 as RawFloat>::EXP_MIN, -1022);
assert_eq!(<f64 as RawFloat>::EXP_SAT, 0x7ff);
assert_eq!(<f64 as RawFloat>::SMALLEST_POWER_OF_TEN, -342);
assert_eq!(<f64 as RawFloat>::LARGEST_POWER_OF_TEN, 308);
assert_eq!(<f64 as RawFloat>::MAX_MANTISSA_FAST_PATH, 9007199254740992);
}

View File

@ -1,13 +1,14 @@
use core::num::dec2flt::float::RawFloat;
use core::num::dec2flt::lemire::compute_float;
fn compute_float32(q: i64, w: u64) -> (i32, u64) {
let fp = compute_float::<f32>(q, w);
(fp.e, fp.f)
(fp.p_biased, fp.m)
}
fn compute_float64(q: i64, w: u64) -> (i32, u64) {
let fp = compute_float::<f64>(q, w);
(fp.e, fp.f)
(fp.p_biased, fp.m)
}
#[test]
@ -27,6 +28,11 @@ fn compute_float_f32_rounding() {
// Let's check the lines to see if anything is different in table...
assert_eq!(compute_float32(-10, 167772190000000000), (151, 2));
assert_eq!(compute_float32(-10, 167772200000000000), (151, 2));
// Check the rounding point between infinity and the next representable number down
assert_eq!(compute_float32(38, 3), (f32::INFINITE_POWER - 1, 6402534));
assert_eq!(compute_float32(38, 4), (f32::INFINITE_POWER, 0)); // infinity
assert_eq!(compute_float32(20, 3402823470000000000), (f32::INFINITE_POWER - 1, 8388607));
}
#[test]

View File

@ -1,5 +1,7 @@
#![allow(overflowing_literals)]
mod decimal;
mod decimal_seq;
mod float;
mod lemire;
mod parse;

View File

@ -1,9 +1,9 @@
use core::num::dec2flt::number::Number;
use core::num::dec2flt::decimal::Decimal;
use core::num::dec2flt::parse::parse_number;
use core::num::dec2flt::{dec2flt, pfe_invalid};
fn new_number(e: i64, m: u64) -> Number {
Number { exponent: e, mantissa: m, negative: false, many_digits: false }
fn new_dec(e: i64, m: u64) -> Decimal {
Decimal { exponent: e, mantissa: m, negative: false, many_digits: false }
}
#[test]
@ -31,23 +31,23 @@ fn invalid_chars() {
}
}
fn parse_positive(s: &[u8]) -> Option<Number> {
fn parse_positive(s: &[u8]) -> Option<Decimal> {
parse_number(s)
}
#[test]
fn valid() {
assert_eq!(parse_positive(b"123.456e789"), Some(new_number(786, 123456)));
assert_eq!(parse_positive(b"123.456e+789"), Some(new_number(786, 123456)));
assert_eq!(parse_positive(b"123.456e-789"), Some(new_number(-792, 123456)));
assert_eq!(parse_positive(b".050"), Some(new_number(-3, 50)));
assert_eq!(parse_positive(b"999"), Some(new_number(0, 999)));
assert_eq!(parse_positive(b"1.e300"), Some(new_number(300, 1)));
assert_eq!(parse_positive(b".1e300"), Some(new_number(299, 1)));
assert_eq!(parse_positive(b"101e-33"), Some(new_number(-33, 101)));
assert_eq!(parse_positive(b"123.456e789"), Some(new_dec(786, 123456)));
assert_eq!(parse_positive(b"123.456e+789"), Some(new_dec(786, 123456)));
assert_eq!(parse_positive(b"123.456e-789"), Some(new_dec(-792, 123456)));
assert_eq!(parse_positive(b".050"), Some(new_dec(-3, 50)));
assert_eq!(parse_positive(b"999"), Some(new_dec(0, 999)));
assert_eq!(parse_positive(b"1.e300"), Some(new_dec(300, 1)));
assert_eq!(parse_positive(b".1e300"), Some(new_dec(299, 1)));
assert_eq!(parse_positive(b"101e-33"), Some(new_dec(-33, 101)));
let zeros = "0".repeat(25);
let s = format!("1.5e{zeros}");
assert_eq!(parse_positive(s.as_bytes()), Some(new_number(-1, 15)));
assert_eq!(parse_positive(s.as_bytes()), Some(new_dec(-1, 15)));
}
macro_rules! assert_float_result_bits_eq {
@ -57,6 +57,21 @@ macro_rules! assert_float_result_bits_eq {
}};
}
#[test]
fn regression() {
// These showed up in fuzz tests when the minimum exponent was incorrect.
assert_float_result_bits_eq!(
0x0,
f64,
"3313756768023998018398807867233977556112078681253148176737587500333136120852692315608454494981109839693784033457129423181787087843504060087613228932431e-475"
);
assert_float_result_bits_eq!(
0x0,
f64,
"5298127456259331337220.92759278003098321644501973966679724599271041396379712108366679824365568578569680024083293475291869842408884554511641179110778276695274832779269225510492006696321279587846006535230380114430977056662212751544508159333199129106162019382177820713609e-346"
);
}
#[test]
fn issue31109() {
// Regression test for #31109.

View File

@ -147,12 +147,12 @@ pub trait Float:
}
macro_rules! impl_float {
($($fty:ty, $ity:ty, $bits:literal);+) => {
($($fty:ty, $ity:ty);+) => {
$(
impl Float for $fty {
type Int = $ity;
type SInt = <Self::Int as Int>::Signed;
const BITS: u32 = $bits;
const BITS: u32 = <$ity>::BITS;
const MAN_BITS: u32 = Self::MANTISSA_DIGITS - 1;
const MAN_MASK: Self::Int = (Self::Int::ONE << Self::MAN_BITS) - Self::Int::ONE;
const SIGN_MASK: Self::Int = Self::Int::ONE << (Self::BITS-1);
@ -168,7 +168,7 @@ macro_rules! impl_float {
}
}
impl_float!(f32, u32, 32; f64, u64, 64);
impl_float!(f32, u32; f64, u64);
/// A test generator. Should provide an iterator that produces unique patterns to parse.
///