user0/rust
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Accurate decimal-to-float parsing routines.

Browse source 
This commit primarily adds implementations of the algorithms from William
Clinger's paper "How to Read Floating Point Numbers Accurately". It also
includes a lot of infrastructure necessary for those algorithms, and some
unit tests.

Since these algorithms reject a few (extreme) inputs that were previously
accepted, this could be seen as a [breaking-change]
This commit is contained in:
Robin Kruppe 2015-07-26 17:50:29 +02:00
parent b7e39a1c2d
commit ba792a4baa
13 changed files with 2787 additions and 15 deletions
Download patch file Download diff file Expand all files Collapse all files

353
src/libcore/num/dec2flt/algorithm.rs Normal file
View file

@ -0,0 +1,353 @@
// Copyright 2015 The Rust Project Developers. See the COPYRIGHT
// file at the top-level directory of this distribution and at
// http://rust-lang.org/COPYRIGHT.
//
// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
// http://www.apache.org/licenses/LICENSE-2.0> or the MIT license
// <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your
// option. This file may not be copied, modified, or distributed
// except according to those terms.
//! The various algorithms from the paper.
use num::flt2dec::strategy::grisu::Fp;
use prelude::v1::*;
use cmp::min;
use cmp::Ordering::{Less, Equal, Greater};
use super::table;
use super::rawfp::{self, Unpacked, RawFloat, fp_to_float, next_float, prev_float};
use super::num::{self, Big};
/// Number of significand bits in Fp
const P: u32 = 64;
// We simply store the best approximation for *all* exponents, so
// the variable "h" and the associated conditions can be omitted.
// This trades performance for space (11 KiB versus... 5 KiB or so?)
fn power_of_ten(e: i16) -> Fp {
assert!(e >= table::MIN_E);
let i = e - table::MIN_E;
let sig = table::POWERS.0[i as usize];
let exp = table::POWERS.1[i as usize];
Fp { f: sig, e: exp }
}
/// The fast path of Bellerophon using machine-sized integers and floats.
///
/// This is extracted into a separate function so that it can be attempted before constructing
/// a bignum.
pub fn fast_path<T: RawFloat>(integral: &[u8], fractional: &[u8], e: i64) -> Option<T> {
let num_digits = integral.len() + fractional.len();
// log_10(f64::max_sig) ~ 15.95. We compare the exact value to max_sig near the end,
// this is just a quick, cheap rejection (and also frees the rest of the code from
// worrying about underflow).
if num_digits > 16 {
return None;
}
if e.abs() >= T::ceil_log5_of_max_sig() as i64 {
return None;
}
let f = num::from_str_unchecked(integral.iter().chain(fractional.iter()));
if f > T::max_sig() {
return None;
}
let e = e as i16; // Can't overflow because e.abs() <= LOG5_OF_EXP_N
// The case e < 0 cannot be folded into the other branch. Negative powers result in
// a repeating fractional part in binary, which are rounded, which causes real
// (and occasioally quite significant!) errors in the final result.
// The case `e == 0`, however, is unnecessary for correctness. It's just measurably faster.
if e == 0 {
Some(T::from_int(f))
} else if e > 0 {
Some(T::from_int(f) * fp_to_float(power_of_ten(e)))
} else {
Some(T::from_int(f) / fp_to_float(power_of_ten(-e)))
}
}
/// Algorithm Bellerophon is trivial code justified by non-trivial numeric analysis.
///
/// It rounds ``f`` to a float with 64 bit significand and multiplies it by the best approximation
/// of `10^e` (in the same floating point format). This is often enough to get the correct result.
/// However, when the result is close to halfway between two adjecent (ordinary) floats, the
/// compound rounding error from multiplying two approximation means the result may be off by a
/// few bits. When this happens, the iterative Algorithm R fixes things up.
///
/// The hand-wavy "close to halfway" is made precise by the numeric analysis in the paper.
/// In the words of Clinger:
///
/// > Slop, expressed in units of the least significant bit, is an inclusive bound for the error
/// > accumulated during the floating point calculation of the approximation to f * 10^e. (Slop is
/// > not a bound for the true error, but bounds the difference between the approximation z and
/// > the best possible approximation that uses p bits of significand.)
pub fn bellerophon<T: RawFloat>(f: &Big, e: i16) -> T {
let slop;
if f <= &Big::from_u64(T::max_sig()) {
// The cases abs(e) < log5(2^N) are in fast_path()
slop = if e >= 0 { 0 } else { 3 };
} else {
slop = if e >= 0 { 1 } else { 4 };
}
let z = rawfp::big_to_fp(f).mul(&power_of_ten(e)).normalize();
let exp_p_n = 1 << (P - T::sig_bits() as u32);
let lowbits: i64 = (z.f % exp_p_n) as i64;
// Is the slop large enough to make a difference when
// rounding to n bits?
if (lowbits - exp_p_n as i64 / 2).abs() <= slop {
algorithm_r(f, e, fp_to_float(z))
} else {
fp_to_float(z)
}
}
/// An iterative algorithm that improves a floating point approximation of `f * 10^e`.
///
/// Each iteration gets one unit in the last place closer, which of course takes terribly long to
/// converge if `z0` is even mildly off. Luckily, when used as fallback for Bellerophon, the
/// starting approximation is off by at most one ULP.
fn algorithm_r<T: RawFloat>(f: &Big, e: i16, z0: T) -> T {
let mut z = z0;
loop {
let raw = z.unpack();
let (m, k) = (raw.sig, raw.k);
let mut x = f.clone();
let mut y = Big::from_u64(m);
// Find positive integers `x`, `y` such that `x / y` is exactly `(f * 10^e) / (m * 2^k)`.
// This not only avoids dealing with the signs of `e` and `k`, we also eliminate the
// power of two common to `10^e` and `2^k` to make the numbers smaller.
make_ratio(&mut x, &mut y, e, k);
let m_digits = [(m & 0xFF_FF_FF_FF) as u32, (m >> 32) as u32];
// This is written a bit awkwardly because our bignums don't support
// negative numbers, so we use the absolute value + sign information.
// The multiplication with m_digits can't overflow. If `x` or `y` are large enough that
// we need to worry about overflow, then they are also large enough that`make_ratio` has
// reduced the fraction by a factor of 2^64 or more.
let (d2, d_negative) = if x >= y {
// Don't need x any more, save a clone().
x.sub(&y).mul_pow2(1).mul_digits(&m_digits);
(x, false)
} else {
// Still need y - make a copy.
let mut y = y.clone();
y.sub(&x).mul_pow2(1).mul_digits(&m_digits);
(y, true)
};
if d2 < y {
let mut d2_double = d2;
d2_double.mul_pow2(1);
if m == T::min_sig() && d_negative && d2_double > y {
z = prev_float(z);
} else {
return z;
}
} else if d2 == y {
if m % 2 == 0 {
if m == T::min_sig() && d_negative {
z = prev_float(z);
} else {
return z;
}
} else if d_negative {
z = prev_float(z);
} else {
z = next_float(z);
}
} else if d_negative {
z = prev_float(z);
} else {
z = next_float(z);
}
}
}
/// Given `x = f` and `y = m` where `f` represent input decimal digits as usual and `m` is the
/// significand of a floating point approximation, make the ratio `x / y` equal to
/// `(f * 10^e) / (m * 2^k)`, possibly reduced by a power of two both have in common.
fn make_ratio(x: &mut Big, y: &mut Big, e: i16, k: i16) {
let (e_abs, k_abs) = (e.abs() as usize, k.abs() as usize);
if e >= 0 {
if k >= 0 {
// x = f * 10^e, y = m * 2^k, except that we reduce the fraction by some power of two.
let common = min(e_abs, k_abs);
x.mul_pow5(e_abs).mul_pow2(e_abs - common);
y.mul_pow2(k_abs - common);
} else {
// x = f * 10^e * 2^abs(k), y = m
// This can't overflow because it requires positive `e` and negative `k`, which can
// only happen for values extremely close to 1, which means that `e` and `k` will be
// comparatively tiny.
x.mul_pow5(e_abs).mul_pow2(e_abs + k_abs);
}
} else {
if k >= 0 {
// x = f, y = m * 10^abs(e) * 2^k
// This can't overflow either, see above.
y.mul_pow5(e_abs).mul_pow2(k_abs + e_abs);
} else {
// x = f * 2^abs(k), y = m * 10^abs(e), again reducing by a common power of two.
let common = min(e_abs, k_abs);
x.mul_pow2(k_abs - common);
y.mul_pow5(e_abs).mul_pow2(e_abs - common);
}
}
}
/// Conceptually, Algorithm M is the simplest way to convert a decimal to a float.
///
/// We form a ratio that is equal to `f * 10^e`, then throwing in powers of two until it gives
/// a valid float significand. The binary exponent `k` is the number of times we multiplied
/// numerator or denominator by two, i.e., at all times `f * 10^e` equals `(u / v) * 2^k`.
/// When we have found out significand, we only need to round by inspecting the remainder of the
/// division, which is done in helper functions further below.
///
/// This algorithm is super slow, even with the optimization described in `quick_start()`.
/// However, it's the simplest of the algorithms to adapt for overflow, underflow, and subnormal
/// results. This implementation takes over when Bellerophon and Algorithm R are overwhelmed.
/// Detecting underflow and overflow is easy: The ratio still isn't an in-range significand,
/// yet the minimum/maximum exponent has been reached. In the case of overflow, we simply return
/// infinity.
///
/// Handling underflow and subnormals is trickier. One big problem is that, with the minimum
/// exponent, the ratio might still be too large for a significand. See underflow() for details.
pub fn algorithm_m<T: RawFloat>(f: &Big, e: i16) -> T {
let mut u;
let mut v;
let e_abs = e.abs() as usize;
let mut k = 0;
if e < 0 {
u = f.clone();
v = Big::from_small(1);
v.mul_pow5(e_abs).mul_pow2(e_abs);
} else {
// FIXME possible optimization: generalize big_to_fp so that we can do the equivalent of
// fp_to_float(big_to_fp(u)) here, only without the double rounding.
u = f.clone();
u.mul_pow5(e_abs).mul_pow2(e_abs);
v = Big::from_small(1);
}
quick_start::<T>(&mut u, &mut v, &mut k);
let mut rem = Big::from_small(0);
let mut x = Big::from_small(0);
let min_sig = Big::from_u64(T::min_sig());
let max_sig = Big::from_u64(T::max_sig());
loop {
u.div_rem(&v, &mut x, &mut rem);
if k == T::min_exp_int() {
// We have to stop at the minimum exponent, if we wait until `k < T::min_exp_int()`,
// then we'd be off by a factor of two. Unfortunately this means we have to special-
// case normal numbers with the minimum exponent.
// FIXME find a more elegant formulation, but run the `tiny-pow10` test to make sure
// that it's actually correct!
if x >= min_sig && x <= max_sig {
break;
}
return underflow(x, v, rem);
}
if k > T::max_exp_int() {
return T::infinity();
}
if x < min_sig {
u.mul_pow2(1);
k -= 1;
} else if x > max_sig {
v.mul_pow2(1);
k += 1;
} else {
break;
}
}
let q = num::to_u64(&x);
let z = rawfp::encode_normal(Unpacked::new(q, k));
round_by_remainder(v, rem, q, z)
}
/// Skip over most AlgorithmM iterations by checking the bit length.
fn quick_start<T: RawFloat>(u: &mut Big, v: &mut Big, k: &mut i16) {
// The bit length is an estimate of the base two logarithm, and log(u / v) = log(u) - log(v).
// The estimate is off by at most 1, but always an under-estimate, so the error on log(u)
// and log(v) are of the same sign and cancel out (if both are large). Therefore the error
// for log(u / v) is at most one as well.
// The target ratio is one where u/v is in an in-range significand. Thus our termination
// condition is log2(u / v) being the significand bits, plus/minus one.
// FIXME Looking at the second bit could improve the estimate and avoid some more divisions.
let target_ratio = f64::sig_bits() as i16;
let log2_u = u.bit_length() as i16;
let log2_v = v.bit_length() as i16;
let mut u_shift: i16 = 0;
let mut v_shift: i16 = 0;
assert!(*k == 0);
loop {
if *k == T::min_exp_int() {
// Underflow or subnormal. Leave it to the main function.
break;
}
if *k == T::max_exp_int() {
// Overflow. Leave it to the main function.
break;
}
let log2_ratio = (log2_u + u_shift) - (log2_v + v_shift);
if log2_ratio < target_ratio - 1 {
u_shift += 1;
*k -= 1;
} else if log2_ratio > target_ratio + 1 {
v_shift += 1;
*k += 1;
} else {
break;
}
}
u.mul_pow2(u_shift as usize);
v.mul_pow2(v_shift as usize);
}
fn underflow<T: RawFloat>(x: Big, v: Big, rem: Big) -> T {
if x < Big::from_u64(T::min_sig()) {
let q = num::to_u64(&x);
let z = rawfp::encode_subnormal(q);
return round_by_remainder(v, rem, q, z);
}
// Ratio isn't an in-range significand with the minimum exponent, so we need to round off
// excess bits and adjust the exponent accordingly. The real value now looks like this:
//
// x lsb
// /--------------\/
// 1010101010101010.10101010101010 * 2^k
// \-----/\-------/ \------------/
// q trunc. (represented by rem)
//
// Therefore, when the rounded-off bits are != 0.5 ULP, they decide the rounding
// on their own. When they are equal and the remainder is non-zero, the value still
// needs to be rounded up. Only when the rounded off bits are 1/2 and the remainer
// is zero, we have a half-to-even situation.
let bits = x.bit_length();
let lsb = bits - T::sig_bits() as usize;
let q = num::get_bits(&x, lsb, bits);
let k = T::min_exp_int() + lsb as i16;
let z = rawfp::encode_normal(Unpacked::new(q, k));
let q_even = q % 2 == 0;
match num::compare_with_half_ulp(&x, lsb) {
Greater => next_float(z),
Less => z,
Equal if rem.is_zero() && q_even => z,
Equal => next_float(z),
}
}
/// Ordinary round-to-even, obfuscated by having to round based on the remainder of a division.
fn round_by_remainder<T: RawFloat>(v: Big, r: Big, q: u64, z: T) -> T {
let mut v_minus_r = v;
v_minus_r.sub(&r);
if r < v_minus_r {
z
} else if r > v_minus_r {
next_float(z)
} else if q % 2 == 0 {
z
} else {
next_float(z)
}
}

