Auto merge of #27307 - rkruppe:dec2flt, r=pnkfelix

Completely rewrite the conversion of decimal strings to `f64` and `f32`. The code is intended to be absolutely positively completely 100% accurate (when it doesn't give up). To the best of my knowledge, it achieves that goal. Any input that is not rejected is converted to the floating point number that is closest to the true value of the input. This includes overflow, subnormal numbers, and underflow to zero. In other words, the rounding error is less than or equal to 0.5 units in the last place. Half-way cases (exactly 0.5 ULP error) are handled with half-to-even rounding, also known as banker's rounding.

This code implements the algorithms from the paper [How to Read Floating Point Numbers Accurately][paper] by William D. Clinger, with extensions to handle underflow, overflow and subnormals, as well as some algorithmic optimizations.

# Correctness

With such a large amount of tricky code, many bugs are to be expected. Indeed tracking down the obscure causes of various rounding errors accounts for the bulk of the development time. Extensive tests (taking in the order of hours to run through to completion) are included in `src/etc/test-float-parse`: Though exhaustively testing all possible inputs is impossible, I've had good success with generating millions of instances from various "classes" of inputs. These tests take far too long to be run by @bors so contributors who touch this code need the discipline to run them. There are `#[test]`s, but they don't even cover every stupid mistake I made in course of writing this.

Another aspect is *integer* overflow. Extreme (or malicious) inputs could cause overflow both in the machine-sized integers used for bookkeeping throughout the algorithms (e.g., the decimal exponent) as well as the arbitrary-precision arithmetic. There is input validation to reject all such cases I know of, and I am quite sure nobody will *accidentally* cause this code to go out of range. Still, no guarantees.

# Limitations

Noticed the weasel words "(when it doesn't give up)" at the beginning? Some otherwise well-formed decimal strings are rejected because spelling out the value of the input requires too many digits, i.e., `digits * 10^abs(exp)` can't be stored in a bignum. This only applies if the value is not "obviously" zero or infinite, i.e., if you take a near-infinity or near-zero value and add many pointless fractional digits. At least with the algorithm used here, computing the precise value would require computing the full value as a fraction, which would overflow. The precise limit is `number_of_digits + abs(exp) > 375` but could be raised almost arbitrarily. In the future, another algorithm might lift this restriction entirely.

This should not be an issue for any realistic inputs. Still, the code does reject inputs that would result in a finite float when evaluated with unlimited precision. Some of these inputs are even regressions that the old code (mostly) handled, such as `0.333...333` with 400+ `3`s. Thus this might qualify as [breaking-change].

# Performance

Benchmarks results are... tolerable. Short numbers that hit the fast paths (`f64` multiplication or shortcuts to zero/inf) have performance in the same order of magnitude as the old code tens of nanoseconds. Numbers that are delegated to Algorithm Bellerophon (using floats with 64 bit significand, implemented in software) are slower, but not drastically so (couple hundred nanoseconds).

Numbers that need the AlgorithmM fallback (for `f64`, roughly everything below 1e-305 and above 1e305) take far, far longer, hundreds of microseconds. Note that my implementation is not quite as naive as the expository version in the paper (it needs one to four division instead of ~1000), but division is fundamentally pretty expensive and my implementation of it is extremely simple and slow.

All benchmarks run on a mediocre laptop with a i5-4200U CPU under light load.

# Binary size

Unfortunately the implementation needs to duplicate almost all code: Once for `f32` and once for `f64`. Before you ask, no, this cannot be avoided, at least not completely (but see the Future Work section). There's also a precomputed table of powers of ten, weighing in at about six kilobytes.

Running a stage1 `rustc` over a stand-alone program that simply parses pi to `f32` and `f64` and outputs both results reveals that the overhead vs. the old parsing code is about 44 KiB normally and about 28 KiB with LTO. It's presumably half of that + 3 KiB when only one of the two code paths is exercised.

| rustc options                 | old       | new       | delta         |
|---------------------------    |---------  |---------  |-----------    |
| [nothing]                     | 2588375   | 2633828   | 44.39 KiB     |
| -O                            | 2585211   | 2630688   | 44.41 KiB     |
| -O -C lto                     | 1026353   | 1054981   | 27.96 KiB     |
| -O -C lto -C link-args=-s     | 414208    | 442368    | 27.5 KiB      |

# Future Work

## Directory layout

The `dec2flt` code uses some types embedded deeply in the `flt2dec` module hierarchy, even though nothing about them it formatting-specific. They should be moved to a more conversion-direction-agnostic location at some point.

## Performance

It could be much better, especially for large inputs. Some low-hanging fruit has been picked but much more work could be done. Some specific ideas are jotted down in `FIXME`s all over the code.

## Binary size

One could try to compress the table further, though I am skeptical. Another avenue would be reducing the code duplication from basically everything being generic over `T: RawFloat`. Perhaps one can reduce the magnitude of the duplication by pushing the parts that don't need to know the target type into separate functions, but this is finicky and probably makes some code read less naturally.

## Other bases

This PR leaves `f{32,64}::from_str_radix` alone. It only replaces `FromStr` (and thus `.parse()`). I am convinced that `from_str_radix` should not exist, and have proposed its [deprecation and speedy removal][deprecate-radix]. Whatever the outcome of that discussion, it is independent from, and out of scope for, this PR.

Fixes #24557
Fixes #14353

r? @pnkfelix

cc @lifthrasiir @huonw 

[paper]: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.45.4152
[deprecate-radix]: https://internals.rust-lang.org/t/deprecate-f-32-64-from-str-radix/2405

3 files changed, 148 insertions, 13 deletions

diff --git a/src/libcore/num/flt2dec/bignum.rs b/src/libcore/num/flt2dec/bignum.rs
index 1e39c53f9e0..ee1f6ffdd0a 100644
--- a/src/libcore/num/flt2dec/bignum.rs
+++ b/src/libcore/num/flt2dec/bignum.rs
@@ -11,9 +11,9 @@
 //! Custom arbitrary-precision number (bignum) implementation.
 //!
 //! This is designed to avoid the heap allocation at expense of stack memory.
-//! The most used bignum type, `Big32x36`, is limited by 32 × 36 = 1,152 bits
-//! and will take at most 152 bytes of stack memory. This is (barely) enough
-//! for handling all possible finite `f64` values.
+//! The most used bignum type, `Big32x40`, is limited by 32 × 40 = 1,280 bits
+//! and will take at most 160 bytes of stack memory. This is more than enough
+//! for round-tripping all possible finite `f64` values.
 //!
 //! In principle it is possible to have multiple bignum types for different
 //! inputs, but we don't do so to avoid the code bloat. Each bignum is still
@@ -92,6 +92,14 @@ impl_full_ops! {
 //  u64: add(intrinsics::u64_add_with_overflow), mul/div(u128); // see RFC #521 for enabling this.
 }
 
+/// Table of powers of 5 representable in digits. Specifically, the largest {u8, u16, u32} value
+/// that's a power of five, plus the corresponding exponent. Used in `mul_pow5`.
+const SMALL_POW5: [(u64, usize); 3] = [
+    (125, 3),
+    (15625, 6),
+    (1_220_703_125, 13),
+];
+
 macro_rules! define_bignum {
     ($name:ident: type=$ty:ty, n=$n:expr) => (
         /// Stack-allocated arbitrary-precision (up to certain limit) integer.
@@ -135,9 +143,52 @@ macro_rules! define_bignum {
                 $name { size: sz, base: base }
             }
 
+            /// Return the internal digits as a slice `[a, b, c, ...]` such that the numeric
+            /// value is `a + b * 2^W + c * 2^(2W) + ...` where `W` is the number of bits in
+            /// the digit type.
+            pub fn digits(&self) -> &[$ty] {
+                &self.base[..self.size]
+            }
+
+            /// Return the `i`-th bit where bit 0 is the least significant one.
+            /// In other words, the bit with weight `2^i`.
+            pub fn get_bit(&self, i: usize) -> u8 {
+                use mem;
+
+                let digitbits = mem::size_of::<$ty>() * 8;
+                let d = i / digitbits;
+                let b = i % digitbits;
+                ((self.base[d] >> b) & 1) as u8
+            }
+
             /// Returns true if the bignum is zero.
             pub fn is_zero(&self) -> bool {
-                self.base[..self.size].iter().all(|&v| v == 0)
+                self.digits().iter().all(|&v| v == 0)
+            }
+
+            /// Returns the number of bits necessary to represent this value. Note that zero
+            /// is considered to need 0 bits.
+            pub fn bit_length(&self) -> usize {
+                use mem;
+
+                // Skip over the most significant digits which are zero.
+                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;
+                }
+                i + 1
             }
 
             /// Adds `other` to itself and returns its own mutable reference.
@@ -160,6 +211,24 @@ macro_rules! define_bignum {
                 self
             }
 
+            pub fn add_small<'a>(&'a mut self, other: $ty) -> &'a mut $name {
+                use num::flt2dec::bignum::FullOps;
+
+                let (mut carry, v) = self.base[0].full_add(other, false);
+                self.base[0] = v;
+                let mut i = 1;
+                while carry {
+                    let (c, v) = self.base[i].full_add(0, carry);
+                    self.base[i] = v;
+                    carry = c;
+                    i += 1;
+                }
+                if i > self.size {
+                    self.size = i;
+                }
+                self
+            }
+
             /// Subtracts `other` from itself and returns its own mutable reference.
             pub fn sub<'a>(&'a mut self, other: &$name) -> &'a mut $name {
                 use cmp;
