]> git.proxmox.com Git - rustc.git/blame - src/libcore/num/dec2flt/mod.rs
Imported Upstream version 1.9.0+dfsg1
[rustc.git] / src / libcore / num / dec2flt / mod.rs
CommitLineData
e9174d1e
SL
1// Copyright 2015 The Rust Project Developers. See the COPYRIGHT
2// file at the top-level directory of this distribution and at
3// http://rust-lang.org/COPYRIGHT.
4//
5// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
6// http://www.apache.org/licenses/LICENSE-2.0> or the MIT license
7// <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your
8// option. This file may not be copied, modified, or distributed
9// except according to those terms.
10
11//! Converting decimal strings into IEEE 754 binary floating point numbers.
12//!
13//! # Problem statement
14//!
15//! We are given a decimal string such as `12.34e56`. This string consists of integral (`12`),
16//! fractional (`45`), and exponent (`56`) parts. All parts are optional and interpreted as zero
17//! when missing.
18//!
19//! We seek the IEEE 754 floating point number that is closest to the exact value of the decimal
20//! string. It is well-known that many decimal strings do not have terminating representations in
21//! base two, so we round to 0.5 units in the last place (in other words, as well as possible).
22//! Ties, decimal values exactly half-way between two consecutive floats, are resolved with the
23//! half-to-even strategy, also known as banker's rounding.
24//!
25//! Needless to say, this is quite hard, both in terms of implementation complexity and in terms
26//! of CPU cycles taken.
27//!
28//! # Implementation
29//!
30//! First, we ignore signs. Or rather, we remove it at the very beginning of the conversion
31//! process and re-apply it at the very end. This is correct in all edge cases since IEEE
32//! floats are symmetric around zero, negating one simply flips the first bit.
33//!
34//! Then we remove the decimal point by adjusting the exponent: Conceptually, `12.34e56` turns
35//! into `1234e54`, which we describe with a positive integer `f = 1234` and an integer `e = 54`.
36//! The `(f, e)` representation is used by almost all code past the parsing stage.
37//!
38//! We then try a long chain of progressively more general and expensive special cases using
39//! machine-sized integers and small, fixed-sized floating point numbers (first `f32`/`f64`, then
40//! a type with 64 bit significand, `Fp`). When all these fail, we bite the bullet and resort to a
41//! simple but very slow algorithm that involved computing `f * 10^e` fully and doing an iterative
42//! search for the best approximation.
43//!
44//! Primarily, this module and its children implement the algorithms described in:
45//! "How to Read Floating Point Numbers Accurately" by William D. Clinger,
46//! available online: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.45.4152
47//!
48//! In addition, there are numerous helper functions that are used in the paper but not available
49//! in Rust (or at least in core). Our version is additionally complicated by the need to handle
50//! overflow and underflow and the desire to handle subnormal numbers. Bellerophon and
51//! Algorithm R have trouble with overflow, subnormals, and underflow. We conservatively switch to
52//! Algorithm M (with the modifications described in section 8 of the paper) well before the
53//! inputs get into the critical region.
54//!
55//! Another aspect that needs attention is the ``RawFloat`` trait by which almost all functions
56//! are parametrized. One might think that it's enough to parse to `f64` and cast the result to
57//! `f32`. Unfortunately this is not the world we live in, and this has nothing to do with using
58//! base two or half-to-even rounding.
59//!
60//! Consider for example two types `d2` and `d4` representing a decimal type with two decimal
61//! digits and four decimal digits each and take "0.01499" as input. Let's use half-up rounding.
62//! Going directly to two decimal digits gives `0.01`, but if we round to four digits first,
63//! we get `0.0150`, which is then rounded up to `0.02`. The same principle applies to other
64//! operations as well, if you want 0.5 ULP accuracy you need to do *everything* in full precision
65//! and round *exactly once, at the end*, by considering all truncated bits at once.
66//!
67//! FIXME Although some code duplication is necessary, perhaps parts of the code could be shuffled
68//! around such that less code is duplicated. Large parts of the algorithms are independent of the
69//! float type to output, or only needs access to a few constants, which could be passed in as
70//! parameters.
71//!
72//! # Other
73//!
74//! The conversion should *never* panic. There are assertions and explicit panics in the code,
75//! but they should never be triggered and only serve as internal sanity checks. Any panics should
76//! be considered a bug.
77//!
78//! There are unit tests but they are woefully inadequate at ensuring correctness, they only cover
79//! a small percentage of possible errors. Far more extensive tests are located in the directory
80//! `src/etc/test-float-parse` as a Python script.
81//!
82//! A note on integer overflow: Many parts of this file perform arithmetic with the decimal
83//! exponent `e`. Primarily, we shift the decimal point around: Before the first decimal digit,
84//! after the last decimal digit, and so on. This could overflow if done carelessly. We rely on
85//! the parsing submodule to only hand out sufficiently small exponents, where "sufficient" means
86//! "such that the exponent +/- the number of decimal digits fits into a 64 bit integer".
87//! Larger exponents are accepted, but we don't do arithmetic with them, they are immediately
88//! turned into {positive,negative} {zero,infinity}.
e9174d1e
SL
89
90#![doc(hidden)]
91#![unstable(feature = "dec2flt",
92 reason = "internal routines only exposed for testing",
93 issue = "0")]
94
95use prelude::v1::*;
96use fmt;
97use str::FromStr;
98
9cc50fc6 99use self::parse::{parse_decimal, Decimal, Sign, ParseResult};
e9174d1e
SL
100use self::num::digits_to_big;
101use self::rawfp::RawFloat;
102
103mod algorithm;
104mod table;
105mod num;
106// These two have their own tests.
107pub mod rawfp;
108pub mod parse;
109
110macro_rules! from_str_float_impl {
9cc50fc6 111 ($t:ty) => {
e9174d1e
SL
112 #[stable(feature = "rust1", since = "1.0.0")]
113 impl FromStr for $t {
114 type Err = ParseFloatError;
115
116 /// Converts a string in base 10 to a float.
117 /// Accepts an optional decimal exponent.
118 ///
119 /// This function accepts strings such as
120 ///
121 /// * '3.14'
122 /// * '-3.14'
123 /// * '2.5E10', or equivalently, '2.5e10'
124 /// * '2.5E-10'
125 /// * '.' (understood as 0)
126 /// * '5.'
127 /// * '.5', or, equivalently, '0.5'
128 /// * 'inf', '-inf', 'NaN'
129 ///
130 /// Leading and trailing whitespace represent an error.
131 ///
132 /// # Arguments
133 ///
134 /// * src - A string
135 ///
136 /// # Return value
137 ///
138 /// `Err(ParseFloatError)` if the string did not represent a valid
139 /// number. Otherwise, `Ok(n)` where `n` is the floating-point
140 /// number represented by `src`.
141 #[inline]
142 fn from_str(src: &str) -> Result<Self, ParseFloatError> {
143 dec2flt(src)
144 }
145 }
146 }
147}
9cc50fc6
SL
148from_str_float_impl!(f32);
149from_str_float_impl!(f64);
e9174d1e
SL
150
151/// An error which can be returned when parsing a float.
7453a54e
SL
152///
153/// This error is used as the error type for the [`FromStr`] implementation
154/// for [`f32`] and [`f64`].
155///
156/// [`FromStr`]: ../str/trait.FromStr.html
54a0048b
SL
157/// [`f32`]: ../../std/primitive.f32.html
158/// [`f64`]: ../../std/primitive.f64.html
e9174d1e
SL
159#[derive(Debug, Clone, PartialEq)]
160#[stable(feature = "rust1", since = "1.0.0")]
161pub struct ParseFloatError {
162 kind: FloatErrorKind
163}
164
165#[derive(Debug, Clone, PartialEq)]
166enum FloatErrorKind {
167 Empty,
168 Invalid,
169}
170
171impl ParseFloatError {
172 #[unstable(feature = "int_error_internals",
173 reason = "available through Error trait and this method should \
174 not be exposed publicly",
175 issue = "0")]
176 #[doc(hidden)]
177 pub fn __description(&self) -> &str {
178 match self.kind {
179 FloatErrorKind::Empty => "cannot parse float from empty string",
180 FloatErrorKind::Invalid => "invalid float literal",
181 }
182 }
183}
184
185#[stable(feature = "rust1", since = "1.0.0")]
186impl fmt::Display for ParseFloatError {
187 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
188 self.__description().fmt(f)
189 }
190}
191
9cc50fc6 192fn pfe_empty() -> ParseFloatError {
e9174d1e
SL
193 ParseFloatError { kind: FloatErrorKind::Empty }
194}
195
9cc50fc6 196fn pfe_invalid() -> ParseFloatError {
e9174d1e
SL
197 ParseFloatError { kind: FloatErrorKind::Invalid }
198}
199
200/// Split decimal string into sign and the rest, without inspecting or validating the rest.
201fn extract_sign(s: &str) -> (Sign, &str) {
202 match s.as_bytes()[0] {
203 b'+' => (Sign::Positive, &s[1..]),
204 b'-' => (Sign::Negative, &s[1..]),
205 // If the string is invalid, we never use the sign, so we don't need to validate here.
206 _ => (Sign::Positive, s),
207 }
208}
209
210/// Convert a decimal string into a floating point number.
211fn dec2flt<T: RawFloat>(s: &str) -> Result<T, ParseFloatError> {
212 if s.is_empty() {
213 return Err(pfe_empty())
214 }
215 let (sign, s) = extract_sign(s);
216 let flt = match parse_decimal(s) {
54a0048b 217 ParseResult::Valid(decimal) => convert(decimal)?,
9cc50fc6
SL
218 ParseResult::ShortcutToInf => T::infinity(),
219 ParseResult::ShortcutToZero => T::zero(),
e9174d1e
SL
220 ParseResult::Invalid => match s {
221 "inf" => T::infinity(),
222 "NaN" => T::nan(),
223 _ => { return Err(pfe_invalid()); }
224 }
225 };
226
227 match sign {
228 Sign::Positive => Ok(flt),
229 Sign::Negative => Ok(-flt),
230 }
231}
232
233/// The main workhorse for the decimal-to-float conversion: Orchestrate all the preprocessing
234/// and figure out which algorithm should do the actual conversion.
235fn convert<T: RawFloat>(mut decimal: Decimal) -> Result<T, ParseFloatError> {
236 simplify(&mut decimal);
237 if let Some(x) = trivial_cases(&decimal) {
238 return Ok(x);
239 }
e9174d1e
SL
240 // Remove/shift out the decimal point.
241 let e = decimal.exp - decimal.fractional.len() as i64;
242 if let Some(x) = algorithm::fast_path(decimal.integral, decimal.fractional, e) {
243 return Ok(x);
244 }
245 // Big32x40 is limited to 1280 bits, which translates to about 385 decimal digits.
7453a54e
SL
246 // If we exceed this, we'll crash, so we error out before getting too close (within 10^10).
247 let upper_bound = bound_intermediate_digits(&decimal, e);
248 if upper_bound > 375 {
e9174d1e
SL
249 return Err(pfe_invalid());
250 }
251 let f = digits_to_big(decimal.integral, decimal.fractional);
252
253 // Now the exponent certainly fits in 16 bit, which is used throughout the main algorithms.
254 let e = e as i16;
255 // FIXME These bounds are rather conservative. A more careful analysis of the failure modes
256 // of Bellerophon could allow using it in more cases for a massive speed up.
257 let exponent_in_range = table::MIN_E <= e && e <= table::MAX_E;
7453a54e 258 let value_in_range = upper_bound <= T::max_normal_digits() as u64;
e9174d1e
SL
259 if exponent_in_range && value_in_range {
260 Ok(algorithm::bellerophon(&f, e))
261 } else {
262 Ok(algorithm::algorithm_m(&f, e))
263 }
264}
265
266// As written, this optimizes badly (see #27130, though it refers to an old version of the code).
267// `inline(always)` is a workaround for that. There are only two call sites overall and it doesn't
268// make code size worse.
269
270/// Strip zeros where possible, even when this requires changing the exponent
271#[inline(always)]
272fn simplify(decimal: &mut Decimal) {
273 let is_zero = &|&&d: &&u8| -> bool { d == b'0' };
274 // Trimming these zeros does not change anything but may enable the fast path (< 15 digits).
275 let leading_zeros = decimal.integral.iter().take_while(is_zero).count();
276 decimal.integral = &decimal.integral[leading_zeros..];
277 let trailing_zeros = decimal.fractional.iter().rev().take_while(is_zero).count();
278 let end = decimal.fractional.len() - trailing_zeros;
279 decimal.fractional = &decimal.fractional[..end];
280 // Simplify numbers of the form 0.0...x and x...0.0, adjusting the exponent accordingly.
281 // This may not always be a win (possibly pushes some numbers out of the fast path), but it
282 // simplifies other parts significantly (notably, approximating the magnitude of the value).
283 if decimal.integral.is_empty() {
284 let leading_zeros = decimal.fractional.iter().take_while(is_zero).count();
285 decimal.fractional = &decimal.fractional[leading_zeros..];
286 decimal.exp -= leading_zeros as i64;
287 } else if decimal.fractional.is_empty() {
288 let trailing_zeros = decimal.integral.iter().rev().take_while(is_zero).count();
289 let end = decimal.integral.len() - trailing_zeros;
290 decimal.integral = &decimal.integral[..end];
291 decimal.exp += trailing_zeros as i64;
292 }
293}
294
7453a54e
SL
295/// Quick and dirty upper bound on the size (log10) of the largest value that Algorithm R and
296/// Algorithm M will compute while working on the given decimal.
297fn bound_intermediate_digits(decimal: &Decimal, e: i64) -> u64 {
298 // We don't need to worry too much about overflow here thanks to trivial_cases() and the
299 // parser, which filter out the most extreme inputs for us.
300 let f_len: u64 = decimal.integral.len() as u64 + decimal.fractional.len() as u64;
301 if e >= 0 {
302 // In the case e >= 0, both algorithms compute about `f * 10^e`. Algorithm R proceeds to
303 // do some complicated calculations with this but we can ignore that for the upper bound
304 // because it also reduces the fraction beforehand, so we have plenty of buffer there.
305 f_len + (e as u64)
306 } else {
307 // If e < 0, Algorithm R does roughly the same thing, but Algorithm M differs:
308 // It tries to find a positive number k such that `f << k / 10^e` is an in-range
309 // significand. This will result in about `2^53 * f * 10^e` < `10^17 * f * 10^e`.
310 // One input that triggers this is 0.33...33 (375 x 3).
311 f_len + (e.abs() as u64) + 17
312 }
313}
314
e9174d1e
SL
315/// Detect obvious overflows and underflows without even looking at the decimal digits.
316fn trivial_cases<T: RawFloat>(decimal: &Decimal) -> Option<T> {
317 // There were zeros but they were stripped by simplify()
318 if decimal.integral.is_empty() && decimal.fractional.is_empty() {
319 return Some(T::zero());
320 }
7453a54e
SL
321 // This is a crude approximation of ceil(log10(the real value)). We don't need to worry too
322 // much about overflow here because the input length is tiny (at least compared to 2^64) and
323 // the parser already handles exponents whose absolute value is greater than 10^18
324 // (which is still 10^19 short of 2^64).
e9174d1e
SL
325 let max_place = decimal.exp + decimal.integral.len() as i64;
326 if max_place > T::inf_cutoff() {
327 return Some(T::infinity());
328 } else if max_place < T::zero_cutoff() {
329 return Some(T::zero());
330 }
331 None
332}