234
src/libcore/num/dec2flt/mod.rs Normal file
View file

@ -0,0 +1,234 @@
// Copyright 2015 The Rust Project Developers. See the COPYRIGHT
// file at the top-level directory of this distribution and at
// http://rust-lang.org/COPYRIGHT.
//
// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
// http://www.apache.org/licenses/LICENSE-2.0> or the MIT license
// <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your
// option. This file may not be copied, modified, or distributed
// except according to those terms.
//! Converting decimal strings into IEEE 754 binary floating point numbers.
//!
//! # Problem statement
//!
//! We are given a decimal string such as `12.34e56`. This string consists of integral (`12`),
//! fractional (`45`), and exponent (`56`) parts. All parts are optional and interpreted as zero
//! 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
//! base two, so we round to 0.5 units in the last place (in other words, as well as possible).
//! Ties, decimal values exactly half-way between two consecutive floats, are resolved with the
//! half-to-even strategy, also known as banker's rounding.
//!
//! Needless to say, this is quite hard, both in terms of implementation complexity and in terms
//! of CPU cycles taken.
//!
//! # Implementation
//!
//! First, we ignore signs. Or rather, we remove it at the very beginning of the conversion
//! process and re-apply it at the very end. This is correct in all edge cases since IEEE
//! floats are symmetric around zero, negating one simply flips the first bit.
//!
//! Then we remove the decimal point by adjusting the exponent: Conceptually, `12.34e56` turns
//! into `1234e54`, which we describe with a positive integer `f = 1234` and an integer `e = 54`.
//! The `(f, e)` representation is used by almost all code past the parsing stage.
//!
//! We then try a long chain of progressively more general and expensive special cases using
//! machine-sized integers and small, fixed-sized floating point numbers (first `f32`/`f64`, then
//! a type with 64 bit significand, `Fp`). When all these fail, we bite the bullet and resort to a
//! simple but very slow algorithm that involved computing `f * 10^e` fully and doing an iterative
//! search for the best approximation.
//!
//! Primarily, this module and its children implement the algorithms described in:
//! "How to Read Floating Point Numbers Accurately" by William D. Clinger,
//! available online: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.45.4152
//!
//! In addition, there are numerous helper functions that are used in the paper but not available
//! in Rust (or at least in core). Our version is additionally complicated by the need to handle
//! overflow and underflow and the desire to handle subnormal numbers. Bellerophon and
//! Algorithm R have trouble with overflow, subnormals, and underflow. We conservatively switch to
//! Algorithm M (with the modifications described in section 8 of the paper) well before the
//! inputs get into the critical region.
//!
//! Another aspect that needs attention is the ``RawFloat`` trait by which almost all functions
//! are parametrized. One might think that it's enough to parse to `f64` and cast the result to
//! `f32`. Unfortunately this is not the world we live in, and this has nothing to do with using
//! base two or half-to-even rounding.
//!
//! Consider for example two types `d2` and `d4` representing a decimal type with two decimal
//! digits and four decimal digits each and take "0.01499" as input. Let's use half-up rounding.
//! Going directly to two decimal digits gives `0.01`, but if we round to four digits first,
//! we get `0.0150`, which is then rounded up to `0.02`. The same principle applies to other
//! operations as well, if you want 0.5 ULP accuracy you need to do *everything* in full precision
//! and round *exactly once, at the end*, by considering all truncated bits at once.
//!
//! FIXME Although some code duplication is necessary, perhaps parts of the code could be shuffled
//! around such that less code is duplicated. Large parts of the algorithms are independent of the
//! float type to output, or only needs access to a few constants, which could be passed in as
//! parameters.
//!
//! # Other
//!
//! The conversion should *never* panic. There are assertions and explicit panics in the code,
//! but they should never be triggered and only serve as internal sanity checks. Any panics should
//! be considered a bug.
//!
//! There are unit tests but they are woefully inadequate at ensuring correctness, they only cover
//! a small percentage of possible errors. Far more extensive tests are located in the directory
//! `src/etc/test-float-parse` as a Python script.
//!
//! A note on integer overflow: Many parts of this file perform arithmetic with the decimal
//! exponent `e`. Primarily, we shift the decimal point around: Before the first decimal digit,
//! after the last decimal digit, and so on. This could overflow if done carelessly. We rely on
//! the parsing submodule to only hand out sufficiently small exponents, where "sufficient" means
//! "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}.
//!
//! FIXME: this uses several things from core::num::flt2dec, which is nonsense. Those things
//! should be moved into core::num::<something else>.
#![doc(hidden)]
#![unstable(feature = "dec2flt",
reason = "internal routines only exposed for testing")]
use prelude::v1::*;
use num::ParseFloatError as PFE;
use num::FloatErrorKind;
use self::parse::{parse_decimal, Decimal, Sign};
use self::parse::ParseResult::{self, Valid, ShortcutToInf, ShortcutToZero};
use self::num::digits_to_big;
use self::rawfp::RawFloat;
mod algorithm;
mod table;
mod num;
// These two have their own tests.
pub mod rawfp;
pub mod parse;
/// Entry point for decimal-to-f32 conversion.
pub fn to_f32(s: &str) -> Result<f32, PFE> {
dec2flt(s)
}
/// Entry point for decimal-to-f64 conversion.
pub fn to_f64(s: &str) -> Result<f64, PFE> {
dec2flt(s)
}
/// Split decimal string into sign and the rest, without inspecting or validating the rest.
fn extract_sign(s: &str) -> (Sign, &str) {
match s.as_bytes()[0] {
b'+' => (Sign::Positive, &s[1..]),
b'-' => (Sign::Negative, &s[1..]),
// If the string is invalid, we never use the sign, so we don't need to validate here.
_ => (Sign::Positive, s),
}
}
/// Convert a decimal string into a floating point number.
fn dec2flt<T: RawFloat>(s: &str) -> Result<T, PFE> {
if s.is_empty() {
return Err(PFE { __kind: FloatErrorKind::Empty });
}
let (sign, s) = extract_sign(s);
let flt = match parse_decimal(s) {
Valid(decimal) => try!(convert(decimal)),
ShortcutToInf => T::infinity(),
ShortcutToZero => T::zero(),
ParseResult::Invalid => match s {
"inf" => T::infinity(),
"NaN" => T::nan(),
_ => { return Err(PFE { __kind: FloatErrorKind::Invalid }); }
}
};
match sign {
Sign::Positive => Ok(flt),
Sign::Negative => Ok(-flt),
}
}
/// The main workhorse for the decimal-to-float conversion: Orchestrate all the preprocessing
/// and figure out which algorithm should do the actual conversion.
fn convert<T: RawFloat>(mut decimal: Decimal) -> Result<T, PFE> {
simplify(&mut decimal);
if let Some(x) = trivial_cases(&decimal) {
return Ok(x);
}
// AlgorithmM and AlgorithmR both compute approximately `f * 10^e`.
let max_digits = decimal.integral.len() + decimal.fractional.len() +
decimal.exp.abs() as usize;
// Remove/shift out the decimal point.
let e = decimal.exp - decimal.fractional.len() as i64;
if let Some(x) = algorithm::fast_path(decimal.integral, decimal.fractional, e) {
return Ok(x);
}
// Big32x40 is limited to 1280 bits, which translates to about 385 decimal digits.
// If we exceed this, perhaps while calculating `f * 10^e` in Algorithm R or Algorithm M,
// we'll crash. So we error out before getting too close, with a generous safety margin.
if max_digits > 375 {
return Err(PFE { __kind: FloatErrorKind::Invalid });
}
let f = digits_to_big(decimal.integral, decimal.fractional);
// Now the exponent certainly fits in 16 bit, which is used throughout the main algorithms.
let e = e as i16;
// FIXME These bounds are rather conservative. A more careful analysis of the failure modes
// of Bellerophon could allow using it in more cases for a massive speed up.
let exponent_in_range = table::MIN_E <= e && e <= table::MAX_E;
let value_in_range = max_digits <= T::max_normal_digits();
if exponent_in_range && value_in_range {
Ok(algorithm::bellerophon(&f, e))
} else {
Ok(algorithm::algorithm_m(&f, e))
}
}
// As written, this optimizes badly (see #27130, though it refers to an old version of the code).
// `inline(always)` is a workaround for that. There are only two call sites overall and it doesn't
// make code size worse.
/// Strip zeros where possible, even when this requires changing the exponent
#[inline(always)]
fn simplify(decimal: &mut Decimal) {
let is_zero = &|&&d: &&u8| -> bool { d == b'0' };
// Trimming these zeros does not change anything but may enable the fast path (< 15 digits).
let leading_zeros = decimal.integral.iter().take_while(is_zero).count();
decimal.integral = &decimal.integral[leading_zeros..];
let trailing_zeros = decimal.fractional.iter().rev().take_while(is_zero).count();
let end = decimal.fractional.len() - trailing_zeros;
decimal.fractional = &decimal.fractional[..end];
// Simplify numbers of the form 0.0...x and x...0.0, adjusting the exponent accordingly.
// This may not always be a win (possibly pushes some numbers out of the fast path), but it
// simplifies other parts significantly (notably, approximating the magnitude of the value).
if decimal.integral.is_empty() {
let leading_zeros = decimal.fractional.iter().take_while(is_zero).count();
decimal.fractional = &decimal.fractional[leading_zeros..];
decimal.exp -= leading_zeros as i64;
} else if decimal.fractional.is_empty() {
let trailing_zeros = decimal.integral.iter().rev().take_while(is_zero).count();
let end = decimal.integral.len() - trailing_zeros;
decimal.integral = &decimal.integral[..end];
decimal.exp += trailing_zeros as i64;
}
}
/// Detect obvious overflows and underflows without even looking at the decimal digits.
fn trivial_cases<T: RawFloat>(decimal: &Decimal) -> Option<T> {
// There were zeros but they were stripped by simplify()
if decimal.integral.is_empty() && decimal.fractional.is_empty() {
return Some(T::zero());
}
// This is a crude approximation of ceil(log10(the real value)).
let max_place = decimal.exp + decimal.integral.len() as i64;
if max_place > T::inf_cutoff() {
return Some(T::infinity());
} else if max_place < T::zero_cutoff() {
return Some(T::zero());
}
None
}