@@ -238,6 +307,34 @@ macro_rules! define_bignum {
                 self
             }
 
+            /// Multiplies itself by `5^e` and returns its own mutable reference.
+            pub fn mul_pow5<'a>(&'a mut self, mut e: usize) -> &'a mut $name {
+                use mem;
+                use num::flt2dec::bignum::SMALL_POW5;
+
+                // There are exactly n trailing zeros on 2^n, and the only relevant digit sizes
+                // are consecutive powers of two, so this is well suited index for the table.
+                let table_index = mem::size_of::<$ty>().trailing_zeros() as usize;
+                let (small_power, small_e) = SMALL_POW5[table_index];
+                let small_power = small_power as $ty;
+
+                // Multiply with the largest single-digit power as long as possible ...
+                while e >= small_e {
+                    self.mul_small(small_power);
+                    e -= small_e;
+                }
+
+                // ... then finish off the remainder.
+                let mut rest_power = 1;
+                for _ in 0..e {
+                    rest_power *= 5;
+                }
+                self.mul_small(rest_power);
+
+                self
+            }
+
+
             /// Multiplies itself by a number described by `other[0] + other[1] * 2^W +
             /// other[2] * 2^(2W) + ...` (where `W` is the number of bits in the digit type)
             /// and returns its own mutable reference.
