From a0ed304c2101b4efdaf9982189e8d6ca73690712 Mon Sep 17 00:00:00 2001 From: Trevor Gross Date: Mon, 9 Dec 2024 09:25:21 +0000 Subject: float: Update some constants to `pub(crate)` These constants can be useful outside of their current module. Make them `pub(crate)` to allow for this. --- library/core/src/num/f32.rs | 6 +++--- library/core/src/num/f64.rs | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/library/core/src/num/f32.rs b/library/core/src/num/f32.rs index a200fd53186..de1557ccc90 100644 --- a/library/core/src/num/f32.rs +++ b/library/core/src/num/f32.rs @@ -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; diff --git a/library/core/src/num/f64.rs b/library/core/src/num/f64.rs index de63a462b61..65b5f3b9af0 100644 --- a/library/core/src/num/f64.rs +++ b/library/core/src/num/f64.rs @@ -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; -- cgit 1.4.1-3-g733a5 From 5a2da96a44abad6be752c7ee1cac173aa1ca2f7c Mon Sep 17 00:00:00 2001 From: Trevor Gross Date: Mon, 9 Dec 2024 09:25:22 +0000 Subject: dec2flt: Update documentation of existing methods Fix or elaborate existing float parsing documentation. This includes introducing a convention that should make naming more consistent. --- library/core/src/num/dec2flt/common.rs | 14 ++++++++------ library/core/src/num/dec2flt/decimal.rs | 22 +++++++++++++--------- library/core/src/num/dec2flt/mod.rs | 16 ++++++++++++++-- 3 files changed, 35 insertions(+), 17 deletions(-) diff --git a/library/core/src/num/dec2flt/common.rs b/library/core/src/num/dec2flt/common.rs index 4dadf406ae8..1646d3a95d3 100644 --- a/library/core/src/num/dec2flt/common.rs +++ b/library/core/src/num/dec2flt/common.rs @@ -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); diff --git a/library/core/src/num/dec2flt/decimal.rs b/library/core/src/num/dec2flt/decimal.rs index b37724ba62d..a84e11ec8e7 100644 --- a/library/core/src/num/dec2flt/decimal.rs +++ b/library/core/src/num/dec2flt/decimal.rs @@ -1,4 +1,4 @@ -//! Arbitrary-precision decimal class for fallback algorithms. +//! 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 @@ -11,6 +11,7 @@ use crate::num::dec2flt::common::{ByteSlice, is_8digits}; +/// A decimal floating-point number. #[derive(Clone)] pub(super) struct Decimal { /// The number of significant digits in the decimal. @@ -30,18 +31,17 @@ impl Default for Decimal { } impl Decimal { - /// The maximum number of digits required to unambiguously round a float. + /// The maximum number of digits required to unambiguously round up to a 64-bit float. /// - /// For a double-precision IEEE 754 float, this required 767 digits, - /// so we store the max digits + 1. + /// 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: + /// 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)⌋` /// @@ -56,11 +56,14 @@ impl Decimal { /// 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. + /// 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. + /// 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; @@ -69,6 +72,7 @@ impl Decimal { } /// Trim trailing zeros from the buffer. + // FIXME(tgross35): this could be `.rev().position()` if perf is okay pub(super) fn trim(&mut self) { // All of the following calls to `Decimal::trim` can't panic because: // @@ -86,7 +90,7 @@ impl Decimal { pub(super) fn round(&self) -> u64 { if self.num_digits == 0 || self.decimal_point < 0 { return 0; - } else if self.decimal_point > 18 { + } 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; diff --git a/library/core/src/num/dec2flt/mod.rs b/library/core/src/num/dec2flt/mod.rs index 6dca7406845..91bfe1bef2e 100644 --- a/library/core/src/num/dec2flt/mod.rs +++ b/library/core/src/num/dec2flt/mod.rs @@ -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( -- cgit 1.4.1-3-g733a5 From 49a2d4c757e27f915d15664ba92bbe0d3f590a8e Mon Sep 17 00:00:00 2001 From: Trevor Gross Date: Mon, 9 Dec 2024 09:25:22 +0000 Subject: dec2flt: Rename `Decimal` to `DecimalSeq` This module currently contains two decimal types, `Decimal` and `Number`. These names don't provide a whole lot of insight into what exactly they are, and `Number` is actually the one that is more like an expected `Decimal` type. In accordance with this, rename the existing `Decimal` to `DecimalSeq`. This highlights that it contains a sequence of decimal digits, rather than representing a base-10 floating point (decimal) number. Additionally, add some tests to validate internal behavior. --- library/core/src/num/dec2flt/decimal.rs | 362 -------------------- library/core/src/num/dec2flt/decimal_seq.rs | 379 +++++++++++++++++++++ library/core/src/num/dec2flt/mod.rs | 2 +- library/core/src/num/dec2flt/slow.rs | 8 +- library/coretests/tests/num/dec2flt/decimal_seq.rs | 30 ++ library/coretests/tests/num/dec2flt/mod.rs | 1 + 6 files changed, 415 insertions(+), 367 deletions(-) delete mode 100644 library/core/src/num/dec2flt/decimal.rs create mode 100644 library/core/src/num/dec2flt/decimal_seq.rs create mode 100644 library/coretests/tests/num/dec2flt/decimal_seq.rs diff --git a/library/core/src/num/dec2flt/decimal.rs b/library/core/src/num/dec2flt/decimal.rs deleted file mode 100644 index a84e11ec8e7..00000000000 --- a/library/core/src/num/dec2flt/decimal.rs +++ /dev/null @@ -1,362 +0,0 @@ -//! 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: . - -use crate::num::dec2flt::common::{ByteSlice, is_8digits}; - -/// A decimal floating-point number. -#[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], -} - -impl Default for Decimal { - fn default() -> Self { - Self { num_digits: 0, decimal_point: 0, truncated: false, digits: [0; Self::MAX_DIGITS] } - } -} - -impl Decimal { - /// 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(super) 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(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; - } - } - - 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(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; - } - } - num_new_digits -} diff --git a/library/core/src/num/dec2flt/decimal_seq.rs b/library/core/src/num/dec2flt/decimal_seq.rs new file mode 100644 index 00000000000..de22280c001 --- /dev/null +++ b/library/core/src/num/dec2flt/decimal_seq.rs @@ -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: . + +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 +} diff --git a/library/core/src/num/dec2flt/mod.rs b/library/core/src/num/dec2flt/mod.rs index 91bfe1bef2e..471505d2494 100644 --- a/library/core/src/num/dec2flt/mod.rs +++ b/library/core/src/num/dec2flt/mod.rs @@ -97,7 +97,7 @@ use crate::fmt; use crate::str::FromStr; mod common; -mod decimal; +pub mod decimal_seq; mod fpu; mod slow; mod table; diff --git a/library/core/src/num/dec2flt/slow.rs b/library/core/src/num/dec2flt/slow.rs index 85d4b13284b..e274e1aa0aa 100644 --- a/library/core/src/num/dec2flt/slow.rs +++ b/library/core/src/num/dec2flt/slow.rs @@ -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(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(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,7 +67,7 @@ pub(crate) fn parse_long_mantissa(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; diff --git a/library/coretests/tests/num/dec2flt/decimal_seq.rs b/library/coretests/tests/num/dec2flt/decimal_seq.rs new file mode 100644 index 00000000000..f46ba7c465a --- /dev/null +++ b/library/coretests/tests/num/dec2flt/decimal_seq.rs @@ -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); + } +} diff --git a/library/coretests/tests/num/dec2flt/mod.rs b/library/coretests/tests/num/dec2flt/mod.rs index 874e0ec7093..51f3017ddc7 100644 --- a/library/coretests/tests/num/dec2flt/mod.rs +++ b/library/coretests/tests/num/dec2flt/mod.rs @@ -1,5 +1,6 @@ #![allow(overflowing_literals)] +mod decimal_seq; mod float; mod lemire; mod parse; -- cgit 1.4.1-3-g733a5 From 626d2c5eed687e2c166d25f29cdfcaaba48f5f23 Mon Sep 17 00:00:00 2001 From: Trevor Gross Date: Mon, 9 Dec 2024 09:25:22 +0000 Subject: dec2flt: Rename `Number` to `Decimal` The previous commit renamed `Decimal` to `DecimalSeq`. Now, rename the type that represents a decimal floating point number to be `Decimal`. Additionally, add some tests for internal behavior. --- library/core/src/num/dec2flt/decimal.rs | 88 ++++++++++++++++++++++++++ library/core/src/num/dec2flt/mod.rs | 2 +- library/core/src/num/dec2flt/number.rs | 88 -------------------------- library/core/src/num/dec2flt/parse.rs | 10 +-- library/coretests/tests/num/dec2flt/decimal.rs | 28 ++++++++ library/coretests/tests/num/dec2flt/mod.rs | 1 + library/coretests/tests/num/dec2flt/parse.rs | 26 ++++---- 7 files changed, 136 insertions(+), 107 deletions(-) create mode 100644 library/core/src/num/dec2flt/decimal.rs delete mode 100644 library/core/src/num/dec2flt/number.rs create mode 100644 library/coretests/tests/num/dec2flt/decimal.rs diff --git a/library/core/src/num/dec2flt/decimal.rs b/library/core/src/num/dec2flt/decimal.rs new file mode 100644 index 00000000000..1ba9150a011 --- /dev/null +++ b/library/core/src/num/dec2flt/decimal.rs @@ -0,0 +1,88 @@ +//! Representation of a float as the significant digits and exponent. + +use crate::num::dec2flt::float::RawFloat; +use crate::num::dec2flt::fpu::set_precision; + +const INT_POW10: [u64; 16] = [ + 1, + 10, + 100, + 1000, + 10000, + 100000, + 1000000, + 10000000, + 100000000, + 1000000000, + 10000000000, + 100000000000, + 1000000000000, + 10000000000000, + 100000000000000, + 1000000000000000, +]; + +/// 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 { + /// Detect if the float can be accurately reconstructed from native floats. + #[inline] + fn is_fast_path(&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(&self) -> Option { + // Here we need to work around . + // 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::(); + + if self.is_fast_path::() { + 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 + } + } +} diff --git a/library/core/src/num/dec2flt/mod.rs b/library/core/src/num/dec2flt/mod.rs index 471505d2494..a82c858df49 100644 --- a/library/core/src/num/dec2flt/mod.rs +++ b/library/core/src/num/dec2flt/mod.rs @@ -97,6 +97,7 @@ use crate::fmt; use crate::str::FromStr; mod common; +pub mod decimal; pub mod decimal_seq; mod fpu; mod slow; @@ -104,7 +105,6 @@ 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 { diff --git a/library/core/src/num/dec2flt/number.rs b/library/core/src/num/dec2flt/number.rs deleted file mode 100644 index 2538991564a..00000000000 --- a/library/core/src/num/dec2flt/number.rs +++ /dev/null @@ -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(&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(&self) -> Option { - // Here we need to work around . - // 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::(); - - if self.is_fast_path::() { - 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 - } - } -} diff --git a/library/core/src/num/dec2flt/parse.rs b/library/core/src/num/dec2flt/parse.rs index 06ee8e95fbc..e38fedc58be 100644 --- a/library/core/src/num/dec2flt/parse.rs +++ b/library/core/src/num/dec2flt/parse.rs @@ -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 { /// /// 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 { +pub fn parse_number(s: &[u8]) -> Option { if let Some((float, rest)) = parse_partial_number(s) { if rest == s.len() { return Some(float); diff --git a/library/coretests/tests/num/dec2flt/decimal.rs b/library/coretests/tests/num/dec2flt/decimal.rs new file mode 100644 index 00000000000..1fa06de692e --- /dev/null +++ b/library/coretests/tests/num/dec2flt/decimal.rs @@ -0,0 +1,28 @@ +use core::num::dec2flt::decimal::Decimal; + +type FPath = ((i64, u64, bool, bool), Option); + +const FPATHS_F32: &[FPath] = + &[((0, 0, false, false), Some(0.0)), ((0, 0, false, false), Some(0.0))]; +const FPATHS_F64: &[FPath] = + &[((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::(); + + 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::(); + + assert_eq!(actual, expected); + } +} diff --git a/library/coretests/tests/num/dec2flt/mod.rs b/library/coretests/tests/num/dec2flt/mod.rs index 51f3017ddc7..a9025be5ca7 100644 --- a/library/coretests/tests/num/dec2flt/mod.rs +++ b/library/coretests/tests/num/dec2flt/mod.rs @@ -1,5 +1,6 @@ #![allow(overflowing_literals)] +mod decimal; mod decimal_seq; mod float; mod lemire; diff --git a/library/coretests/tests/num/dec2flt/parse.rs b/library/coretests/tests/num/dec2flt/parse.rs index 4a5d24ba7d5..ec705a61e13 100644 --- a/library/coretests/tests/num/dec2flt/parse.rs +++ b/library/coretests/tests/num/dec2flt/parse.rs @@ -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 { +fn parse_positive(s: &[u8]) -> Option { 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 { -- cgit 1.4.1-3-g733a5 From 6c34daff57101922dc97e1860bf5940ea88d3906 Mon Sep 17 00:00:00 2001 From: Trevor Gross Date: Mon, 9 Dec 2024 09:25:22 +0000 Subject: dec2flt: Rename fields to be consistent with documented notation --- library/core/src/num/dec2flt/common.rs | 13 +++++++------ library/core/src/num/dec2flt/lemire.rs | 4 ++-- library/core/src/num/dec2flt/mod.rs | 13 ++++++++----- library/core/src/num/dec2flt/slow.rs | 2 +- library/coretests/tests/num/dec2flt/float.rs | 8 ++++---- library/coretests/tests/num/dec2flt/lemire.rs | 4 ++-- 6 files changed, 24 insertions(+), 20 deletions(-) diff --git a/library/core/src/num/dec2flt/common.rs b/library/core/src/num/dec2flt/common.rs index 1646d3a95d3..a140a311c45 100644 --- a/library/core/src/num/dec2flt/common.rs +++ b/library/core/src/num/dec2flt/common.rs @@ -63,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 } } } diff --git a/library/core/src/num/dec2flt/lemire.rs b/library/core/src/num/dec2flt/lemire.rs index 01642e1b111..0cc39cb9b62 100644 --- a/library/core/src/num/dec2flt/lemire.rs +++ b/library/core/src/num/dec2flt/lemire.rs @@ -73,7 +73,7 @@ pub fn compute_float(q: i64, mut w: u64) -> BiasedFp { mantissa += mantissa & 1; mantissa >>= 1; power2 = (mantissa >= (1_u64 << F::MANTISSA_EXPLICIT_BITS)) as i32; - return BiasedFp { f: mantissa, e: power2 }; + 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 @@ -111,7 +111,7 @@ pub fn compute_float(q: i64, mut w: u64) -> BiasedFp { // 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. diff --git a/library/core/src/num/dec2flt/mod.rs b/library/core/src/num/dec2flt/mod.rs index a82c858df49..9d1989c4b81 100644 --- a/library/core/src/num/dec2flt/mod.rs +++ b/library/core/src/num/dec2flt/mod.rs @@ -233,8 +233,8 @@ pub fn pfe_invalid() -> ParseFloatError { /// Converts a `BiasedFp` to the closest machine float type. fn biased_fp_to_float(x: BiasedFp) -> T { - let mut word = x.f; - word |= (x.e as u64) << T::MANTISSA_EXPLICIT_BITS; + let mut word = x.m; + word |= (x.p_biased as u64) << T::MANTISSA_EXPLICIT_BITS; T::from_u64_bits(word) } @@ -272,12 +272,15 @@ pub fn dec2flt(s: &str) -> Result { // redundantly using the Eisel-Lemire algorithm if it was unable to // correctly round on the first pass. let mut fp = compute_float::(num.exponent, num.mantissa); - if num.many_digits && fp.e >= 0 && fp != compute_float::(num.exponent, num.mantissa + 1) { - fp.e = -1; + if num.many_digits + && fp.p_biased >= 0 + && fp != compute_float::(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::(s); } diff --git a/library/core/src/num/dec2flt/slow.rs b/library/core/src/num/dec2flt/slow.rs index e274e1aa0aa..e0c4ae117da 100644 --- a/library/core/src/num/dec2flt/slow.rs +++ b/library/core/src/num/dec2flt/slow.rs @@ -105,5 +105,5 @@ pub(crate) fn parse_long_mantissa(s: &[u8]) -> BiasedFp { } // Zero out all the bits above the explicit mantissa bits. mantissa &= (1_u64 << F::MANTISSA_EXPLICIT_BITS) - 1; - BiasedFp { f: mantissa, e: power2 } + BiasedFp { m: mantissa, p_biased: power2 } } diff --git a/library/coretests/tests/num/dec2flt/float.rs b/library/coretests/tests/num/dec2flt/float.rs index 7a9587a18d0..40f0776e021 100644 --- a/library/coretests/tests/num/dec2flt/float.rs +++ b/library/coretests/tests/num/dec2flt/float.rs @@ -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,6 @@ 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)); } diff --git a/library/coretests/tests/num/dec2flt/lemire.rs b/library/coretests/tests/num/dec2flt/lemire.rs index f71bbb7c7a3..9f228a25e46 100644 --- a/library/coretests/tests/num/dec2flt/lemire.rs +++ b/library/coretests/tests/num/dec2flt/lemire.rs @@ -2,12 +2,12 @@ use core::num::dec2flt::lemire::compute_float; fn compute_float32(q: i64, w: u64) -> (i32, u64) { let fp = compute_float::(q, w); - (fp.e, fp.f) + (fp.p_biased, fp.m) } fn compute_float64(q: i64, w: u64) -> (i32, u64) { let fp = compute_float::(q, w); - (fp.e, fp.f) + (fp.p_biased, fp.m) } #[test] -- cgit 1.4.1-3-g733a5 From 19a909ae0e9f923db2da0fc10264a75b7aa5ad99 Mon Sep 17 00:00:00 2001 From: Trevor Gross Date: Mon, 9 Dec 2024 09:25:22 +0000 Subject: dec2flt: Refactor float traits A lot of the magic constants can be turned into expressions. This reduces some code duplication. Additionally, add traits to make these operations fully generic. This will make it easier to support `f16` and `f128`. --- library/core/src/num/dec2flt/float.rs | 272 +++++++++++++++++--------- library/core/src/num/dec2flt/lemire.rs | 18 +- library/core/src/num/dec2flt/mod.rs | 6 +- library/core/src/num/dec2flt/slow.rs | 18 +- library/coretests/tests/num/dec2flt/float.rs | 40 ++++ library/coretests/tests/num/dec2flt/lemire.rs | 6 + library/coretests/tests/num/dec2flt/parse.rs | 15 ++ src/etc/test-float-parse/src/traits.rs | 6 +- 8 files changed, 267 insertions(+), 114 deletions(-) diff --git a/library/core/src/num/dec2flt/float.rs b/library/core/src/num/dec2flt/float.rs index da57aa9a546..b8a28a67569 100644 --- a/library/core/src/num/dec2flt/float.rs +++ b/library/core/src/num/dec2flt/float.rs @@ -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}; + +/// Lossy `as` casting between two types. +pub trait CastInto: 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 + + ops::Shl + + ops::BitAnd + + ops::BitOr + + PartialEq + + CastInto +{ + const ZERO: Self; + const ONE: Self; +} + +macro_rules! int { + ($($ty:ty),+) => { + $( + impl CastInto 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 `f32` and `f64`. +/// 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; + + /* 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; - - // 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; + /// Bit width of the float + const BITS: u32; - // Minimum exponent that for a fast path case, or `-⌊(MANTISSA_EXPLICIT_BITS+1)/log2(5)⌋` - const MIN_EXPONENT_FAST_PATH: i64; + /// The number of bits in the significand, *including* the hidden bit. + const SIG_TOTAL_BITS: u32; - // Maximum exponent that for a fast path case, or `⌊(MANTISSA_EXPLICIT_BITS+1)/log2(5)⌋` - const MAX_EXPONENT_FAST_PATH: i64; + const EXP_MASK: Self::Int; + const SIG_MASK: Self::Int; - // 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; + /// 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; - // Largest exponent value `(1 << EXP_BITS) - 1`. - const INFINITE_POWER: i32; + /* limits related to Fast pathing */ + + /// 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 + }; + + /// 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; - // Index (in bits) of the sign. - const SIGN_INDEX: usize; + /// 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 + }; - // Smallest decimal exponent for a non-zero value. - const SMALLEST_POWER_OF_TEN: i32; + /// Minimum exponent for a fast path case, or `-⌊(SIG_BITS+1)/log2(5)⌋` + const MIN_EXPONENT_FAST_PATH: i64 = -Self::MAX_EXPONENT_FAST_PATH; - // Largest decimal exponent for a non-infinite value. - const LARGEST_POWER_OF_TEN: i32; + /// 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 = 2_u64 << Self::MANTISSA_EXPLICIT_BITS; + /// 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 { diff --git a/library/core/src/num/dec2flt/lemire.rs b/library/core/src/num/dec2flt/lemire.rs index 0cc39cb9b62..f84929a03c1 100644 --- a/library/core/src/num/dec2flt/lemire.rs +++ b/library/core/src/num/dec2flt/lemire.rs @@ -38,7 +38,7 @@ pub fn compute_float(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(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,7 +72,7 @@ pub fn compute_float(q: i64, mut w: u64) -> BiasedFp { mantissa >>= -power2 + 1; mantissa += mantissa & 1; mantissa >>= 1; - power2 = (mantissa >= (1_u64 << F::MANTISSA_EXPLICIT_BITS)) as i32; + 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, @@ -89,8 +89,8 @@ pub fn compute_float(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,15 +98,15 @@ pub fn compute_float(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; diff --git a/library/core/src/num/dec2flt/mod.rs b/library/core/src/num/dec2flt/mod.rs index 9d1989c4b81..d1a0e1db313 100644 --- a/library/core/src/num/dec2flt/mod.rs +++ b/library/core/src/num/dec2flt/mod.rs @@ -232,10 +232,10 @@ pub fn pfe_invalid() -> ParseFloatError { } /// Converts a `BiasedFp` to the closest machine float type. -fn biased_fp_to_float(x: BiasedFp) -> T { +fn biased_fp_to_float(x: BiasedFp) -> F { let mut word = x.m; - word |= (x.p_biased as u64) << T::MANTISSA_EXPLICIT_BITS; - T::from_u64_bits(word) + word |= (x.p_biased as u64) << F::SIG_BITS; + F::from_u64_bits(word) } /// Converts a decimal string into a floating point number. diff --git a/library/core/src/num/dec2flt/slow.rs b/library/core/src/num/dec2flt/slow.rs index e0c4ae117da..3baed426523 100644 --- a/library/core/src/num/dec2flt/slow.rs +++ b/library/core/src/num/dec2flt/slow.rs @@ -74,36 +74,36 @@ pub(crate) fn parse_long_mantissa(s: &[u8]) -> BiasedFp { } // 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; + mantissa &= (1_u64 << F::SIG_BITS) - 1; BiasedFp { m: mantissa, p_biased: power2 } } diff --git a/library/coretests/tests/num/dec2flt/float.rs b/library/coretests/tests/num/dec2flt/float.rs index 40f0776e021..b5afd3e3b24 100644 --- a/library/coretests/tests/num/dec2flt/float.rs +++ b/library/coretests/tests/num/dec2flt/float.rs @@ -31,3 +31,43 @@ fn test_f64_integer_decode() { 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!(::INFINITY, f32::INFINITY); + assert_eq!(::NEG_INFINITY, -f32::INFINITY); + assert_eq!(::NAN.to_bits(), f32::NAN.to_bits()); + assert_eq!(::NEG_NAN.to_bits(), (-f32::NAN).to_bits()); + assert_eq!(::SIG_BITS, 23); + assert_eq!(::MIN_EXPONENT_ROUND_TO_EVEN, -17); + assert_eq!(::MAX_EXPONENT_ROUND_TO_EVEN, 10); + assert_eq!(::MIN_EXPONENT_FAST_PATH, -10); + assert_eq!(::MAX_EXPONENT_FAST_PATH, 10); + assert_eq!(::MAX_EXPONENT_DISGUISED_FAST_PATH, 17); + assert_eq!(::EXP_MIN, -126); + assert_eq!(::EXP_SAT, 0xff); + assert_eq!(::SMALLEST_POWER_OF_TEN, -65); + assert_eq!(::LARGEST_POWER_OF_TEN, 38); + assert_eq!(::MAX_MANTISSA_FAST_PATH, 16777216); +} + +#[test] +fn test_f64_consts() { + assert_eq!(::INFINITY, f64::INFINITY); + assert_eq!(::NEG_INFINITY, -f64::INFINITY); + assert_eq!(::NAN.to_bits(), f64::NAN.to_bits()); + assert_eq!(::NEG_NAN.to_bits(), (-f64::NAN).to_bits()); + assert_eq!(::SIG_BITS, 52); + assert_eq!(::MIN_EXPONENT_ROUND_TO_EVEN, -4); + assert_eq!(::MAX_EXPONENT_ROUND_TO_EVEN, 23); + assert_eq!(::MIN_EXPONENT_FAST_PATH, -22); + assert_eq!(::MAX_EXPONENT_FAST_PATH, 22); + assert_eq!(::MAX_EXPONENT_DISGUISED_FAST_PATH, 37); + assert_eq!(::EXP_MIN, -1022); + assert_eq!(::EXP_SAT, 0x7ff); + assert_eq!(::SMALLEST_POWER_OF_TEN, -342); + assert_eq!(::LARGEST_POWER_OF_TEN, 308); + assert_eq!(::MAX_MANTISSA_FAST_PATH, 9007199254740992); +} diff --git a/library/coretests/tests/num/dec2flt/lemire.rs b/library/coretests/tests/num/dec2flt/lemire.rs index 9f228a25e46..0db80fbd525 100644 --- a/library/coretests/tests/num/dec2flt/lemire.rs +++ b/library/coretests/tests/num/dec2flt/lemire.rs @@ -1,3 +1,4 @@ +use core::num::dec2flt::float::RawFloat; use core::num::dec2flt::lemire::compute_float; fn compute_float32(q: i64, w: u64) -> (i32, u64) { @@ -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] diff --git a/library/coretests/tests/num/dec2flt/parse.rs b/library/coretests/tests/num/dec2flt/parse.rs index ec705a61e13..59be3915052 100644 --- a/library/coretests/tests/num/dec2flt/parse.rs +++ b/library/coretests/tests/num/dec2flt/parse.rs @@ -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. diff --git a/src/etc/test-float-parse/src/traits.rs b/src/etc/test-float-parse/src/traits.rs index f5333d63b36..f7f5f5dd35b 100644 --- a/src/etc/test-float-parse/src/traits.rs +++ b/src/etc/test-float-parse/src/traits.rs @@ -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 = ::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. /// -- cgit 1.4.1-3-g733a5 From 37e223ccaa1d0036540e6e3de5d4ed4b721afd60 Mon Sep 17 00:00:00 2001 From: Trevor Gross Date: Mon, 9 Dec 2024 09:25:22 +0000 Subject: dec2flt: Refactor the fast path This is just a bit of code cleanup to make use of returning early. --- library/core/src/num/dec2flt/decimal.rs | 47 ++++++++++++++++----------------- 1 file changed, 23 insertions(+), 24 deletions(-) diff --git a/library/core/src/num/dec2flt/decimal.rs b/library/core/src/num/dec2flt/decimal.rs index 1ba9150a011..db7176c1243 100644 --- a/library/core/src/num/dec2flt/decimal.rs +++ b/library/core/src/num/dec2flt/decimal.rs @@ -34,14 +34,15 @@ pub struct Decimal { impl Decimal { /// Detect if the float can be accurately reconstructed from native floats. #[inline] - fn is_fast_path(&self) -> bool { + fn can_use_fast_path(&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. + /// 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 @@ -59,30 +60,28 @@ impl Decimal { // require setting it by changing the global state (like the control word of the x87 FPU). let _cw = set_precision::(); - if self.is_fast_path::() { - 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 _) - } + if !self.can_use_fast_path::() { + return None; + } + + 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 { - // 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; + value * F::pow10_fast_path(self.exponent as _) } - Some(value) } else { - None - } + // 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) } } } -- cgit 1.4.1-3-g733a5