95
src/libcore/num/dec2flt/num.rs Normal file
View file

@ -0,0 +1,95 @@
// Copyright 2015 The Rust Project Developers. See the COPYRIGHT
// file at the top-level directory of this distribution and at
// http://rust-lang.org/COPYRIGHT.
//
// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
// http://www.apache.org/licenses/LICENSE-2.0> or the MIT license
// <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your
// option. This file may not be copied, modified, or distributed
// except according to those terms.
//! Utility functions for bignums that don't make too much sense to turn into methods.
// FIXME This module's name is a bit unfortunate, since other modules also import `core::num`.
use prelude::v1::*;
use cmp::Ordering::{self, Less, Equal, Greater};
use num::flt2dec::bignum::Big32x40;
pub type Big = Big32x40;
/// Test whether truncating all bits less significant than `ones_place` introduces
/// a relative error less, equal, or greater than 0.5 ULP.
pub fn compare_with_half_ulp(f: &Big, ones_place: usize) -> Ordering {
if ones_place == 0 {
return Less;
}
let half_bit = ones_place - 1;
if f.get_bit(half_bit) == 0 {
// < 0.5 ULP
return Less;
}
// If all remaining bits are zero, it's = 0.5 ULP, otherwise > 0.5
// If there are no more bits (half_bit == 0), the below also correctly returns Equal.
for i in 0..half_bit {
if f.get_bit(i) == 1 {
return Greater;
}
}
Equal
}
/// Convert an ASCII string containing only decimal digits to a `u64`.
///
/// Does not perform checks for overflow or invalid characters, so if the caller is not careful,
/// the result is bogus and can panic (though it won't be `unsafe`). Additionally, empty strings
/// are treated as zero. This function exists because
///
/// 1. using `FromStr` on `&[u8]` requires `from_utf8_unchecked`, which is bad, and
/// 2. piecing together the results of `integral.parse()` and `fractional.parse()` is
/// more complicated than this entire function.
pub fn from_str_unchecked<'a, T>(bytes: T) -> u64 where T : IntoIterator<Item=&'a u8> {
let mut result = 0;
for &c in bytes {
result = result * 10 + (c - b'0') as u64;
}
result
}
/// Convert a string of ASCII digits into a bignum.
///
/// Like `from_str_unchecked`, this function relies on the parser to weed out non-digits.
pub fn digits_to_big(integral: &[u8], fractional: &[u8]) -> Big {
let mut f = Big::from_small(0);
for &c in integral.iter().chain(fractional) {
let n = (c - b'0') as u32;
f.mul_small(10);
f.add_small(n);
}
f
}
/// Unwraps a bignum into a 64 bit integer. Panics if the number is too large.
pub fn to_u64(x: &Big) -> u64 {
assert!(x.bit_length() < 64);
let d = x.digits();
if d.len() < 2 {
d[0] as u64
} else {
(d[1] as u64) << 32 | d[0] as u64
}
}
/// Extract a range of bits.
/// Index 0 is the least significant bit and the range is half-open as usual.
/// Panics if asked to extract more bits than fit into the return type.
pub fn get_bits(x: &Big, start: usize, end: usize) -> u64 {
assert!(end - start <= 64);
let mut result: u64 = 0;
for i in (start..end).rev() {
result = result << 1 | x.get_bit(i) as u64;
}
result
}