@@ -269,9 +366,9 @@ macro_rules! define_bignum {
 
                 let mut ret = [0; $n];
                 let retsz = if self.size < other.len() {
-                    mul_inner(&mut ret, &self.base[..self.size], other)
+                    mul_inner(&mut ret, &self.digits(), other)
                 } else {
-                    mul_inner(&mut ret, other, &self.base[..self.size])
+                    mul_inner(&mut ret, other, &self.digits())
                 };
                 self.base = ret;
                 self.size = retsz;
@@ -294,6 +391,45 @@ macro_rules! define_bignum {
                 }
                 (self, borrow)
             }
+
+            /// Divide self by another bignum, overwriting `q` with the quotient and `r` with the
+            /// remainder.
+            pub fn div_rem(&self, d: &$name, q: &mut $name, r: &mut $name) {
+                use mem;
+
+                // Stupid slow base-2 long division taken from
+                // https://en.wikipedia.org/wiki/Division_algorithm
+                // FIXME use a greater base ($ty) for the long division.
+                assert!(!d.is_zero());
+                let digitbits = mem::size_of::<$ty>() * 8;
+                for digit in &mut q.base[..] {
+                    *digit = 0;
+                }
+                for digit in &mut r.base[..] {
+                    *digit = 0;
+                }
+                r.size = d.size;
+                q.size = 1;
+                let mut q_is_zero = true;
+                let end = self.bit_length();
+                for i in (0..end).rev() {
+                    r.mul_pow2(1);
+                    r.base[0] |= self.get_bit(i) as $ty;
+                    if &*r >= d {
+                        r.sub(d);
+                        // Set bit `i` of q to 1.
+                        let digit_idx = i / digitbits;
+                        let bit_idx = i % digitbits;
+                        if q_is_zero {
+                            q.size = digit_idx + 1;
+                            q_is_zero = false;
+                        }
+                        q.base[digit_idx] |= 1 << bit_idx;
+                    }
+                }
+                debug_assert!(q.base[q.size..].iter().all(|&d| d == 0));
+                debug_assert!(r.base[r.size..].iter().all(|&d| d == 0));
+            }
         }
 
         impl ::cmp::PartialEq for $name {
@@ -344,10 +480,10 @@ macro_rules! define_bignum {
     )
 }
 
