      1
      2
      3
      4
      5
      6
      7
      8
      9
     10
     11
     12
     13
     14
     15
     16
     17
     18
     19
     20
     21
     22
     23
     24
     25
     26
     27
     28
     29
     30
     31
     32
     33
     34
     35
     36
     37
     38
     39
     40
     41
     42
     43
     44
     45
     46
     47
     48
     49
     50
     51
     52
     53
     54
     55
     56
     57
     58
     59
     60
     61
     62
     63
     64
     65
     66
     67
     68
     69
     70
     71
     72
     73
     74
     75
     76
     77
     78
     79
     80
     81
     82
     83
     84
     85
     86
     87
     88
     89
     90
     91
     92
     93
     94
     95
     96
     97
     98
     99
    100
    101
    102
    103
    104
    105
    106
    107
    108
    109
    110
    111
    112
    113
    114
    115
    116
    117
    118
    119
    120
    121
    122
    123
    124
    125
    126
    127
    128
    129
    130
    131
    132
    133
    134
    135
    136
    137
    138
    139
    140
    141
    142
    143
    144
    145
    146
    147
    148
    149
    150
    151
    152
    153
    154
    155
    156
    157
    158
    159
    160
    161
    162
    163
    164
    165
    166
    167
    168
    169
    170
    171
    172
    173
    174
    175
    176
    177
    178
    179
    180
    181
    182
    183
    184
    185
    186
    187
    188
    189
    190
    191
    192
    193
    194
    195
    196
    197
    198
    199
    200
    201
    202
    203
    204
    205
    206
    207
    208
    209
    210
    211
    212
    213
    214
    215
    216
    217
    218
    219
    220
    221
    222
    223
    224
    225
    226
    227
    228
    229
    230
    231
    232
    233
    234
    235
    236
    237
    238
    239
    240
    241
    242
    243
    244
    245
    246
    247
    248
    249
    250
    251
    252
    253
    254
    255
    256
    257
    258
    259
    260
    261
    262
    263
    264
    265
    266
    267
    268
    269
    270
    271
    272
    273
    274
    275
    276
    277
    278
    279
    280
    281
    282
    283
    284
    285
    286
    287
    288
    289
    290
    291
    292
    293
    294
    295
    296
    297
    298
    299
    300
    301
    302
    303
    304
    305
    306
    307
    308
    309
    310
    311
    312
    313
    314
    315
    316
    317
    318
    319
    320
    321
    322
    323
    324
    325
    326
    327
    328
    329
    330
    331
    332
    333
    334
    335
    336
    337
    338
    339
    340
    341
    342
    343
    344
    345
    346
    347
    348
    349
    350

    //! 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}.

    #![doc(hidden)]
    #![unstable(
        feature = "dec2flt",
        reason = "internal routines only exposed for testing",
        issue = "none"
    )]

    use crate::fmt;
    use crate::str::FromStr;

    use self::num::digits_to_big;
    use self::parse::{parse_decimal, Decimal, ParseResult, Sign};
    use self::rawfp::RawFloat;

    mod algorithm;
    mod num;
    mod table;
    // These two have their own tests.
    pub mod parse;
    pub mod rawfp;

    macro_rules! from_str_float_impl {
        ($t:ty) => {
            #[stable(feature = "rust1", since = "1.0.0")]
            impl FromStr for $t {
                type Err = ParseFloatError;

                /// Converts a string in base 10 to a float.
                /// Accepts an optional decimal exponent.
                ///
                /// This function accepts strings such as
                ///
                /// * '3.14'
                /// * '-3.14'
                /// * '2.5E10', or equivalently, '2.5e10'
                /// * '2.5E-10'
                /// * '5.'
                /// * '.5', or, equivalently, '0.5'
                /// * 'inf', '-inf', 'NaN'
                ///
                /// Leading and trailing whitespace represent an error.
                ///
                /// # Grammar
                ///
                /// All strings that adhere to the following [EBNF] grammar
                /// will result in an [`Ok`] being returned:
                ///
                /// ```txt
                /// Float  ::= Sign? ( 'inf' | 'NaN' | Number )
                /// Number ::= ( Digit+ |
                ///              Digit+ '.' Digit* |
                ///              Digit* '.' Digit+ ) Exp?
                /// Exp    ::= [eE] Sign? Digit+
                /// Sign   ::= [+-]
                /// Digit  ::= [0-9]
                /// ```
                ///
                /// [EBNF]: https://www.w3.org/TR/REC-xml/#sec-notation
                ///
                /// # Known bugs
                ///
                /// In some situations, some strings that should create a valid float
                /// instead return an error. See [issue #31407] for details.
                ///
                /// [issue #31407]: https://github.com/rust-lang/rust/issues/31407
                ///
                /// # Arguments
                ///
                /// * src - A string
                ///
                /// # Return value
                ///
                /// `Err(ParseFloatError)` if the string did not represent a valid
                /// number. Otherwise, `Ok(n)` where `n` is the floating-point
                /// number represented by `src`.
                #[inline]
                fn from_str(src: &str) -> Result<Self, ParseFloatError> {
                    dec2flt(src)
                }
            }
        };
    }
    from_str_float_impl!(f32);
    from_str_float_impl!(f64);

    /// An error which can be returned when parsing a float.
    ///
    /// This error is used as the error type for the [`FromStr`] implementation
    /// for [`f32`] and [`f64`].
    ///
    /// [`FromStr`]: ../str/trait.FromStr.html
    /// [`f32`]: ../../std/primitive.f32.html
    /// [`f64`]: ../../std/primitive.f64.html
    #[derive(Debug, Clone, PartialEq, Eq)]
    #[stable(feature = "rust1", since = "1.0.0")]
    pub struct ParseFloatError {
        kind: FloatErrorKind,
    }

    #[derive(Debug, Clone, PartialEq, Eq)]
    enum FloatErrorKind {
        Empty,
        Invalid,
    }

    impl ParseFloatError {
        #[unstable(
            feature = "int_error_internals",
            reason = "available through Error trait and this method should \
                      not be exposed publicly",
            issue = "none"
        )]
        #[doc(hidden)]
        pub fn __description(&self) -> &str {
            match self.kind {
                FloatErrorKind::Empty => "cannot parse float from empty string",
                FloatErrorKind::Invalid => "invalid float literal",
            }
        }
    }

    #[stable(feature = "rust1", since = "1.0.0")]
    impl fmt::Display for ParseFloatError {
        fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
            self.__description().fmt(f)
        }
    }

    fn pfe_empty() -> ParseFloatError {
        ParseFloatError { kind: FloatErrorKind::Empty }
    }

    fn pfe_invalid() -> ParseFloatError {
        ParseFloatError { kind: FloatErrorKind::Invalid }
    }

    /// Splits a 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),
        }
    }

    /// Converts a decimal string into a floating point number.
    fn dec2flt<T: RawFloat>(s: &str) -> Result<T, ParseFloatError> {
        if s.is_empty() {
            return Err(pfe_empty());
        }
        let (sign, s) = extract_sign(s);
        let flt = match parse_decimal(s) {
            ParseResult::Valid(decimal) => convert(decimal)?,
            ParseResult::ShortcutToInf => T::INFINITY,
            ParseResult::ShortcutToZero => T::ZERO,
            ParseResult::Invalid => match s {
                "inf" => T::INFINITY,
                "NaN" => T::NAN,
                _ => {
                    return Err(pfe_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, ParseFloatError> {
        simplify(&mut decimal);
        if let Some(x) = trivial_cases(&decimal) {
            return Ok(x);
        }
        // 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, we'll crash, so we error out before getting too close (within 10^10).
        let upper_bound = bound_intermediate_digits(&decimal, e);
        if upper_bound > 375 {
            return Err(pfe_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 = upper_bound <= T::MAX_NORMAL_DIGITS as u64;
        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;
        }
    }

    /// Returns a quick-an-dirty upper bound on the size (log10) of the largest value that Algorithm R
    /// and Algorithm M will compute while working on the given decimal.
    fn bound_intermediate_digits(decimal: &Decimal<'_>, e: i64) -> u64 {
        // We don't need to worry too much about overflow here thanks to trivial_cases() and the
        // parser, which filter out the most extreme inputs for us.
        let f_len: u64 = decimal.integral.len() as u64 + decimal.fractional.len() as u64;
        if e >= 0 {
            // In the case e >= 0, both algorithms compute about `f * 10^e`. Algorithm R proceeds to
            // do some complicated calculations with this but we can ignore that for the upper bound
            // because it also reduces the fraction beforehand, so we have plenty of buffer there.
            f_len + (e as u64)
        } else {
            // If e < 0, Algorithm R does roughly the same thing, but Algorithm M differs:
            // It tries to find a positive number k such that `f << k / 10^e` is an in-range
            // significand. This will result in about `2^53 * f * 10^e` < `10^17 * f * 10^e`.
            // One input that triggers this is 0.33...33 (375 x 3).
            f_len + (e.abs() as u64) + 17
        }
    }

    /// Detects 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)). We don't need to worry too
        // much about overflow here because the input length is tiny (at least compared to 2^64) and
        // the parser already handles exponents whose absolute value is greater than 10^18
        // (which is still 10^19 short of 2^64).
        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
    }