128
src/libcore/num/dec2flt/parse.rs Normal file
View file

@ -0,0 +1,128 @@
// Copyright 2015 The Rust Project Developers. See the COPYRIGHT
// file at the top-level directory of this distribution and at
// http://rust-lang.org/COPYRIGHT.
//
// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
// http://www.apache.org/licenses/LICENSE-2.0> or the MIT license
// <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your
// option. This file may not be copied, modified, or distributed
// except according to those terms.
//! Validating and decomposing a decimal string of the form:
//!
//! `(digits | digits? '.'? digits?) (('e' | 'E') ('+' | '-')? digits)?`
//!
//! In other words, standard floating-point syntax, with two exceptions: No sign, and no
//! handling of "inf" and "NaN". These are handled by the driver function (super::dec2flt).
//!
//! Although recognizing valid inputs is relatively easy, this module also has to reject the
//! countless invalid variations, never panic, and perform numerous checks that the other
//! modules rely on to not panic (or overflow) in turn.
//! To make matters worse, all that happens in a single pass over the input.
//! So, be careful when modifying anything, and double-check with the other modules.
use prelude::v1::*;
use super::num;
use self::ParseResult::{Valid, ShortcutToInf, ShortcutToZero, Invalid};
#[derive(Debug)]
pub enum Sign {
Positive,
Negative,
}
#[derive(Debug, PartialEq, Eq)]
/// The interesting parts of a decimal string.
pub struct Decimal<'a> {
pub integral: &'a [u8],
pub fractional: &'a [u8],
/// The decimal exponent, guaranteed to have fewer than 18 decimal digits.
pub exp: i64,
}
impl<'a> Decimal<'a> {
pub fn new(integral: &'a [u8], fractional: &'a [u8], exp: i64) -> Decimal<'a> {
Decimal { integral: integral, fractional: fractional, exp: exp }
}
}
#[derive(Debug, PartialEq, Eq)]
pub enum ParseResult<'a> {
Valid(Decimal<'a>),
ShortcutToInf,
ShortcutToZero,
Invalid,
}
/// Check if the input string is a valid floating point number and if so, locate the integral
/// part, the fractional part, and the exponent in it. Does not handle signs.
pub fn parse_decimal(s: &str) -> ParseResult {
let s = s.as_bytes();
let (integral, s) = eat_digits(s);
match s.first() {
None => Valid(Decimal::new(integral, b"", 0)),
Some(&b'e') | Some(&b'E') => {
if integral.is_empty() {
return Invalid; // No digits before 'e'
}
parse_exp(integral, b"", &s[1..])
}
Some(&b'.') => {
let (fractional, s) = eat_digits(&s[1..]);
if integral.is_empty() && fractional.is_empty() && s.is_empty() {
// For historic reasons "." is a valid input.
return Valid(Decimal::new(b"", b"", 0));
}
match s.first() {
None => Valid(Decimal::new(integral, fractional, 0)),
Some(&b'e') | Some(&b'E') => parse_exp(integral, fractional, &s[1..]),
_ => Invalid, // Trailing junk after fractional part
}
}
_ => Invalid, // Trailing junk after first digit string
}
}
/// Carve off decimal digits up to the first non-digit character.
fn eat_digits(s: &[u8]) -> (&[u8], &[u8]) {
let mut i = 0;
while i < s.len() && b'0' <= s[i] && s[i] <= b'9' {
i += 1;
}
(&s[..i], &s[i..])
}
/// Exponent extraction and error checking.
fn parse_exp<'a>(integral: &'a [u8], fractional: &'a [u8], rest: &'a [u8]) -> ParseResult<'a> {
let (sign, rest) = match rest.first() {
Some(&b'-') => (Sign::Negative, &rest[1..]),
Some(&b'+') => (Sign::Positive, &rest[1..]),
_ => (Sign::Positive, rest),
};
let (mut number, trailing) = eat_digits(rest);
if !trailing.is_empty() {
return Invalid; // Trailing junk after exponent
}
if number.is_empty() {
return Invalid; // Empty exponent
}
// At this point, we certainly have a valid string of digits. It may be too long to put into
// an `i64`, but if it's that huge, the input is certainly zero or infinity. Since each zero
// in the decimal digits only adjusts the exponent by +/- 1, at exp = 10^18 the input would
// have to be 17 exabyte (!) of zeros to get even remotely close to being finite.
// This is not exactly a use case we need to cater to.
while number.first() == Some(&b'0') {
number = &number[1..];
}
if number.len() >= 18 {
return match sign {
Sign::Positive => ShortcutToInf,
Sign::Negative => ShortcutToZero,
};
}
let abs_exp = num::from_str_unchecked(number);
let e = match sign {
Sign::Positive => abs_exp as i64,
Sign::Negative => -(abs_exp as i64),
};
Valid(Decimal::new(integral, fractional, e))
}