-/// The digit type for `Big32x36`.
+/// The digit type for `Big32x40`.
 pub type Digit32 = u32;
 
-define_bignum!(Big32x36: type=Digit32, n=36);
+define_bignum!(Big32x40: type=Digit32, n=40);
 
 // this one is used for testing only.
 #[doc(hidden)]
@@ -355,4 +491,3 @@ pub mod tests {
     use prelude::v1::*;
     define_bignum!(Big8x3: type=u8, n=3);
 }
-
diff --git a/src/libcore/num/flt2dec/strategy/dragon.rs b/src/libcore/num/flt2dec/strategy/dragon.rs
index b03286ddd0d..cdc23c45fa0 100644
--- a/src/libcore/num/flt2dec/strategy/dragon.rs
+++ b/src/libcore/num/flt2dec/strategy/dragon.rs
@@ -23,7 +23,7 @@ use cmp::Ordering;
 use num::flt2dec::{Decoded, MAX_SIG_DIGITS, round_up};
 use num::flt2dec::estimator::estimate_scaling_factor;
 use num::flt2dec::bignum::Digit32 as Digit;
-use num::flt2dec::bignum::Big32x36 as Big;
+use num::flt2dec::bignum::Big32x40 as Big;
 
 static POW10: [Digit; 10] = [1, 10, 100, 1000, 10000, 100000,
                              1000000, 10000000, 100000000, 1000000000];
diff --git a/src/libcore/num/flt2dec/strategy/grisu.rs b/src/libcore/num/flt2dec/strategy/grisu.rs
index 390920a354c..52eafcec184 100644
--- a/src/libcore/num/flt2dec/strategy/grisu.rs
+++ b/src/libcore/num/flt2dec/strategy/grisu.rs
@@ -34,7 +34,7 @@ pub struct Fp {
 
 impl Fp {
     /// Returns a correctly rounded product of itself and `other`.
-    fn mul(&self, other: &Fp) -> Fp {
+    pub fn mul(&self, other: &Fp) -> Fp {
         const MASK: u64 = 0xffffffff;
         let a = self.f >> 32;
         let b = self.f & MASK;
@@ -51,7 +51,7 @@ impl Fp {
     }
 
     /// Normalizes itself so that the resulting mantissa is at least `2^63`.
-    fn normalize(&self) -> Fp {
+    pub fn normalize(&self) -> Fp {
         let mut f = self.f;
         let mut e = self.e;
         if f >> (64 - 32) == 0 { f <<= 32; e -= 32; }
@@ -66,7 +66,7 @@ impl Fp {
 
     /// Normalizes itself to have the shared exponent.
     /// It can only decrease the exponent (and thus increase the mantissa).
-    fn normalize_to(&self, e: i16) -> Fp {
+    pub fn normalize_to(&self, e: i16) -> Fp {
         let edelta = self.e - e;
         assert!(edelta >= 0);
         let edelta = edelta as usize;