356
src/libcore/num/dec2flt/rawfp.rs Normal file
View file

@ -0,0 +1,356 @@
// Copyright 2015 The Rust Project Developers. See the COPYRIGHT
// file at the top-level directory of this distribution and at
// http://rust-lang.org/COPYRIGHT.
//
// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
// http://www.apache.org/licenses/LICENSE-2.0> or the MIT license
// <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your
// option. This file may not be copied, modified, or distributed
// except according to those terms.
//! Bit fiddling on positive IEEE 754 floats. Negative numbers aren't and needn't be handled.
//! Normal floating point numbers have a canonical representation as (frac, exp) such that the
//! value is 2^exp * (1 + sum(frac[N-i] / 2^i)) where N is the number of bits. Subnormals are
//! slightly different and weird, but the same principle applies.
//!
//! Here, however, we represent them as (sig, k) with f positive, such that the value is f * 2^e.
//! Besides making the "hidden bit" explicit, this changes the exponent by the so-called
//! mantissa shift.
//!
//! Put another way, normally floats are written as (1) but here they are written as (2):
//!
//! 1. `1.101100...11 * 2^m`
//! 2. `1101100...11 * 2^n`
//!
//! We call (1) the **fractional representation** and (2) the **integral representation**.
//!
//! Many functions in this module only handle normal numbers. The dec2flt routines conservatively
//! take the universally-correct slow path (Algorithm M) for very small and very large numbers.
//! That algorithm needs only next_float() which does handle subnormals and zeros.
use prelude::v1::*;
use u32;
use cmp::Ordering::{Less, Equal, Greater};
use ops::{Mul, Div, Neg};
use fmt::{Debug, LowerExp};
use mem::transmute;
use num::flt2dec::strategy::grisu::Fp;
use num::FpCategory::{Infinite, Zero, Subnormal, Normal, Nan};
use num::Float;
use super::num::{self, Big};
#[derive(Copy, Clone, Debug)]
pub struct Unpacked {
pub sig: u64,
pub k: i16,
}
impl Unpacked {
pub fn new(sig: u64, k: i16) -> Self {
Unpacked { sig: sig, k: k }
}
}
/// A helper trait to avoid duplicating basically all the conversion code for `f32` and `f64`.
///
/// 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.
/// Inherits from `Float` because there is some overlap, but all the reused methods are trivial.
/// The "methods" (pseudo-constants) with default implementation should not be overriden.
pub trait RawFloat : Float + Copy + Debug + LowerExp
+ Mul<Output=Self> + Div<Output=Self> + Neg<Output=Self>
{
/// Get the raw binary representation of the float.
fn transmute(self) -> u64;
/// Transmute the raw binary representation into a float.
fn from_bits(bits: u64) -> Self;
/// Decode the float.
fn unpack(self) -> Unpacked;
/// Cast from a small integer that can be represented exactly. Panic if the integer can't be
/// represented, the other code in this module makes sure to never let that happen.
fn from_int(x: u64) -> Self;
// FIXME Everything that follows should be associated constants, but taking the value of an
// associated constant from a type parameter does not work (yet?)
// A possible workaround is having a `FloatInfo` struct for all the constants, but so far
// the methods aren't painful enough to rewrite.
/// What the name says. It's easier to hard code than juggling intrinsics and
/// hoping LLVM constant folds it.
fn ceil_log5_of_max_sig() -> i16;
// A conservative bound on the decimal digits of inputs that can't produce overflow or zero or
/// subnormals. Probably the decimal exponent of the maximum normal value, hence the name.
fn max_normal_digits() -> usize;
/// When the most significant decimal digit has a place value greater than this, the number
/// is certainly rounded to infinity.
fn inf_cutoff() -> i64;
/// When the most significant decimal digit has a place value less than this, the number
/// is certainly rounded to zero.
fn zero_cutoff() -> i64;
/// The number of bits in the exponent.
fn exp_bits() -> u8;
/// The number of bits in the singificand, *including* the hidden bit.
fn sig_bits() -> u8;
/// The number of bits in the singificand, *excluding* the hidden bit.
fn explicit_sig_bits() -> u8 {
Self::sig_bits() - 1
}
/// The maximum legal exponent in fractional representation.
fn max_exp() -> i16 {
(1 << (Self::exp_bits() - 1)) - 1
}
/// The minimum legal exponent in fractional representation, excluding subnormals.
fn min_exp() -> i16 {
-Self::max_exp() + 1
}
/// `MAX_EXP` for integral representation, i.e., with the shift applied.
fn max_exp_int() -> i16 {
Self::max_exp() - (Self::sig_bits() as i16 - 1)
}
/// `MAX_EXP` encoded (i.e., with offset bias)
fn max_encoded_exp() -> i16 {
(1 << Self::exp_bits()) - 1
}
/// `MIN_EXP` for integral representation, i.e., with the shift applied.
fn min_exp_int() -> i16 {
Self::min_exp() - (Self::sig_bits() as i16 - 1)
}
/// The maximum normalized singificand in integral representation.
fn max_sig() -> u64 {
(1 << Self::sig_bits()) - 1
}
/// The minimal normalized significand in integral representation.
fn min_sig() -> u64 {
1 << (Self::sig_bits() - 1)
}
}
impl RawFloat for f32 {
fn sig_bits() -> u8 {
24
}
fn exp_bits() -> u8 {
8
}
fn ceil_log5_of_max_sig() -> i16 {
11
}
fn transmute(self) -> u64 {
let bits: u32 = unsafe { transmute(self) };
bits as u64
}
fn from_bits(bits: u64) -> f32 {
assert!(bits < u32::MAX as u64, "f32::from_bits: too many bits");
unsafe { transmute(bits as u32) }
}
fn unpack(self) -> Unpacked {
let (sig, exp, _sig) = self.integer_decode();
Unpacked::new(sig, exp)
}
fn from_int(x: u64) -> f32 {
// rkruppe is uncertain whether `as` rounds correctly on all platforms.
debug_assert!(x as f32 == fp_to_float(Fp { f: x, e: 0 }));
x as f32
}
fn max_normal_digits() -> usize {
35
}
fn inf_cutoff() -> i64 {
40
}
fn zero_cutoff() -> i64 {
-48
}
}
impl RawFloat for f64 {
fn sig_bits() -> u8 {
53
}
fn exp_bits() -> u8 {
11
}
fn ceil_log5_of_max_sig() -> i16 {
23
}
fn transmute(self) -> u64 {
let bits: u64 = unsafe { transmute(self) };
bits
}
fn from_bits(bits: u64) -> f64 {
unsafe { transmute(bits) }
}
fn unpack(self) -> Unpacked {
let (sig, exp, _sig) = self.integer_decode();
Unpacked::new(sig, exp)
}
fn from_int(x: u64) -> f64 {
// rkruppe is uncertain whether `as` rounds correctly on all platforms.
debug_assert!(x as f64 == fp_to_float(Fp { f: x, e: 0 }));
x as f64
}
fn max_normal_digits() -> usize {
305
}
fn inf_cutoff() -> i64 {
310
}
fn zero_cutoff() -> i64 {
-326
}
}
/// Convert an Fp to the closest f64. Only handles number that fit into a normalized f64.
pub fn fp_to_float<T: RawFloat>(x: Fp) -> T {
let x = x.normalize();
// x.f is 64 bit, so x.e has a mantissa shift of 63
let e = x.e + 63;
if e > T::max_exp() {
panic!("fp_to_float: exponent {} too large", e)
} else if e > T::min_exp() {
encode_normal(round_normal::<T>(x))
} else {
panic!("fp_to_float: exponent {} too small", e)
}
}
/// Round the 64-bit significand to 53 bit with half-to-even. Does not handle exponent overflow.
pub fn round_normal<T: RawFloat>(x: Fp) -> Unpacked {
let excess = 64 - T::sig_bits() as i16;
let half: u64 = 1 << (excess - 1);
let (q, rem) = (x.f >> excess, x.f & ((1 << excess) - 1));
assert_eq!(q << excess | rem, x.f);
// Adjust mantissa shift
let k = x.e + excess;
if rem < half {
Unpacked::new(q, k)
} else if rem == half && (q % 2) == 0 {
Unpacked::new(q, k)
} else if q == T::max_sig() {
Unpacked::new(T::min_sig(), k + 1)
} else {
Unpacked::new(q + 1, k)
}
}
/// Inverse of `RawFloat::unpack()` for normalized numbers.
/// Panics if the significand or exponent are not valid for normalized numbers.
pub fn encode_normal<T: RawFloat>(x: Unpacked) -> T {
debug_assert!(T::min_sig() <= x.sig && x.sig <= T::max_sig(),
"encode_normal: significand not normalized");
// Remove the hidden bit
let sig_enc = x.sig & !(1 << T::explicit_sig_bits());
// Adjust the exponent for exponent bias and mantissa shift
let k_enc = x.k + T::max_exp() + T::explicit_sig_bits() as i16;
debug_assert!(k_enc != 0 && k_enc < T::max_encoded_exp(),
"encode_normal: exponent out of range");
// Leave sign bit at 0 ("+"), our numbers are all positive
let bits = (k_enc as u64) << T::explicit_sig_bits() | sig_enc;
T::from_bits(bits)
}
/// Construct the subnormal. A mantissa of 0 is allowed and constructs zero.
pub fn encode_subnormal<T: RawFloat>(significand: u64) -> T {
assert!(significand < T::min_sig(), "encode_subnormal: not actually subnormal");
// Êncoded exponent is 0, the sign bit is 0, so we just have to reinterpret the bits.
T::from_bits(significand)
}
/// Approximate a bignum with an Fp. Rounds within 0.5 ULP with half-to-even.
pub fn big_to_fp(f: &Big) -> Fp {
let end = f.bit_length();
assert!(end != 0, "big_to_fp: unexpectedly, input is zero");
let start = end.saturating_sub(64);
let leading = num::get_bits(f, start, end);
// We cut off all bits prior to the index `start`, i.e., we effectively right-shift by
// an amount of `start`, so this is also the exponent we need.
let e = start as i16;
let rounded_down = Fp { f: leading, e: e }.normalize();
// Round (half-to-even) depending on the truncated bits.
match num::compare_with_half_ulp(f, start) {
Less => rounded_down,
Equal if leading % 2 == 0 => rounded_down,
Equal | Greater => match leading.checked_add(1) {
Some(f) => Fp { f: f, e: e }.normalize(),
None => Fp { f: 1 << 63, e: e + 1 },
}
}
}
/// Find the largest floating point number strictly smaller than the argument.
/// Does not handle subnormals, zero, or exponent underflow.
pub fn prev_float<T: RawFloat>(x: T) -> T {
match x.classify() {
Infinite => panic!("prev_float: argument is infinite"),
Nan => panic!("prev_float: argument is NaN"),
Subnormal => panic!("prev_float: argument is subnormal"),
Zero => panic!("prev_float: argument is zero"),
Normal => {
let Unpacked { sig, k } = x.unpack();
if sig == T::min_sig() {
encode_normal(Unpacked::new(T::max_sig(), k - 1))
} else {
encode_normal(Unpacked::new(sig - 1, k))
}
}
}
}
// Find the smallest floating point number strictly larger than the argument.
// This operation is saturating, i.e. next_float(inf) == inf.
// Unlike most code in this module, this function does handle zero, subnormals, and infinities.
// However, like all other code here, it does not deal with NaN and negative numbers.
pub fn next_float<T: RawFloat>(x: T) -> T {
match x.classify() {
Nan => panic!("next_float: argument is NaN"),
Infinite => T::infinity(),
// This seems too good to be true, but it works.
// 0.0 is encoded as the all-zero word. Subnormals are 0x000m...m where m is the mantissa.
// In particular, the smallest subnormal is 0x0...01 and the largest is 0x000F...F.
// The smallest normal number is 0x0010...0, so this corner case works as well.
// If the increment overflows the mantissa, the carry bit increments the exponent as we
// want, and the mantissa bits become zero. Because of the hidden bit convention, this
// too is exactly what we want!
// Finally, f64::MAX + 1 = 7eff...f + 1 = 7ff0...0 = f64::INFINITY.
Zero | Subnormal | Normal => {
let bits: u64 = x.transmute();
T::from_bits(bits + 1)
}
}
}

1239
src/libcore/num/dec2flt/table.rs Normal file
View file

File diff suppressed because it is too large Load diff

21
src/libcore/num/flt2dec/bignum.rs
View file

@ -171,20 +171,19 @@ macro_rules! define_bignum {
pub fn bit_length(&self) -> usize {
use mem;
let digitbits = mem::size_of::<$ty>()* 8;
// Skip over the most significant digits which are zero.
let nonzero = match self.digits().iter().rposition(|&x| x != 0) {
Some(n) => {
let end = self.size - n;
&self.digits()[..end]
}
None => {
// There are no non-zero digits, i.e. the number is zero.
return 0;
}
};
let digits = self.digits();
let zeros = digits.iter().rev().take_while(|&&x| x == 0).count();
let end = digits.len() - zeros;
let nonzero = &digits[..end];
if nonzero.is_empty() {
// There are no non-zero digits, i.e. the number is zero.
return 0;
}
// This could be optimized with leading_zeros() and bit shifts, but that's
// probably not worth the hassle.
let digitbits = mem::size_of::<$ty>()* 8;
let mut i = nonzero.len() * digitbits - 1;
while self.get_bit(i) == 0 {
i -= 1;

9
src/libcore/num/mod.rs
View file

@ -44,6 +44,7 @@ pub struct Wrapping<T>(#[stable(feature = "rust1", since = "1.0.0")] pub T);
pub mod wrapping;
pub mod flt2dec;
pub mod dec2flt;
/// Types that have a "zero" value.
///
@ -1364,7 +1365,7 @@ pub trait Float {
}
macro_rules! from_str_float_impl {
($t:ty) => {
($t:ty, $func:ident) => {
#[stable(feature = "rust1", since = "1.0.0")]
impl FromStr for $t {
type Err = ParseFloatError;
@ -1397,13 +1398,13 @@ macro_rules! from_str_float_impl {
#[inline]
#[allow(deprecated)]
fn from_str(src: &str) -> Result<Self, ParseFloatError> {
Self::from_str_radix(src, 10)
dec2flt::$func(src)
}
}
}
}
from_str_float_impl!(f32);
from_str_float_impl!(f64);
from_str_float_impl!(f32, to_f32);
from_str_float_impl!(f64, to_f64);
macro_rules! from_str_radix_int_impl {
($($t:ty)*) => {$(

1
src/libcoretest/lib.rs
View file

@ -19,6 +19,7 @@
#![feature(float_extras)]
#![feature(float_from_str_radix)]
#![feature(flt2dec)]
#![feature(dec2flt)]
#![feature(fmt_radix)]
#![feature(hash_default)]
#![feature(hasher_write)]

174
src/libcoretest/num/dec2flt/mod.rs Normal file
View file

@ -0,0 +1,174 @@
// Copyright 2015 The Rust Project Developers. See the COPYRIGHT
// file at the top-level directory of this distribution and at
// http://rust-lang.org/COPYRIGHT.
//
// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
// http://www.apache.org/licenses/LICENSE-2.0> or the MIT license
// <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your
// option. This file may not be copied, modified, or distributed
// except according to those terms.
#![allow(overflowing_literals)]
use std::{i64, f32, f64};
use test;
use core::num::dec2flt::{to_f32, to_f64};
mod parse;
mod rawfp;
// Take an float literal, turn it into a string in various ways (that are all trusted
// to be correct) and see if those strings are parsed back to the value of the literal.
// Requires a *polymorphic literal*, i.e. one that can serve as f64 as well as f32.
macro_rules! test_literal {
($x: expr) => ({
let x32: f32 = $x;
let x64: f64 = $x;
let inputs = &[stringify!($x).into(), format!("{:?}", x64), format!("{:e}", x64)];
for input in inputs {
if input != "inf" {
assert_eq!(to_f64(input), Ok(x64));
assert_eq!(to_f32(input), Ok(x32));
let neg_input = &format!("-{}", input);
assert_eq!(to_f64(neg_input), Ok(-x64));
assert_eq!(to_f32(neg_input), Ok(-x32));
}
}
})
}
#[test]
fn ordinary() {
test_literal!(1.0);
test_literal!(3e-5);
test_literal!(0.1);
test_literal!(12345.);
test_literal!(0.9999999);
test_literal!(2.2250738585072014e-308);
}
#[test]
fn special_code_paths() {
test_literal!(36893488147419103229.0); // 2^65 - 3, triggers half-to-even with even significand
test_literal!(101e-33); // Triggers the tricky underflow case in AlgorithmM (for f32)
test_literal!(1e23); // Triggers AlgorithmR
test_literal!(2075e23); // Triggers another path through AlgorithmR
test_literal!(8713e-23); // ... and yet another.
}
#[test]
fn large() {
test_literal!(1e300);
test_literal!(123456789.34567e250);
test_literal!(943794359898089732078308743689303290943794359843568973207830874368930329.);
}
#[test]
fn subnormals() {
test_literal!(5e-324);
test_literal!(91e-324);
test_literal!(1e-322);
test_literal!(13245643e-320);
test_literal!(2.22507385851e-308);
test_literal!(2.1e-308);
test_literal!(4.9406564584124654e-324);
}
#[test]
fn infinity() {
test_literal!(1e400);
test_literal!(1e309);
test_literal!(2e308);
test_literal!(1.7976931348624e308);
}
#[test]
fn zero() {
test_literal!(0.0);
test_literal!(1e-325);
test_literal!(1e-326);
test_literal!(1e-500);
}
#[test]
fn lonely_dot() {
assert_eq!(to_f64("."), Ok(0.0));
}
#[test]
fn nan() {
assert!(to_f64("NaN").unwrap().is_nan());
assert!(to_f32("NaN").unwrap().is_nan());
}
#[test]
fn inf() {
assert_eq!(to_f64("inf"), Ok(f64::INFINITY));
assert_eq!(to_f64("-inf"), Ok(f64::NEG_INFINITY));
assert_eq!(to_f32("inf"), Ok(f32::INFINITY));
assert_eq!(to_f32("-inf"), Ok(f32::NEG_INFINITY));
}
#[test]
fn massive_exponent() {
let max = i64::MAX;
assert_eq!(to_f64(&format!("1e{}000", max)), Ok(f64::INFINITY));
assert_eq!(to_f64(&format!("1e-{}000", max)), Ok(0.0));
assert_eq!(to_f64(&format!("1e{}000", max)), Ok(f64::INFINITY));
}
#[bench]
fn bench_0(b: &mut test::Bencher) {
b.iter(|| to_f64("0.0"));
}
#[bench]
fn bench_42(b: &mut test::Bencher) {
b.iter(|| to_f64("42"));
}
#[bench]
fn bench_huge_int(b: &mut test::Bencher) {
// 2^128 - 1
b.iter(|| to_f64("170141183460469231731687303715884105727"));
}
#[bench]
fn bench_short_decimal(b: &mut test::Bencher) {
b.iter(|| to_f64("1234.5678"));
}
#[bench]
fn bench_pi_long(b: &mut test::Bencher) {
b.iter(|| to_f64("3.14159265358979323846264338327950288"));
}
#[bench]
fn bench_pi_short(b: &mut test::Bencher) {
b.iter(|| to_f64("3.141592653589793"))
}
#[bench]
fn bench_1e150(b: &mut test::Bencher) {
b.iter(|| to_f64("1e150"));
}
#[bench]
fn bench_long_decimal_and_exp(b: &mut test::Bencher) {
b.iter(|| to_f64("727501488517303786137132964064381141071e-123"));
}
#[bench]
fn bench_min_subnormal(b: &mut test::Bencher) {
b.iter(|| to_f64("5e-324"));
}
#[bench]
fn bench_min_normal(b: &mut test::Bencher) {
b.iter(|| to_f64("2.2250738585072014e-308"));
}
#[bench]
fn bench_max(b: &mut test::Bencher) {
b.iter(|| to_f64("1.7976931348623157e308"));
}

52
src/libcoretest/num/dec2flt/parse.rs Normal file
View file

@ -0,0 +1,52 @@
// Copyright 2015 The Rust Project Developers. See the COPYRIGHT
// file at the top-level directory of this distribution and at
// http://rust-lang.org/COPYRIGHT.
//
// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
// http://www.apache.org/licenses/LICENSE-2.0> or the MIT license
// <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your
// option. This file may not be copied, modified, or distributed
// except according to those terms.
use std::iter;
use core::num::dec2flt::parse::{Decimal, parse_decimal};
use core::num::dec2flt::parse::ParseResult::{Valid, Invalid};
#[test]
fn missing_pieces() {
let permutations = &[".e", "1e", "e4", "e", ".12e", "321.e", "32.12e+", "12.32e-"];
for &s in permutations {
assert_eq!(parse_decimal(s), Invalid);
}
}
#[test]
fn invalid_chars() {
let invalid = "r,?<j";
let valid_strings = &["123", "666.", ".1", "5e1", "7e-3", "0.0e+1"];
for c in invalid.chars() {
for s in valid_strings {
for i in 0..s.len() {
let mut input = String::new();
input.push_str(s);
input.insert(i, c);
assert!(parse_decimal(&input) == Invalid, "did not reject invalid {:?}", input);
}
}
}
}
#[test]
fn valid() {
assert_eq!(parse_decimal("123.456e789"), Valid(Decimal::new(b"123", b"456", 789)));
assert_eq!(parse_decimal("123.456e+789"), Valid(Decimal::new(b"123", b"456", 789)));
assert_eq!(parse_decimal("123.456e-789"), Valid(Decimal::new(b"123", b"456", -789)));
assert_eq!(parse_decimal(".050"), Valid(Decimal::new(b"", b"050", 0)));
assert_eq!(parse_decimal("999"), Valid(Decimal::new(b"999", b"", 0)));
assert_eq!(parse_decimal("1.e300"), Valid(Decimal::new(b"1", b"", 300)));
assert_eq!(parse_decimal(".1e300"), Valid(Decimal::new(b"", b"1", 300)));
assert_eq!(parse_decimal("101e-33"), Valid(Decimal::new(b"101", b"", -33)));
let zeros: String = iter::repeat('0').take(25).collect();
let s = format!("1.5e{}", zeros);
assert_eq!(parse_decimal(&s), Valid(Decimal::new(b"1", b"5", 0)));
}

139
src/libcoretest/num/dec2flt/rawfp.rs Normal file
View file

@ -0,0 +1,139 @@
// Copyright 2015 The Rust Project Developers. See the COPYRIGHT
// file at the top-level directory of this distribution and at
// http://rust-lang.org/COPYRIGHT.
//
// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
// http://www.apache.org/licenses/LICENSE-2.0> or the MIT license
// <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your
// option. This file may not be copied, modified, or distributed
// except according to those terms.
use std::f64;
use core::num::flt2dec::strategy::grisu::Fp;
use core::num::dec2flt::rawfp::{fp_to_float, prev_float, next_float, round_normal};
#[test]
fn fp_to_float_half_to_even() {
fn is_normalized(sig: u64) -> bool {
// intentionally written without {min,max}_sig() as a sanity check
sig >> 52 == 1 && sig >> 53 == 0
}
fn conv(sig: u64) -> u64 {
// The significands are perfectly in range, so the exponent should not matter
let (m1, e1, _) = fp_to_float::<f64>(Fp { f: sig, e: 0 }).integer_decode();
assert_eq!(e1, 0 + 64 - 53);
let (m2, e2, _) = fp_to_float::<f64>(Fp { f: sig, e: 55 }).integer_decode();
assert_eq!(e2, 55 + 64 - 53);
assert_eq!(m2, m1);
let (m3, e3, _) = fp_to_float::<f64>(Fp { f: sig, e: -78 }).integer_decode();
assert_eq!(e3, -78 + 64 - 53);
assert_eq!(m3, m2);
m3
}
let odd = 0x1F_EDCB_A012_345F;
let even = odd - 1;
assert!(is_normalized(odd));
assert!(is_normalized(even));
assert_eq!(conv(odd << 11), odd);
assert_eq!(conv(even << 11), even);
assert_eq!(conv(odd << 11 | 1 << 10), odd + 1);
assert_eq!(conv(even << 11 | 1 << 10), even);
assert_eq!(conv(even << 11 | 1 << 10 | 1), even + 1);
assert_eq!(conv(odd << 11 | 1 << 9), odd);
assert_eq!(conv(even << 11 | 1 << 9), even);
assert_eq!(conv(odd << 11 | 0x7FF), odd + 1);
assert_eq!(conv(even << 11 | 0x7FF), even + 1);
assert_eq!(conv(odd << 11 | 0x3FF), odd);
assert_eq!(conv(even << 11 | 0x3FF), even);
}
#[test]
fn integers_to_f64() {
assert_eq!(fp_to_float::<f64>(Fp { f: 1, e: 0 }), 1.0);
assert_eq!(fp_to_float::<f64>(Fp { f: 42, e: 7 }), (42 << 7) as f64);
assert_eq!(fp_to_float::<f64>(Fp { f: 1 << 20, e: 30 }), (1u64 << 50) as f64);
assert_eq!(fp_to_float::<f64>(Fp { f: 4, e: -3 }), 0.5);
}
const SOME_FLOATS: [f64; 9] =
[0.1f64, 33.568, 42.1e-5, 777.0e9, 1.1111, 0.347997,
9843579834.35892, 12456.0e-150, 54389573.0e-150];
#[test]
fn human_f64_roundtrip() {
for &x in &SOME_FLOATS {
let (f, e, _) = x.integer_decode();
let fp = Fp { f: f, e: e};
assert_eq!(fp_to_float::<f64>(fp), x);
}
}
#[test]
fn rounding_overflow() {
let x = Fp { f: 0xFF_FF_FF_FF_FF_FF_FF_00u64, e: 42 };
let rounded = round_normal::<f64>(x);
let adjusted_k = x.e + 64 - 53;
assert_eq!(rounded.sig, 1 << 52);
assert_eq!(rounded.k, adjusted_k + 1);
}
#[test]
fn prev_float_monotonic() {
let mut x = 1.0;
for _ in 0..100 {
let x1 = prev_float(x);
assert!(x1 < x);
assert!(x - x1 < 1e-15);
x = x1;
}
}
const MIN_SUBNORMAL: f64 = 5e-324;
#[test]
fn next_float_zero() {
let tiny = next_float(0.0);
assert_eq!(tiny, MIN_SUBNORMAL);
assert!(tiny != 0.0);
}
#[test]
fn next_float_subnormal() {
let second = next_float(MIN_SUBNORMAL);
// For subnormals, MIN_SUBNORMAL is the ULP
assert!(second != MIN_SUBNORMAL);
assert!(second > 0.0);
assert_eq!(second - MIN_SUBNORMAL, MIN_SUBNORMAL);
}
#[test]
fn next_float_inf() {
assert_eq!(next_float(f64::MAX), f64::INFINITY);
assert_eq!(next_float(f64::INFINITY), f64::INFINITY);
}
#[test]
fn next_prev_identity() {
for &x in &SOME_FLOATS {
assert_eq!(prev_float(next_float(x)), x);
assert_eq!(prev_float(prev_float(next_float(next_float(x)))), x);
assert_eq!(next_float(prev_float(x)), x);
assert_eq!(next_float(next_float(prev_float(prev_float(x)))), x);
}
}
#[test]
fn next_float_monotonic() {
let mut x = 0.49999999999999;
assert!(x < 0.5);
for _ in 0..200 {
let x1 = next_float(x);
assert!(x1 > x);
assert!(x1 - x < 1e-15, "next_float_monotonic: delta = {:?}", x1 - x);
x = x1;
}
assert!(x > 0.5);
}

1
src/libcoretest/num/mod.rs
View file

@ -30,6 +30,7 @@ mod u32;
mod u64;
mod flt2dec;
mod dec2flt;
/// Helper function for testing numeric operations
pub fn test_num<T>(ten: T, two: T) where