From b55806ca8ff97ee89f77f9c784619ace3034c32c Mon Sep 17 00:00:00 2001 From: Robin Kruppe Date: Sat, 27 Jun 2015 21:59:31 +0200 Subject: [PATCH 1/7] Make core::num::dec2flt::strategy::grisu::Fp methods public. The intent is to allow decimal-to-float parsing to use Fp in its fast path. That code is added in a later commit. --- src/libcore/num/flt2dec/strategy/grisu.rs | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/libcore/num/flt2dec/strategy/grisu.rs b/src/libcore/num/flt2dec/strategy/grisu.rs index 390920a354cfb..52eafcec184d2 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; From 7ff10209aa9b8da6d6d4ceea0161757048126d2d Mon Sep 17 00:00:00 2001 From: Robin Kruppe Date: Thu, 23 Jul 2015 22:18:44 +0200 Subject: [PATCH 2/7] Enlarge Bignum type from 1152 to 1280 bits. This is necessary for decimal-to-float code (in a later commit) to handle inputs such as 4.9406564584124654e-324 (the smallest subnormal f64). According to the benchmarks for flt2dec::dragon, this does not affect performance measurably. It probably uses slightly more stack space though. --- src/libcore/num/flt2dec/bignum.rs | 10 +++++----- src/libcore/num/flt2dec/strategy/dragon.rs | 2 +- src/libcoretest/num/flt2dec/strategy/dragon.rs | 2 +- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/libcore/num/flt2dec/bignum.rs b/src/libcore/num/flt2dec/bignum.rs index 1e39c53f9e06a..33c5a3ded980f 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 formatting and parsing 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 @@ -344,10 +344,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)] diff --git a/src/libcore/num/flt2dec/strategy/dragon.rs b/src/libcore/num/flt2dec/strategy/dragon.rs index b03286ddd0dc2..cdc23c45fa0b6 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/libcoretest/num/flt2dec/strategy/dragon.rs b/src/libcoretest/num/flt2dec/strategy/dragon.rs index f2397f6b48037..26bb5b9d9c5bc 100644 --- a/src/libcoretest/num/flt2dec/strategy/dragon.rs +++ b/src/libcoretest/num/flt2dec/strategy/dragon.rs @@ -12,7 +12,7 @@ use std::prelude::v1::*; use std::{i16, f64}; use super::super::*; use core::num::flt2dec::*; -use core::num::flt2dec::bignum::Big32x36 as Big; +use core::num::flt2dec::bignum::Big32x40 as Big; use core::num::flt2dec::strategy::dragon::*; #[test] From 7ebd7f3b9a7ebc020663a13b29b1e50446b3c262 Mon Sep 17 00:00:00 2001 From: Robin Kruppe Date: Wed, 1 Jul 2015 21:38:43 +0000 Subject: [PATCH 3/7] Add various methods to Bignum: - Exposing digits and individual bits - Counting the number of bits - Add small (digit-sized) values - Multiplication by power of 5 - Division with remainder All are necessary for decimal to floating point conversions. All but the most trivial ones come with tests. --- src/libcore/num/flt2dec/bignum.rs | 146 +++++++++++++++++++++++++- src/libcoretest/num/flt2dec/bignum.rs | 89 ++++++++++++++++ 2 files changed, 230 insertions(+), 5 deletions(-) diff --git a/src/libcore/num/flt2dec/bignum.rs b/src/libcore/num/flt2dec/bignum.rs index 33c5a3ded980f..15980f925fa03 100644 --- a/src/libcore/num/flt2dec/bignum.rs +++ b/src/libcore/num/flt2dec/bignum.rs @@ -13,7 +13,7 @@ //! This is designed to avoid the heap allocation at expense of stack memory. //! 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 formatting and parsing all possible finite `f64` values. +//! 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,53 @@ 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; + + 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; + } + }; + // This could be optimized with leading_zeros() and bit shifts, but that's + // probably not worth the hassle. + 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 +212,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 +308,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 +367,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 +392,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 { @@ -355,4 +492,3 @@ pub mod tests { use prelude::v1::*; define_bignum!(Big8x3: type=u8, n=3); } - diff --git a/src/libcoretest/num/flt2dec/bignum.rs b/src/libcoretest/num/flt2dec/bignum.rs index 09a1ed41dadd1..31065b2898f82 100644 --- a/src/libcoretest/num/flt2dec/bignum.rs +++ b/src/libcoretest/num/flt2dec/bignum.rs @@ -39,6 +39,23 @@ fn test_add_overflow_2() { Big::from_u64(0xffffff).add(&Big::from_small(1)); } +#[test] +fn test_add_small() { + assert_eq!(*Big::from_small(3).add_small(4), Big::from_small(7)); + assert_eq!(*Big::from_small(3).add_small(0), Big::from_small(3)); + assert_eq!(*Big::from_small(0).add_small(3), Big::from_small(3)); + assert_eq!(*Big::from_small(7).add_small(250), Big::from_u64(257)); + assert_eq!(*Big::from_u64(0x7fff).add_small(1), Big::from_u64(0x8000)); + assert_eq!(*Big::from_u64(0x2ffe).add_small(0x35), Big::from_u64(0x3033)); + assert_eq!(*Big::from_small(0xdc).add_small(0x89), Big::from_u64(0x165)); +} + +#[test] +#[should_panic] +fn test_add_small_overflow() { + Big::from_u64(0xffffff).add_small(1); +} + #[test] fn test_sub() { assert_eq!(*Big::from_small(7).sub(&Big::from_small(4)), Big::from_small(3)); @@ -97,6 +114,30 @@ fn test_mul_pow2_overflow_2() { Big::from_u64(0x123).mul_pow2(16); } +#[test] +fn test_mul_pow5() { + assert_eq!(*Big::from_small(42).mul_pow5(0), Big::from_small(42)); + assert_eq!(*Big::from_small(1).mul_pow5(2), Big::from_small(25)); + assert_eq!(*Big::from_small(1).mul_pow5(4), Big::from_u64(25 * 25)); + assert_eq!(*Big::from_small(4).mul_pow5(3), Big::from_u64(500)); + assert_eq!(*Big::from_small(140).mul_pow5(2), Big::from_u64(25 * 140)); + assert_eq!(*Big::from_small(25).mul_pow5(1), Big::from_small(125)); + assert_eq!(*Big::from_small(125).mul_pow5(7), Big::from_u64(9765625)); + assert_eq!(*Big::from_small(0).mul_pow5(127), Big::from_small(0)); +} + +#[test] +#[should_panic] +fn test_mul_pow5_overflow_1() { + Big::from_small(1).mul_pow5(12); +} + +#[test] +#[should_panic] +fn test_mul_pow5_overflow_2() { + Big::from_small(230).mul_pow5(8); +} + #[test] fn test_mul_digits() { assert_eq!(*Big::from_small(3).mul_digits(&[5]), Big::from_small(15)); @@ -132,6 +173,25 @@ fn test_div_rem_small() { (Big::from_u64(0x10000 / 123), (0x10000u64 % 123) as u8)); } +#[test] +fn test_div_rem() { + fn div_rem(n: u64, d: u64) -> (Big, Big) { + let mut q = Big::from_small(42); + let mut r = Big::from_small(42); + Big::from_u64(n).div_rem(&Big::from_u64(d), &mut q, &mut r); + (q, r) + } + assert_eq!(div_rem(1, 1), (Big::from_small(1), Big::from_small(0))); + assert_eq!(div_rem(4, 3), (Big::from_small(1), Big::from_small(1))); + assert_eq!(div_rem(1, 7), (Big::from_small(0), Big::from_small(1))); + assert_eq!(div_rem(45, 9), (Big::from_small(5), Big::from_small(0))); + assert_eq!(div_rem(103, 9), (Big::from_small(11), Big::from_small(4))); + assert_eq!(div_rem(123456, 77), (Big::from_u64(1603), Big::from_small(25))); + assert_eq!(div_rem(0xffff, 1), (Big::from_u64(0xffff), Big::from_small(0))); + assert_eq!(div_rem(0xeeee, 0xffff), (Big::from_small(0), Big::from_u64(0xeeee))); + assert_eq!(div_rem(2_000_000, 2), (Big::from_u64(1_000_000), Big::from_u64(0))); +} + #[test] fn test_is_zero() { assert!(Big::from_small(0).is_zero()); @@ -141,6 +201,35 @@ fn test_is_zero() { assert!(Big::from_u64(0xffffff).sub(&Big::from_u64(0xffffff)).is_zero()); } +#[test] +fn test_get_bit() { + let x = Big::from_small(0b1101); + assert_eq!(x.get_bit(0), 1); + assert_eq!(x.get_bit(1), 0); + assert_eq!(x.get_bit(2), 1); + assert_eq!(x.get_bit(3), 1); + let y = Big::from_u64(1 << 15); + assert_eq!(y.get_bit(14), 0); + assert_eq!(y.get_bit(15), 1); + assert_eq!(y.get_bit(16), 0); +} + +#[test] +#[should_panic] +fn test_get_bit_out_of_range() { + Big::from_small(42).get_bit(24); +} + +#[test] +fn test_bit_length() { + assert_eq!(Big::from_small(0).bit_length(), 0); + assert_eq!(Big::from_small(1).bit_length(), 1); + assert_eq!(Big::from_small(5).bit_length(), 3); + assert_eq!(Big::from_small(0x18).bit_length(), 5); + assert_eq!(Big::from_u64(0x4073).bit_length(), 15); + assert_eq!(Big::from_u64(0xffffff).bit_length(), 24); +} + #[test] fn test_ord() { assert!(Big::from_u64(0) < Big::from_u64(0xffffff)); From b7e39a1c2dd24fd4110c22c70cad254365b0ffd3 Mon Sep 17 00:00:00 2001 From: Robin Kruppe Date: Sun, 26 Jul 2015 15:44:24 +0200 Subject: [PATCH 4/7] Script for generating the powers-of-ten tables necessary for correct and fast decimal-to-float conversions. --- src/etc/dec2flt_table.py | 134 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 134 insertions(+) create mode 100644 src/etc/dec2flt_table.py diff --git a/src/etc/dec2flt_table.py b/src/etc/dec2flt_table.py new file mode 100644 index 0000000000000..b0140fb24559d --- /dev/null +++ b/src/etc/dec2flt_table.py @@ -0,0 +1,134 @@ +#!/usr/bin/env python2.7 +# +# 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 or the MIT license +# , at your +# option. This file may not be copied, modified, or distributed +# except according to those terms. + +""" +Generate powers of ten using William Clinger's ``AlgorithmM`` for use in +decimal to floating point conversions. + +Specifically, computes and outputs (as Rust code) a table of 10^e for some +range of exponents e. The output is one array of 64 bit significands and +another array of corresponding base two exponents. The approximations are +normalized and rounded perfectly, i.e., within 0.5 ULP of the true value. + +The representation ([u64], [i16]) instead of the more natural [(u64, i16)] +is used because (u64, i16) has a ton of padding which would make the table +even larger, and it's already uncomfortably large (6 KiB). +""" +from __future__ import print_function +import sys +from fractions import Fraction +from collections import namedtuple + + +N = 64 # Size of the significand field in bits +MIN_SIG = 2 ** (N - 1) +MAX_SIG = (2 ** N) - 1 + + +# Hand-rolled fp representation without arithmetic or any other operations. +# The significand is normalized and always N bit, but the exponent is +# unrestricted in range. +Fp = namedtuple('Fp', 'sig exp') + + +def algorithm_m(f, e): + assert f > 0 + if e < 0: + u = f + v = 10 ** abs(e) + else: + u = f * 10 ** e + v = 1 + k = 0 + x = u // v + while True: + if x < MIN_SIG: + u <<= 1 + k -= 1 + elif x >= MAX_SIG: + v <<= 1 + k += 1 + else: + break + x = u // v + return ratio_to_float(u, v, k) + + +def ratio_to_float(u, v, k): + q, r = divmod(u, v) + v_r = v - r + z = Fp(q, k) + if r < v_r: + return z + elif r > v_r: + return next_float(z) + elif q % 2 == 0: + return z + else: + return next_float(z) + + +def next_float(z): + if z.sig == MAX_SIG: + return Fp(MIN_SIG, z.exp + 1) + else: + return Fp(z.sig + 1, z.exp) + + +def error(f, e, z): + decimal = f * Fraction(10) ** e + binary = z.sig * Fraction(2) ** z.exp + abs_err = abs(decimal - binary) + # The unit in the last place has value z.exp + ulp_err = abs_err / Fraction(2) ** z.exp + return float(ulp_err) + +LICENSE = """ +// 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 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. +""" + +def main(): + MIN_E = -305 + MAX_E = 305 + e_range = range(MIN_E, MAX_E+1) + powers = [] + for e in e_range: + z = algorithm_m(1, e) + err = error(1, e, z) + assert err < 0.5 + powers.append(z) + typ = "([u64; {0}], [i16; {0}])".format(len(e_range)) + print(LICENSE.strip()) + print("// Table of approximations of powers of ten.") + print("// DO NOT MODIFY: Generated by a src/etc/dec2flt_table.py") + print("pub const MIN_E: i16 = {};".format(MIN_E)) + print("pub const MAX_E: i16 = {};".format(MAX_E)) + print() + print("pub const POWERS: ", typ, " = ([", sep='') + for z in powers: + print(" 0x{:x},".format(z.sig)) + print("], [") + for z in powers: + print(" {},".format(z.exp)) + print("]);") + + +if __name__ == '__main__': + main() From ba792a4baa856d83c3001afa181db91c5b4c9732 Mon Sep 17 00:00:00 2001 From: Robin Kruppe Date: Sun, 26 Jul 2015 17:50:29 +0200 Subject: [PATCH 5/7] Accurate decimal-to-float parsing routines. 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] --- src/libcore/num/dec2flt/algorithm.rs | 353 ++++++++ src/libcore/num/dec2flt/mod.rs | 234 +++++ src/libcore/num/dec2flt/num.rs | 95 ++ src/libcore/num/dec2flt/parse.rs | 128 +++ src/libcore/num/dec2flt/rawfp.rs | 356 ++++++++ src/libcore/num/dec2flt/table.rs | 1239 ++++++++++++++++++++++++++ src/libcore/num/flt2dec/bignum.rs | 21 +- src/libcore/num/mod.rs | 9 +- src/libcoretest/lib.rs | 1 + src/libcoretest/num/dec2flt/mod.rs | 174 ++++ src/libcoretest/num/dec2flt/parse.rs | 52 ++ src/libcoretest/num/dec2flt/rawfp.rs | 139 +++ src/libcoretest/num/mod.rs | 1 + 13 files changed, 2787 insertions(+), 15 deletions(-) create mode 100644 src/libcore/num/dec2flt/algorithm.rs create mode 100644 src/libcore/num/dec2flt/mod.rs create mode 100644 src/libcore/num/dec2flt/num.rs create mode 100644 src/libcore/num/dec2flt/parse.rs create mode 100644 src/libcore/num/dec2flt/rawfp.rs create mode 100644 src/libcore/num/dec2flt/table.rs create mode 100644 src/libcoretest/num/dec2flt/mod.rs create mode 100644 src/libcoretest/num/dec2flt/parse.rs create mode 100644 src/libcoretest/num/dec2flt/rawfp.rs diff --git a/src/libcore/num/dec2flt/algorithm.rs b/src/libcore/num/dec2flt/algorithm.rs new file mode 100644 index 0000000000000..97019090b56c7 --- /dev/null +++ b/src/libcore/num/dec2flt/algorithm.rs @@ -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 or the MIT license +// , 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(integral: &[u8], fractional: &[u8], e: i64) -> Option { + 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(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(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(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::(&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(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(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(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) + } +} diff --git a/src/libcore/num/dec2flt/mod.rs b/src/libcore/num/dec2flt/mod.rs new file mode 100644 index 0000000000000..413a67b256321 --- /dev/null +++ b/src/libcore/num/dec2flt/mod.rs @@ -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 or the MIT license +// , 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::. + +#![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 { + dec2flt(s) +} + +/// Entry point for decimal-to-f64 conversion. +pub fn to_f64(s: &str) -> Result { + 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(s: &str) -> Result { + 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(mut decimal: Decimal) -> Result { + 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(decimal: &Decimal) -> Option { + // 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 +} diff --git a/src/libcore/num/dec2flt/num.rs b/src/libcore/num/dec2flt/num.rs new file mode 100644 index 0000000000000..dcba73d7c9322 --- /dev/null +++ b/src/libcore/num/dec2flt/num.rs @@ -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 or the MIT license +// , 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 { + 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 +} diff --git a/src/libcore/num/dec2flt/parse.rs b/src/libcore/num/dec2flt/parse.rs new file mode 100644 index 0000000000000..58e2a6e9bba46 --- /dev/null +++ b/src/libcore/num/dec2flt/parse.rs @@ -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 or the MIT license +// , 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)) +} diff --git a/src/libcore/num/dec2flt/rawfp.rs b/src/libcore/num/dec2flt/rawfp.rs new file mode 100644 index 0000000000000..830d2dad42fe5 --- /dev/null +++ b/src/libcore/num/dec2flt/rawfp.rs @@ -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 or the MIT license +// , 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 + Div + Neg +{ + /// 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(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::(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(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(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(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(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(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) + } + } +} diff --git a/src/libcore/num/dec2flt/table.rs b/src/libcore/num/dec2flt/table.rs new file mode 100644 index 0000000000000..dd985fd155b85 --- /dev/null +++ b/src/libcore/num/dec2flt/table.rs @@ -0,0 +1,1239 @@ +// 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 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. +// Table of approximations of powers of ten. +// DO NOT MODIFY: Generated by a src/etc/dec2flt_table.py +pub const MIN_E: i16 = -305; +pub const MAX_E: i16 = 305; + +pub const POWERS: ([u64; 611], [i16; 611]) = ([ + 0xe0b62e2929aba83c, + 0x8c71dcd9ba0b4926, + 0xaf8e5410288e1b6f, + 0xdb71e91432b1a24b, + 0x892731ac9faf056f, + 0xab70fe17c79ac6ca, + 0xd64d3d9db981787d, + 0x85f0468293f0eb4e, + 0xa76c582338ed2622, + 0xd1476e2c07286faa, + 0x82cca4db847945ca, + 0xa37fce126597973d, + 0xcc5fc196fefd7d0c, + 0xff77b1fcbebcdc4f, + 0x9faacf3df73609b1, + 0xc795830d75038c1e, + 0xf97ae3d0d2446f25, + 0x9becce62836ac577, + 0xc2e801fb244576d5, + 0xf3a20279ed56d48a, + 0x9845418c345644d7, + 0xbe5691ef416bd60c, + 0xedec366b11c6cb8f, + 0x94b3a202eb1c3f39, + 0xb9e08a83a5e34f08, + 0xe858ad248f5c22ca, + 0x91376c36d99995be, + 0xb58547448ffffb2e, + 0xe2e69915b3fff9f9, + 0x8dd01fad907ffc3c, + 0xb1442798f49ffb4b, + 0xdd95317f31c7fa1d, + 0x8a7d3eef7f1cfc52, + 0xad1c8eab5ee43b67, + 0xd863b256369d4a41, + 0x873e4f75e2224e68, + 0xa90de3535aaae202, + 0xd3515c2831559a83, + 0x8412d9991ed58092, + 0xa5178fff668ae0b6, + 0xce5d73ff402d98e4, + 0x80fa687f881c7f8e, + 0xa139029f6a239f72, + 0xc987434744ac874f, + 0xfbe9141915d7a922, + 0x9d71ac8fada6c9b5, + 0xc4ce17b399107c23, + 0xf6019da07f549b2b, + 0x99c102844f94e0fb, + 0xc0314325637a193a, + 0xf03d93eebc589f88, + 0x96267c7535b763b5, + 0xbbb01b9283253ca3, + 0xea9c227723ee8bcb, + 0x92a1958a7675175f, + 0xb749faed14125d37, + 0xe51c79a85916f485, + 0x8f31cc0937ae58d3, + 0xb2fe3f0b8599ef08, + 0xdfbdcece67006ac9, + 0x8bd6a141006042be, + 0xaecc49914078536d, + 0xda7f5bf590966849, + 0x888f99797a5e012d, + 0xaab37fd7d8f58179, + 0xd5605fcdcf32e1d7, + 0x855c3be0a17fcd26, + 0xa6b34ad8c9dfc070, + 0xd0601d8efc57b08c, + 0x823c12795db6ce57, + 0xa2cb1717b52481ed, + 0xcb7ddcdda26da269, + 0xfe5d54150b090b03, + 0x9efa548d26e5a6e2, + 0xc6b8e9b0709f109a, + 0xf867241c8cc6d4c1, + 0x9b407691d7fc44f8, + 0xc21094364dfb5637, + 0xf294b943e17a2bc4, + 0x979cf3ca6cec5b5b, + 0xbd8430bd08277231, + 0xece53cec4a314ebe, + 0x940f4613ae5ed137, + 0xb913179899f68584, + 0xe757dd7ec07426e5, + 0x9096ea6f3848984f, + 0xb4bca50b065abe63, + 0xe1ebce4dc7f16dfc, + 0x8d3360f09cf6e4bd, + 0xb080392cc4349ded, + 0xdca04777f541c568, + 0x89e42caaf9491b61, + 0xac5d37d5b79b6239, + 0xd77485cb25823ac7, + 0x86a8d39ef77164bd, + 0xa8530886b54dbdec, + 0xd267caa862a12d67, + 0x8380dea93da4bc60, + 0xa46116538d0deb78, + 0xcd795be870516656, + 0x806bd9714632dff6, + 0xa086cfcd97bf97f4, + 0xc8a883c0fdaf7df0, + 0xfad2a4b13d1b5d6c, + 0x9cc3a6eec6311a64, + 0xc3f490aa77bd60fd, + 0xf4f1b4d515acb93c, + 0x991711052d8bf3c5, + 0xbf5cd54678eef0b7, + 0xef340a98172aace5, + 0x9580869f0e7aac0f, + 0xbae0a846d2195713, + 0xe998d258869facd7, + 0x91ff83775423cc06, + 0xb67f6455292cbf08, + 0xe41f3d6a7377eeca, + 0x8e938662882af53e, + 0xb23867fb2a35b28e, + 0xdec681f9f4c31f31, + 0x8b3c113c38f9f37f, + 0xae0b158b4738705f, + 0xd98ddaee19068c76, + 0x87f8a8d4cfa417ca, + 0xa9f6d30a038d1dbc, + 0xd47487cc8470652b, + 0x84c8d4dfd2c63f3b, + 0xa5fb0a17c777cf0a, + 0xcf79cc9db955c2cc, + 0x81ac1fe293d599c0, + 0xa21727db38cb0030, + 0xca9cf1d206fdc03c, + 0xfd442e4688bd304b, + 0x9e4a9cec15763e2f, + 0xc5dd44271ad3cdba, + 0xf7549530e188c129, + 0x9a94dd3e8cf578ba, + 0xc13a148e3032d6e8, + 0xf18899b1bc3f8ca2, + 0x96f5600f15a7b7e5, + 0xbcb2b812db11a5de, + 0xebdf661791d60f56, + 0x936b9fcebb25c996, + 0xb84687c269ef3bfb, + 0xe65829b3046b0afa, + 0x8ff71a0fe2c2e6dc, + 0xb3f4e093db73a093, + 0xe0f218b8d25088b8, + 0x8c974f7383725573, + 0xafbd2350644eead0, + 0xdbac6c247d62a584, + 0x894bc396ce5da772, + 0xab9eb47c81f5114f, + 0xd686619ba27255a3, + 0x8613fd0145877586, + 0xa798fc4196e952e7, + 0xd17f3b51fca3a7a1, + 0x82ef85133de648c5, + 0xa3ab66580d5fdaf6, + 0xcc963fee10b7d1b3, + 0xffbbcfe994e5c620, + 0x9fd561f1fd0f9bd4, + 0xc7caba6e7c5382c9, + 0xf9bd690a1b68637b, + 0x9c1661a651213e2d, + 0xc31bfa0fe5698db8, + 0xf3e2f893dec3f126, + 0x986ddb5c6b3a76b8, + 0xbe89523386091466, + 0xee2ba6c0678b597f, + 0x94db483840b717f0, + 0xba121a4650e4ddec, + 0xe896a0d7e51e1566, + 0x915e2486ef32cd60, + 0xb5b5ada8aaff80b8, + 0xe3231912d5bf60e6, + 0x8df5efabc5979c90, + 0xb1736b96b6fd83b4, + 0xddd0467c64bce4a1, + 0x8aa22c0dbef60ee4, + 0xad4ab7112eb3929e, + 0xd89d64d57a607745, + 0x87625f056c7c4a8b, + 0xa93af6c6c79b5d2e, + 0xd389b47879823479, + 0x843610cb4bf160cc, + 0xa54394fe1eedb8ff, + 0xce947a3da6a9273e, + 0x811ccc668829b887, + 0xa163ff802a3426a9, + 0xc9bcff6034c13053, + 0xfc2c3f3841f17c68, + 0x9d9ba7832936edc1, + 0xc5029163f384a931, + 0xf64335bcf065d37d, + 0x99ea0196163fa42e, + 0xc06481fb9bcf8d3a, + 0xf07da27a82c37088, + 0x964e858c91ba2655, + 0xbbe226efb628afeb, + 0xeadab0aba3b2dbe5, + 0x92c8ae6b464fc96f, + 0xb77ada0617e3bbcb, + 0xe55990879ddcaabe, + 0x8f57fa54c2a9eab7, + 0xb32df8e9f3546564, + 0xdff9772470297ebd, + 0x8bfbea76c619ef36, + 0xaefae51477a06b04, + 0xdab99e59958885c5, + 0x88b402f7fd75539b, + 0xaae103b5fcd2a882, + 0xd59944a37c0752a2, + 0x857fcae62d8493a5, + 0xa6dfbd9fb8e5b88f, + 0xd097ad07a71f26b2, + 0x825ecc24c8737830, + 0xa2f67f2dfa90563b, + 0xcbb41ef979346bca, + 0xfea126b7d78186bd, + 0x9f24b832e6b0f436, + 0xc6ede63fa05d3144, + 0xf8a95fcf88747d94, + 0x9b69dbe1b548ce7d, + 0xc24452da229b021c, + 0xf2d56790ab41c2a3, + 0x97c560ba6b0919a6, + 0xbdb6b8e905cb600f, + 0xed246723473e3813, + 0x9436c0760c86e30c, + 0xb94470938fa89bcf, + 0xe7958cb87392c2c3, + 0x90bd77f3483bb9ba, + 0xb4ecd5f01a4aa828, + 0xe2280b6c20dd5232, + 0x8d590723948a535f, + 0xb0af48ec79ace837, + 0xdcdb1b2798182245, + 0x8a08f0f8bf0f156b, + 0xac8b2d36eed2dac6, + 0xd7adf884aa879177, + 0x86ccbb52ea94baeb, + 0xa87fea27a539e9a5, + 0xd29fe4b18e88640f, + 0x83a3eeeef9153e89, + 0xa48ceaaab75a8e2b, + 0xcdb02555653131b6, + 0x808e17555f3ebf12, + 0xa0b19d2ab70e6ed6, + 0xc8de047564d20a8c, + 0xfb158592be068d2f, + 0x9ced737bb6c4183d, + 0xc428d05aa4751e4d, + 0xf53304714d9265e0, + 0x993fe2c6d07b7fac, + 0xbf8fdb78849a5f97, + 0xef73d256a5c0f77d, + 0x95a8637627989aae, + 0xbb127c53b17ec159, + 0xe9d71b689dde71b0, + 0x9226712162ab070e, + 0xb6b00d69bb55c8d1, + 0xe45c10c42a2b3b06, + 0x8eb98a7a9a5b04e3, + 0xb267ed1940f1c61c, + 0xdf01e85f912e37a3, + 0x8b61313bbabce2c6, + 0xae397d8aa96c1b78, + 0xd9c7dced53c72256, + 0x881cea14545c7575, + 0xaa242499697392d3, + 0xd4ad2dbfc3d07788, + 0x84ec3c97da624ab5, + 0xa6274bbdd0fadd62, + 0xcfb11ead453994ba, + 0x81ceb32c4b43fcf5, + 0xa2425ff75e14fc32, + 0xcad2f7f5359a3b3e, + 0xfd87b5f28300ca0e, + 0x9e74d1b791e07e48, + 0xc612062576589ddb, + 0xf79687aed3eec551, + 0x9abe14cd44753b53, + 0xc16d9a0095928a27, + 0xf1c90080baf72cb1, + 0x971da05074da7bef, + 0xbce5086492111aeb, + 0xec1e4a7db69561a5, + 0x9392ee8e921d5d07, + 0xb877aa3236a4b449, + 0xe69594bec44de15b, + 0x901d7cf73ab0acd9, + 0xb424dc35095cd80f, + 0xe12e13424bb40e13, + 0x8cbccc096f5088cc, + 0xafebff0bcb24aaff, + 0xdbe6fecebdedd5bf, + 0x89705f4136b4a597, + 0xabcc77118461cefd, + 0xd6bf94d5e57a42bc, + 0x8637bd05af6c69b6, + 0xa7c5ac471b478423, + 0xd1b71758e219652c, + 0x83126e978d4fdf3b, + 0xa3d70a3d70a3d70a, + 0xcccccccccccccccd, + 0x8000000000000000, + 0xa000000000000000, + 0xc800000000000000, + 0xfa00000000000000, + 0x9c40000000000000, + 0xc350000000000000, + 0xf424000000000000, + 0x9896800000000000, + 0xbebc200000000000, + 0xee6b280000000000, + 0x9502f90000000000, + 0xba43b74000000000, + 0xe8d4a51000000000, + 0x9184e72a00000000, + 0xb5e620f480000000, + 0xe35fa931a0000000, + 0x8e1bc9bf04000000, + 0xb1a2bc2ec5000000, + 0xde0b6b3a76400000, + 0x8ac7230489e80000, + 0xad78ebc5ac620000, + 0xd8d726b7177a8000, + 0x878678326eac9000, + 0xa968163f0a57b400, + 0xd3c21bcecceda100, + 0x84595161401484a0, + 0xa56fa5b99019a5c8, + 0xcecb8f27f4200f3a, + 0x813f3978f8940984, + 0xa18f07d736b90be5, + 0xc9f2c9cd04674edf, + 0xfc6f7c4045812296, + 0x9dc5ada82b70b59e, + 0xc5371912364ce305, + 0xf684df56c3e01bc7, + 0x9a130b963a6c115c, + 0xc097ce7bc90715b3, + 0xf0bdc21abb48db20, + 0x96769950b50d88f4, + 0xbc143fa4e250eb31, + 0xeb194f8e1ae525fd, + 0x92efd1b8d0cf37be, + 0xb7abc627050305ae, + 0xe596b7b0c643c719, + 0x8f7e32ce7bea5c70, + 0xb35dbf821ae4f38c, + 0xe0352f62a19e306f, + 0x8c213d9da502de45, + 0xaf298d050e4395d7, + 0xdaf3f04651d47b4c, + 0x88d8762bf324cd10, + 0xab0e93b6efee0054, + 0xd5d238a4abe98068, + 0x85a36366eb71f041, + 0xa70c3c40a64e6c52, + 0xd0cf4b50cfe20766, + 0x82818f1281ed44a0, + 0xa321f2d7226895c8, + 0xcbea6f8ceb02bb3a, + 0xfee50b7025c36a08, + 0x9f4f2726179a2245, + 0xc722f0ef9d80aad6, + 0xf8ebad2b84e0d58c, + 0x9b934c3b330c8577, + 0xc2781f49ffcfa6d5, + 0xf316271c7fc3908b, + 0x97edd871cfda3a57, + 0xbde94e8e43d0c8ec, + 0xed63a231d4c4fb27, + 0x945e455f24fb1cf9, + 0xb975d6b6ee39e437, + 0xe7d34c64a9c85d44, + 0x90e40fbeea1d3a4b, + 0xb51d13aea4a488dd, + 0xe264589a4dcdab15, + 0x8d7eb76070a08aed, + 0xb0de65388cc8ada8, + 0xdd15fe86affad912, + 0x8a2dbf142dfcc7ab, + 0xacb92ed9397bf996, + 0xd7e77a8f87daf7fc, + 0x86f0ac99b4e8dafd, + 0xa8acd7c0222311bd, + 0xd2d80db02aabd62c, + 0x83c7088e1aab65db, + 0xa4b8cab1a1563f52, + 0xcde6fd5e09abcf27, + 0x80b05e5ac60b6178, + 0xa0dc75f1778e39d6, + 0xc913936dd571c84c, + 0xfb5878494ace3a5f, + 0x9d174b2dcec0e47b, + 0xc45d1df942711d9a, + 0xf5746577930d6501, + 0x9968bf6abbe85f20, + 0xbfc2ef456ae276e9, + 0xefb3ab16c59b14a3, + 0x95d04aee3b80ece6, + 0xbb445da9ca61281f, + 0xea1575143cf97227, + 0x924d692ca61be758, + 0xb6e0c377cfa2e12e, + 0xe498f455c38b997a, + 0x8edf98b59a373fec, + 0xb2977ee300c50fe7, + 0xdf3d5e9bc0f653e1, + 0x8b865b215899f46d, + 0xae67f1e9aec07188, + 0xda01ee641a708dea, + 0x884134fe908658b2, + 0xaa51823e34a7eedf, + 0xd4e5e2cdc1d1ea96, + 0x850fadc09923329e, + 0xa6539930bf6bff46, + 0xcfe87f7cef46ff17, + 0x81f14fae158c5f6e, + 0xa26da3999aef774a, + 0xcb090c8001ab551c, + 0xfdcb4fa002162a63, + 0x9e9f11c4014dda7e, + 0xc646d63501a1511e, + 0xf7d88bc24209a565, + 0x9ae757596946075f, + 0xc1a12d2fc3978937, + 0xf209787bb47d6b85, + 0x9745eb4d50ce6333, + 0xbd176620a501fc00, + 0xec5d3fa8ce427b00, + 0x93ba47c980e98ce0, + 0xb8a8d9bbe123f018, + 0xe6d3102ad96cec1e, + 0x9043ea1ac7e41393, + 0xb454e4a179dd1877, + 0xe16a1dc9d8545e95, + 0x8ce2529e2734bb1d, + 0xb01ae745b101e9e4, + 0xdc21a1171d42645d, + 0x899504ae72497eba, + 0xabfa45da0edbde69, + 0xd6f8d7509292d603, + 0x865b86925b9bc5c2, + 0xa7f26836f282b733, + 0xd1ef0244af2364ff, + 0x8335616aed761f1f, + 0xa402b9c5a8d3a6e7, + 0xcd036837130890a1, + 0x802221226be55a65, + 0xa02aa96b06deb0fe, + 0xc83553c5c8965d3d, + 0xfa42a8b73abbf48d, + 0x9c69a97284b578d8, + 0xc38413cf25e2d70e, + 0xf46518c2ef5b8cd1, + 0x98bf2f79d5993803, + 0xbeeefb584aff8604, + 0xeeaaba2e5dbf6785, + 0x952ab45cfa97a0b3, + 0xba756174393d88e0, + 0xe912b9d1478ceb17, + 0x91abb422ccb812ef, + 0xb616a12b7fe617aa, + 0xe39c49765fdf9d95, + 0x8e41ade9fbebc27d, + 0xb1d219647ae6b31c, + 0xde469fbd99a05fe3, + 0x8aec23d680043bee, + 0xada72ccc20054aea, + 0xd910f7ff28069da4, + 0x87aa9aff79042287, + 0xa99541bf57452b28, + 0xd3fa922f2d1675f2, + 0x847c9b5d7c2e09b7, + 0xa59bc234db398c25, + 0xcf02b2c21207ef2f, + 0x8161afb94b44f57d, + 0xa1ba1ba79e1632dc, + 0xca28a291859bbf93, + 0xfcb2cb35e702af78, + 0x9defbf01b061adab, + 0xc56baec21c7a1916, + 0xf6c69a72a3989f5c, + 0x9a3c2087a63f6399, + 0xc0cb28a98fcf3c80, + 0xf0fdf2d3f3c30b9f, + 0x969eb7c47859e744, + 0xbc4665b596706115, + 0xeb57ff22fc0c795a, + 0x9316ff75dd87cbd8, + 0xb7dcbf5354e9bece, + 0xe5d3ef282a242e82, + 0x8fa475791a569d11, + 0xb38d92d760ec4455, + 0xe070f78d3927556b, + 0x8c469ab843b89563, + 0xaf58416654a6babb, + 0xdb2e51bfe9d0696a, + 0x88fcf317f22241e2, + 0xab3c2fddeeaad25b, + 0xd60b3bd56a5586f2, + 0x85c7056562757457, + 0xa738c6bebb12d16d, + 0xd106f86e69d785c8, + 0x82a45b450226b39d, + 0xa34d721642b06084, + 0xcc20ce9bd35c78a5, + 0xff290242c83396ce, + 0x9f79a169bd203e41, + 0xc75809c42c684dd1, + 0xf92e0c3537826146, + 0x9bbcc7a142b17ccc, + 0xc2abf989935ddbfe, + 0xf356f7ebf83552fe, + 0x98165af37b2153df, + 0xbe1bf1b059e9a8d6, + 0xeda2ee1c7064130c, + 0x9485d4d1c63e8be8, + 0xb9a74a0637ce2ee1, + 0xe8111c87c5c1ba9a, + 0x910ab1d4db9914a0, + 0xb54d5e4a127f59c8, + 0xe2a0b5dc971f303a, + 0x8da471a9de737e24, + 0xb10d8e1456105dad, + 0xdd50f1996b947519, + 0x8a5296ffe33cc930, + 0xace73cbfdc0bfb7b, + 0xd8210befd30efa5a, + 0x8714a775e3e95c78, + 0xa8d9d1535ce3b396, + 0xd31045a8341ca07c, + 0x83ea2b892091e44e, + 0xa4e4b66b68b65d61, + 0xce1de40642e3f4b9, + 0x80d2ae83e9ce78f4, + 0xa1075a24e4421731, + 0xc94930ae1d529cfd, + 0xfb9b7cd9a4a7443c, + 0x9d412e0806e88aa6, + 0xc491798a08a2ad4f, + 0xf5b5d7ec8acb58a3, + 0x9991a6f3d6bf1766, + 0xbff610b0cc6edd3f, + 0xeff394dcff8a948f, + 0x95f83d0a1fb69cd9, + 0xbb764c4ca7a44410, + 0xea53df5fd18d5514, + 0x92746b9be2f8552c, + 0xb7118682dbb66a77, + 0xe4d5e82392a40515, + 0x8f05b1163ba6832d, + 0xb2c71d5bca9023f8, + 0xdf78e4b2bd342cf7, + 0x8bab8eefb6409c1a, + 0xae9672aba3d0c321, + 0xda3c0f568cc4f3e9, + 0x8865899617fb1871, + 0xaa7eebfb9df9de8e, + 0xd51ea6fa85785631, + 0x8533285c936b35df, + 0xa67ff273b8460357, + 0xd01fef10a657842c, + 0x8213f56a67f6b29c, + 0xa298f2c501f45f43, + 0xcb3f2f7642717713, + 0xfe0efb53d30dd4d8, + 0x9ec95d1463e8a507, + 0xc67bb4597ce2ce49, + 0xf81aa16fdc1b81db, + 0x9b10a4e5e9913129, + 0xc1d4ce1f63f57d73, + 0xf24a01a73cf2dcd0, + 0x976e41088617ca02, + 0xbd49d14aa79dbc82, + 0xec9c459d51852ba3, + 0x93e1ab8252f33b46, + 0xb8da1662e7b00a17, + 0xe7109bfba19c0c9d, + 0x906a617d450187e2, + 0xb484f9dc9641e9db, + 0xe1a63853bbd26451, + 0x8d07e33455637eb3, + 0xb049dc016abc5e60, + 0xdc5c5301c56b75f7, + 0x89b9b3e11b6329bb, + 0xac2820d9623bf429, + 0xd732290fbacaf134, + 0x867f59a9d4bed6c0, + 0xa81f301449ee8c70, + 0xd226fc195c6a2f8c, + 0x83585d8fd9c25db8, + 0xa42e74f3d032f526, + 0xcd3a1230c43fb26f, + 0x80444b5e7aa7cf85, + 0xa0555e361951c367, + 0xc86ab5c39fa63441, + 0xfa856334878fc151, + 0x9c935e00d4b9d8d2, + 0xc3b8358109e84f07, + 0xf4a642e14c6262c9, + 0x98e7e9cccfbd7dbe, + 0xbf21e44003acdd2d, + 0xeeea5d5004981478, + 0x95527a5202df0ccb, + 0xbaa718e68396cffe, + 0xe950df20247c83fd, + 0x91d28b7416cdd27e, +], [ + -1077, + -1073, + -1070, + -1067, + -1063, + -1060, + -1057, + -1053, + -1050, + -1047, + -1043, + -1040, + -1037, + -1034, + -1030, + -1027, + -1024, + -1020, + -1017, + -1014, + -1010, + -1007, + -1004, + -1000, + -997, + -994, + -990, + -987, + -984, + -980, + -977, + -974, + -970, + -967, + -964, + -960, + -957, + -954, + -950, + -947, + -944, + -940, + -937, + -934, + -931, + -927, + -924, + -921, + -917, + -914, + -911, + -907, + -904, + -901, + -897, + -894, + -891, + -887, + -884, + -881, + -877, + -874, + -871, + -867, + -864, + -861, + -857, + -854, + -851, + -847, + -844, + -841, + -838, + -834, + -831, + -828, + -824, + -821, + -818, + -814, + -811, + -808, + -804, + -801, + -798, + -794, + -791, + -788, + -784, + -781, + -778, + -774, + -771, + -768, + -764, + -761, + -758, + -754, + -751, + -748, + -744, + -741, + -738, + -735, + -731, + -728, + -725, + -721, + -718, + -715, + -711, + -708, + -705, + -701, + -698, + -695, + -691, + -688, + -685, + -681, + -678, + -675, + -671, + -668, + -665, + -661, + -658, + -655, + -651, + -648, + -645, + -642, + -638, + -635, + -632, + -628, + -625, + -622, + -618, + -615, + -612, + -608, + -605, + -602, + -598, + -595, + -592, + -588, + -585, + -582, + -578, + -575, + -572, + -568, + -565, + -562, + -558, + -555, + -552, + -549, + -545, + -542, + -539, + -535, + -532, + -529, + -525, + -522, + -519, + -515, + -512, + -509, + -505, + -502, + -499, + -495, + -492, + -489, + -485, + -482, + -479, + -475, + -472, + -469, + -465, + -462, + -459, + -455, + -452, + -449, + -446, + -442, + -439, + -436, + -432, + -429, + -426, + -422, + -419, + -416, + -412, + -409, + -406, + -402, + -399, + -396, + -392, + -389, + -386, + -382, + -379, + -376, + -372, + -369, + -366, + -362, + -359, + -356, + -353, + -349, + -346, + -343, + -339, + -336, + -333, + -329, + -326, + -323, + -319, + -316, + -313, + -309, + -306, + -303, + -299, + -296, + -293, + -289, + -286, + -283, + -279, + -276, + -273, + -269, + -266, + -263, + -259, + -256, + -253, + -250, + -246, + -243, + -240, + -236, + -233, + -230, + -226, + -223, + -220, + -216, + -213, + -210, + -206, + -203, + -200, + -196, + -193, + -190, + -186, + -183, + -180, + -176, + -173, + -170, + -166, + -163, + -160, + -157, + -153, + -150, + -147, + -143, + -140, + -137, + -133, + -130, + -127, + -123, + -120, + -117, + -113, + -110, + -107, + -103, + -100, + -97, + -93, + -90, + -87, + -83, + -80, + -77, + -73, + -70, + -67, + -63, + -60, + -57, + -54, + -50, + -47, + -44, + -40, + -37, + -34, + -30, + -27, + -24, + -20, + -17, + -14, + -10, + -7, + -4, + 0, + 3, + 6, + 10, + 13, + 16, + 20, + 23, + 26, + 30, + 33, + 36, + 39, + 43, + 46, + 49, + 53, + 56, + 59, + 63, + 66, + 69, + 73, + 76, + 79, + 83, + 86, + 89, + 93, + 96, + 99, + 103, + 106, + 109, + 113, + 116, + 119, + 123, + 126, + 129, + 132, + 136, + 139, + 142, + 146, + 149, + 152, + 156, + 159, + 162, + 166, + 169, + 172, + 176, + 179, + 182, + 186, + 189, + 192, + 196, + 199, + 202, + 206, + 209, + 212, + 216, + 219, + 222, + 226, + 229, + 232, + 235, + 239, + 242, + 245, + 249, + 252, + 255, + 259, + 262, + 265, + 269, + 272, + 275, + 279, + 282, + 285, + 289, + 292, + 295, + 299, + 302, + 305, + 309, + 312, + 315, + 319, + 322, + 325, + 328, + 332, + 335, + 338, + 342, + 345, + 348, + 352, + 355, + 358, + 362, + 365, + 368, + 372, + 375, + 378, + 382, + 385, + 388, + 392, + 395, + 398, + 402, + 405, + 408, + 412, + 415, + 418, + 422, + 425, + 428, + 431, + 435, + 438, + 441, + 445, + 448, + 451, + 455, + 458, + 461, + 465, + 468, + 471, + 475, + 478, + 481, + 485, + 488, + 491, + 495, + 498, + 501, + 505, + 508, + 511, + 515, + 518, + 521, + 524, + 528, + 531, + 534, + 538, + 541, + 544, + 548, + 551, + 554, + 558, + 561, + 564, + 568, + 571, + 574, + 578, + 581, + 584, + 588, + 591, + 594, + 598, + 601, + 604, + 608, + 611, + 614, + 617, + 621, + 624, + 627, + 631, + 634, + 637, + 641, + 644, + 647, + 651, + 654, + 657, + 661, + 664, + 667, + 671, + 674, + 677, + 681, + 684, + 687, + 691, + 694, + 697, + 701, + 704, + 707, + 711, + 714, + 717, + 720, + 724, + 727, + 730, + 734, + 737, + 740, + 744, + 747, + 750, + 754, + 757, + 760, + 764, + 767, + 770, + 774, + 777, + 780, + 784, + 787, + 790, + 794, + 797, + 800, + 804, + 807, + 810, + 813, + 817, + 820, + 823, + 827, + 830, + 833, + 837, + 840, + 843, + 847, + 850, + 853, + 857, + 860, + 863, + 867, + 870, + 873, + 877, + 880, + 883, + 887, + 890, + 893, + 897, + 900, + 903, + 907, + 910, + 913, + 916, + 920, + 923, + 926, + 930, + 933, + 936, + 940, + 943, + 946, + 950, +]); diff --git a/src/libcore/num/flt2dec/bignum.rs b/src/libcore/num/flt2dec/bignum.rs index 15980f925fa03..ee1f6ffdd0aef 100644 --- a/src/libcore/num/flt2dec/bignum.rs +++ b/src/libcore/num/flt2dec/bignum.rs @@ -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; diff --git a/src/libcore/num/mod.rs b/src/libcore/num/mod.rs index ad891bf8fa623..19ca83b6ea40a 100644 --- a/src/libcore/num/mod.rs +++ b/src/libcore/num/mod.rs @@ -44,6 +44,7 @@ pub struct Wrapping(#[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::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)*) => {$( diff --git a/src/libcoretest/lib.rs b/src/libcoretest/lib.rs index 08536e6320475..4262be83361a3 100644 --- a/src/libcoretest/lib.rs +++ b/src/libcoretest/lib.rs @@ -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)] diff --git a/src/libcoretest/num/dec2flt/mod.rs b/src/libcoretest/num/dec2flt/mod.rs new file mode 100644 index 0000000000000..bd8cfc74f0c63 --- /dev/null +++ b/src/libcoretest/num/dec2flt/mod.rs @@ -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 or the MIT license +// , 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")); +} diff --git a/src/libcoretest/num/dec2flt/parse.rs b/src/libcoretest/num/dec2flt/parse.rs new file mode 100644 index 0000000000000..09acf2bc517b0 --- /dev/null +++ b/src/libcoretest/num/dec2flt/parse.rs @@ -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 or the MIT license +// , 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,? or the MIT license +// , 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::(Fp { f: sig, e: 0 }).integer_decode(); + assert_eq!(e1, 0 + 64 - 53); + let (m2, e2, _) = fp_to_float::(Fp { f: sig, e: 55 }).integer_decode(); + assert_eq!(e2, 55 + 64 - 53); + assert_eq!(m2, m1); + let (m3, e3, _) = fp_to_float::(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::(Fp { f: 1, e: 0 }), 1.0); + assert_eq!(fp_to_float::(Fp { f: 42, e: 7 }), (42 << 7) as f64); + assert_eq!(fp_to_float::(Fp { f: 1 << 20, e: 30 }), (1u64 << 50) as f64); + assert_eq!(fp_to_float::(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::(fp), x); + } +} + +#[test] +fn rounding_overflow() { + let x = Fp { f: 0xFF_FF_FF_FF_FF_FF_FF_00u64, e: 42 }; + let rounded = round_normal::(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); +} diff --git a/src/libcoretest/num/mod.rs b/src/libcoretest/num/mod.rs index a9baa2cc477f6..9f9d2a4ca1659 100644 --- a/src/libcoretest/num/mod.rs +++ b/src/libcoretest/num/mod.rs @@ -30,6 +30,7 @@ mod u32; mod u64; mod flt2dec; +mod dec2flt; /// Helper function for testing numeric operations pub fn test_num(ten: T, two: T) where From 82dbc2ea619cbfc98ca9ad2f9e06a5acd294cbe3 Mon Sep 17 00:00:00 2001 From: Robin Kruppe Date: Sun, 26 Jul 2015 15:24:08 +0200 Subject: [PATCH 6/7] Add optional, external tests for floating point parsing. Running these tests takes hours, so they are not run by @bors. --- src/etc/test-float-parse/_common.rs | 26 ++ src/etc/test-float-parse/few-ones.rs | 27 ++ src/etc/test-float-parse/huge-pow10.rs | 21 ++ src/etc/test-float-parse/many-digits.rs | 39 ++ src/etc/test-float-parse/rand-f64.rs | 32 ++ src/etc/test-float-parse/runtests.py | 399 +++++++++++++++++++++ src/etc/test-float-parse/short-decimals.rs | 29 ++ src/etc/test-float-parse/subnorm.rs | 23 ++ src/etc/test-float-parse/tiny-pow10.rs | 21 ++ src/etc/test-float-parse/u32-small.rs | 19 + src/etc/test-float-parse/u64-pow2.rs | 28 ++ 11 files changed, 664 insertions(+) create mode 100644 src/etc/test-float-parse/_common.rs create mode 100644 src/etc/test-float-parse/few-ones.rs create mode 100644 src/etc/test-float-parse/huge-pow10.rs create mode 100644 src/etc/test-float-parse/many-digits.rs create mode 100644 src/etc/test-float-parse/rand-f64.rs create mode 100644 src/etc/test-float-parse/runtests.py create mode 100644 src/etc/test-float-parse/short-decimals.rs create mode 100644 src/etc/test-float-parse/subnorm.rs create mode 100644 src/etc/test-float-parse/tiny-pow10.rs create mode 100644 src/etc/test-float-parse/u32-small.rs create mode 100644 src/etc/test-float-parse/u64-pow2.rs diff --git a/src/etc/test-float-parse/_common.rs b/src/etc/test-float-parse/_common.rs new file mode 100644 index 0000000000000..b4a2a593e0118 --- /dev/null +++ b/src/etc/test-float-parse/_common.rs @@ -0,0 +1,26 @@ +// 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 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. + +use std::io; +use std::io::prelude::*; +use std::mem::transmute; + +// Nothing up my sleeve: Just (PI - 3) in base 16. +#[allow(dead_code)] +pub const SEED: [u32; 3] = [0x243f_6a88, 0x85a3_08d3, 0x1319_8a2e]; + +pub fn validate(text: String) { + let mut out = io::stdout(); + let x: f64 = text.parse().unwrap(); + let f64_bytes: u64 = unsafe { transmute(x) }; + let x: f32 = text.parse().unwrap(); + let f32_bytes: u32 = unsafe { transmute(x) }; + writeln!(&mut out, "{:016x} {:08x} {}", f64_bytes, f32_bytes, text).unwrap(); +} diff --git a/src/etc/test-float-parse/few-ones.rs b/src/etc/test-float-parse/few-ones.rs new file mode 100644 index 0000000000000..390f4da6f63ae --- /dev/null +++ b/src/etc/test-float-parse/few-ones.rs @@ -0,0 +1,27 @@ +// 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 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. + +mod _common; + +use _common::validate; + +fn main() { + let mut pow = vec![]; + for i in 0..63 { + pow.push(1u64 << i); + } + for a in &pow { + for b in &pow { + for c in &pow { + validate((a | b | c).to_string()); + } + } + } +} diff --git a/src/etc/test-float-parse/huge-pow10.rs b/src/etc/test-float-parse/huge-pow10.rs new file mode 100644 index 0000000000000..e62afc7851564 --- /dev/null +++ b/src/etc/test-float-parse/huge-pow10.rs @@ -0,0 +1,21 @@ +// 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 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. + +mod _common; + +use _common::validate; + +fn main() { + for e in 300..310 { + for i in 0..100000 { + validate(format!("{}e{}", i, e)); + } + } +} diff --git a/src/etc/test-float-parse/many-digits.rs b/src/etc/test-float-parse/many-digits.rs new file mode 100644 index 0000000000000..0cbf57183dfc2 --- /dev/null +++ b/src/etc/test-float-parse/many-digits.rs @@ -0,0 +1,39 @@ +// 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 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. + +#![feature(rand)] + +extern crate rand; + +mod _common; + +use std::char; +use rand::{IsaacRng, Rng, SeedableRng}; +use rand::distributions::{Range, Sample}; +use _common::{validate, SEED}; + +fn main() { + let mut rnd = IsaacRng::from_seed(&SEED); + let mut range = Range::new(0, 10); + for _ in 0..5_000_000u64 { + let num_digits = rnd.gen_range(100, 300); + let digits = gen_digits(num_digits, &mut range, &mut rnd); + validate(digits); + } +} + +fn gen_digits(n: u32, range: &mut Range, rnd: &mut R) -> String { + let mut s = String::new(); + for _ in 0..n { + let digit = char::from_digit(range.sample(rnd), 10).unwrap(); + s.push(digit); + } + s +} diff --git a/src/etc/test-float-parse/rand-f64.rs b/src/etc/test-float-parse/rand-f64.rs new file mode 100644 index 0000000000000..762c3d95ec6eb --- /dev/null +++ b/src/etc/test-float-parse/rand-f64.rs @@ -0,0 +1,32 @@ +// 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 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. + +#![feature(rand)] + +extern crate rand; + +mod _common; + +use _common::{validate, SEED}; +use rand::{IsaacRng, Rng, SeedableRng}; +use std::mem::transmute; + +fn main() { + let mut rnd = IsaacRng::from_seed(&SEED); + let mut i = 0; + while i < 10_000_000 { + let bits = rnd.next_u64(); + let x: f64 = unsafe { transmute(bits) }; + if x.is_finite() { + validate(format!("{:e}", x)); + i += 1; + } + } +} diff --git a/src/etc/test-float-parse/runtests.py b/src/etc/test-float-parse/runtests.py new file mode 100644 index 0000000000000..17a1b769bd63f --- /dev/null +++ b/src/etc/test-float-parse/runtests.py @@ -0,0 +1,399 @@ +#!/usr/bin/python2.7 +# +# 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 or the MIT license +# , at your +# option. This file may not be copied, modified, or distributed +# except according to those terms. + +""" +Testing dec2flt +=============== +These are *really* extensive tests. Expect them to run for hours. Due to the +nature of the problem (the input is a string of arbitrary length), exhaustive +testing is not really possible. Instead, there are exhaustive tests for some +classes of inputs for which that is feasible and a bunch of deterministic and +random non-exhaustive tests for covering everything else. + +The actual tests (generating decimal strings and feeding them to dec2flt) is +performed by a set of stand-along rust programs. This script compiles, runs, +and supervises them. In particular, the programs report the strings they +generate and the floating point numbers they converted those strings to. + +You can run specific tests rather than all of them by giving their names +(without .rs extension) as command line parameters. + +Verification +------------ +The tricky part is not generating those inputs but verifying the outputs. +Comparing with the result of Python's float() does not cut it because +(and this is apparently undocumented) although Python includes a version of +Martin Gay's code including the decimal-to-float part, it doesn't actually use +it for float() (only for round()) instead relying on the system scanf() which +is not necessarily completely accurate. + +Instead, we take the input and compute the true value with bignum arithmetic +(as a fraction, using the ``fractions`` module). + +Given an input string and the corresponding float computed via Rust, simply +decode the float into f * 2^k (for intergers f, k) and the ULP. +We can now easily compute the error and check if it is within 0.5 ULP as it +should be. Zero and infinites are handled similarly: + +- If the approximation is 0.0, the exact value should be *less or equal* + half the smallest denormal float: the smallest denormal floating point + number has an odd mantissa (00...001) and thus half of that is rounded + to 00...00, i.e., zero. +- If the approximation is Inf, the exact value should be *greater or equal* + to the largest finite float + 0.5 ULP: the largest finite float has an odd + mantissa (11...11), so that plus half an ULP is rounded up to the nearest + even number, which overflows. + +Implementation details +---------------------- +This directory contains a set of single-file Rust programs that perform +tests with a particular class of inputs. Each is compiled and run without +parameters, outputs (f64, f32, decimal) pairs to verify externally, and +in any case either exits gracefully or with a panic. + +If a test binary writes *anything at all* to stderr or exits with an +exit code that's not 0, the test fails. +The output on stdout is treated as (f64, f32, decimal) record, encoded thusly: + +- The first eight bytes are a binary64 (native endianness). +- The following four bytes are a binary32 (native endianness). +- Then the corresponding string input follows, in ASCII (no newline). +- The record is terminated with a newline. + +Incomplete records are an error. Not-a-Number bit patterns are invalid too. + +The tests run serially but the validaition for a a single test is parallelized +with ``multiprocessing``. Each test is launched as a subprocess. +One thread supervises it: Accepts and enqueues records to validate, observe +stderr, and waits for the process to exit. A set of worker processes perform +the validation work for the outputs enqueued there. Another thread listens +for progress updates from the workers. + +Known issues +------------ +Some errors (e.g., NaN outputs) aren't handled very gracefully. +Also, if there is an exception or the process is interrupted (at least on +Windows) the worker processes are leaked and stick around forever. +They're only a few megabytes each, but still, this script should not be run +if you aren't prepared to manually kill a lot of orphaned processes. +""" +from __future__ import print_function +import sys +import os.path +import time +import struct +from fractions import Fraction +from collections import namedtuple +from subprocess import Popen, check_call, PIPE +from glob import glob +import multiprocessing +import Queue +import threading +import ctypes +import binascii + +NUM_WORKERS = 2 +UPDATE_EVERY_N = 50000 +INF = namedtuple('INF', '')() +NEG_INF = namedtuple('NEG_INF', '')() +ZERO = namedtuple('ZERO', '')() +MAILBOX = None # The queue for reporting errors to the main process. +STDOUT_LOCK = threading.Lock() +test_name = None +child_processes = [] +exit_status = 0 + +def msg(*args): + with STDOUT_LOCK: + print("[" + test_name + "]", *args) + sys.stdout.flush() + + +def write_errors(): + global exit_status + f = open("errors.txt", 'w') + have_seen_error = False + while True: + args = MAILBOX.get() + if args is None: + f.close() + break + print(*args, file=f) + f.flush() + if not have_seen_error: + have_seen_error = True + msg("Something is broken:", *args) + msg("Future errors logged to errors.txt") + exit_status = 101 + + +def rustc(test): + rs = test + '.rs' + exe = test + '.exe' # hopefully this makes it work on *nix + print("compiling", test) + sys.stdout.flush() + check_call(['rustc', rs, '-o', exe]) + + +def run(test): + global test_name + test_name = test + + t0 = time.clock() + msg("setting up supervisor") + exe = test + '.exe' + proc = Popen(exe, bufsize=1<<20 , stdin=PIPE, stdout=PIPE, stderr=PIPE) + done = multiprocessing.Value(ctypes.c_bool) + queue = multiprocessing.Queue(maxsize=5)#(maxsize=1024) + workers = [] + for n in range(NUM_WORKERS): + worker = multiprocessing.Process(name='Worker-' + str(n + 1), + target=init_worker, + args=[test, MAILBOX, queue, done]) + workers.append(worker) + child_processes.append(worker) + for worker in workers: + worker.start() + msg("running test") + interact(proc, queue) + with done.get_lock(): + done.value = True + for worker in workers: + worker.join() + msg("python is done") + assert queue.empty(), "did not validate everything" + dt = time.clock() - t0 + msg("took", round(dt, 3), "seconds") + + +def interact(proc, queue): + line = "" + n = 0 + while proc.poll() is None: + line = proc.stdout.readline() + if not line: + continue + assert line.endswith('\n'), "incomplete line: " + repr(line) + queue.put(line) + line = "" + n += 1 + if n % UPDATE_EVERY_N == 0: + msg("got", str(n // 1000) + "k", "records") + msg("rust is done. exit code:", proc.returncode) + rest, stderr = proc.communicate() + if stderr: + msg("rust stderr output:", stderr) + for line in rest.split('\n'): + if not line: + continue + queue.put(line) + + +def main(): + global MAILBOX + tests = [os.path.splitext(f)[0] for f in glob('*.rs') + if not f.startswith('_')] + whitelist = sys.argv[1:] + if whitelist: + tests = [test for test in tests if test in whitelist] + if not tests: + print("Error: No tests to run") + sys.exit(1) + # Compile first for quicker feedback + for test in tests: + rustc(test) + # Set up mailbox once for all tests + MAILBOX = multiprocessing.Queue() + mailman = threading.Thread(target=write_errors) + mailman.daemon = True + mailman.start() + for test in tests: + if whitelist and test not in whitelist: + continue + run(test) + MAILBOX.put(None) + mailman.join() + + +# ---- Worker thread code ---- + + +POW2 = { e: Fraction(2) ** e for e in range(-1100, 1100) } +HALF_ULP = { e: (Fraction(2) ** e)/2 for e in range(-1100, 1100) } +DONE_FLAG = None + + +def send_error_to_supervisor(*args): + MAILBOX.put(args) + + +def init_worker(test, mailbox, queue, done): + global test_name, MAILBOX, DONE_FLAG + test_name = test + MAILBOX = mailbox + DONE_FLAG = done + do_work(queue) + + +def is_done(): + with DONE_FLAG.get_lock(): + return DONE_FLAG.value + + +def do_work(queue): + while True: + try: + line = queue.get(timeout=0.01) + except Queue.Empty: + if queue.empty() and is_done(): + return + else: + continue + bin64, bin32, text = line.rstrip().split() + validate(bin64, bin32, text) + + +def decode_binary64(x): + """ + Turn a IEEE 754 binary64 into (mantissa, exponent), except 0.0 and + infinity (positive and negative), which return ZERO, INF, and NEG_INF + respectively. + """ + x = binascii.unhexlify(x) + assert len(x) == 8, repr(x) + [bits] = struct.unpack(b'>Q', x) + if bits == 0: + return ZERO + exponent = (bits >> 52) & 0x7FF + negative = bits >> 63 + low_bits = bits & 0xFFFFFFFFFFFFF + if exponent == 0: + mantissa = low_bits + exponent += 1 + if mantissa == 0: + return ZERO + elif exponent == 0x7FF: + assert low_bits == 0, "NaN" + if negative: + return NEG_INF + else: + return INF + else: + mantissa = low_bits | (1 << 52) + exponent -= 1023 + 52 + if negative: + mantissa = -mantissa + return (mantissa, exponent) + + +def decode_binary32(x): + """ + Turn a IEEE 754 binary32 into (mantissa, exponent), except 0.0 and + infinity (positive and negative), which return ZERO, INF, and NEG_INF + respectively. + """ + x = binascii.unhexlify(x) + assert len(x) == 4, repr(x) + [bits] = struct.unpack(b'>I', x) + if bits == 0: + return ZERO + exponent = (bits >> 23) & 0xFF + negative = bits >> 31 + low_bits = bits & 0x7FFFFF + if exponent == 0: + mantissa = low_bits + exponent += 1 + if mantissa == 0: + return ZERO + elif exponent == 0xFF: + if negative: + return NEG_INF + else: + return INF + else: + mantissa = low_bits | (1 << 23) + exponent -= 127 + 23 + if negative: + mantissa = -mantissa + return (mantissa, exponent) + + +MIN_SUBNORMAL_DOUBLE = Fraction(2) ** -1074 +MIN_SUBNORMAL_SINGLE = Fraction(2) ** -149 # XXX unsure +MAX_DOUBLE = (2 - Fraction(2) ** -52) * (2 ** 1023) +MAX_SINGLE = (2 - Fraction(2) ** -23) * (2 ** 127) +MAX_ULP_DOUBLE = 1023 - 52 +MAX_ULP_SINGLE = 127 - 23 +DOUBLE_ZERO_CUTOFF = MIN_SUBNORMAL_DOUBLE / 2 +DOUBLE_INF_CUTOFF = MAX_DOUBLE + 2 ** (MAX_ULP_DOUBLE - 1) +SINGLE_ZERO_CUTOFF = MIN_SUBNORMAL_SINGLE / 2 +SINGLE_INF_CUTOFF = MAX_SINGLE + 2 ** (MAX_ULP_SINGLE - 1) + +def validate(bin64, bin32, text): + double = decode_binary64(bin64) + single = decode_binary32(bin32) + real = Fraction(text) + + if double is ZERO: + if real > DOUBLE_ZERO_CUTOFF: + record_special_error(text, "f64 zero") + elif double is INF: + if real < DOUBLE_INF_CUTOFF: + record_special_error(text, "f64 inf") + elif double is NEG_INF: + if -real < DOUBLE_INF_CUTOFF: + record_special_error(text, "f64 -inf") + elif len(double) == 2: + sig, k = double + validate_normal(text, real, sig, k, "f64") + else: + assert 0, "didn't handle binary64" + if single is ZERO: + if real > SINGLE_ZERO_CUTOFF: + record_special_error(text, "f32 zero") + elif single is INF: + if real < SINGLE_INF_CUTOFF: + record_special_error(text, "f32 inf") + elif single is NEG_INF: + if -real < SINGLE_INF_CUTOFF: + record_special_error(text, "f32 -inf") + elif len(single) == 2: + sig, k = single + validate_normal(text, real, sig, k, "f32") + else: + assert 0, "didn't handle binary32" + +def record_special_error(text, descr): + send_error_to_supervisor(text.strip(), "wrongly rounded to", descr) + + +def validate_normal(text, real, sig, k, kind): + approx = sig * POW2[k] + error = abs(approx - real) + if error > HALF_ULP[k]: + record_normal_error(text, error, k, kind) + + +def record_normal_error(text, error, k, kind): + one_ulp = HALF_ULP[k + 1] + assert one_ulp == 2 * HALF_ULP[k] + relative_error = error / one_ulp + text = text.strip() + try: + err_repr = float(relative_error) + except ValueError: + err_repr = str(err_repr).replace('/', ' / ') + send_error_to_supervisor(err_repr, "ULP error on", text, "(" + kind + ")") + + +if __name__ == '__main__': + main() diff --git a/src/etc/test-float-parse/short-decimals.rs b/src/etc/test-float-parse/short-decimals.rs new file mode 100644 index 0000000000000..baefb9c930543 --- /dev/null +++ b/src/etc/test-float-parse/short-decimals.rs @@ -0,0 +1,29 @@ +// 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 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. + +mod _common; + +use _common::validate; + +fn main() { + // Skip e = 0 because small-u32 already does those. + for e in 1..301 { + for i in 0..10000 { + // If it ends in zeros, the parser will strip those (and adjust the exponent), + // which almost always (except for exponents near +/- 300) result in an input + // equivalent to something we already generate in a different way. + if i % 10 == 0 { + continue; + } + validate(format!("{}e{}", i, e)); + validate(format!("{}e-{}", i, e)); + } + } +} diff --git a/src/etc/test-float-parse/subnorm.rs b/src/etc/test-float-parse/subnorm.rs new file mode 100644 index 0000000000000..70682c9b21810 --- /dev/null +++ b/src/etc/test-float-parse/subnorm.rs @@ -0,0 +1,23 @@ +// 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 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. + +mod _common; + +use std::mem::transmute; +use _common::validate; + +fn main() { + for bits in 0u32..(1 << 21) { + let single: f32 = unsafe { transmute(bits) }; + validate(format!("{:e}", single)); + let double: f64 = unsafe { transmute(bits as u64) }; + validate(format!("{:e}", double)); + } +} diff --git a/src/etc/test-float-parse/tiny-pow10.rs b/src/etc/test-float-parse/tiny-pow10.rs new file mode 100644 index 0000000000000..a01c6d5a07893 --- /dev/null +++ b/src/etc/test-float-parse/tiny-pow10.rs @@ -0,0 +1,21 @@ +// 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 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. + +mod _common; + +use _common::validate; + +fn main() { + for e in 301..327 { + for i in 0..100000 { + validate(format!("{}e-{}", i, e)); + } + } +} diff --git a/src/etc/test-float-parse/u32-small.rs b/src/etc/test-float-parse/u32-small.rs new file mode 100644 index 0000000000000..a4e8488e74529 --- /dev/null +++ b/src/etc/test-float-parse/u32-small.rs @@ -0,0 +1,19 @@ +// 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 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. + +mod _common; + +use _common::validate; + +fn main() { + for i in 0..(1 << 19) { + validate(i.to_string()); + } +} diff --git a/src/etc/test-float-parse/u64-pow2.rs b/src/etc/test-float-parse/u64-pow2.rs new file mode 100644 index 0000000000000..a31304d3f68aa --- /dev/null +++ b/src/etc/test-float-parse/u64-pow2.rs @@ -0,0 +1,28 @@ +// 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 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. + +mod _common; + +use _common::validate; +use std::u64; + +fn main() { + for exp in 19..64 { + let power: u64 = 1 << exp; + validate(power.to_string()); + for offset in 1..123 { + validate((power + offset).to_string()); + validate((power - offset).to_string()); + } + } + for offset in 0..123 { + validate((u64::MAX - offset).to_string()); + } +} From 15518a9c0c1c5f1774503c72072c73e64411d130 Mon Sep 17 00:00:00 2001 From: Robin Kruppe Date: Mon, 10 Aug 2015 23:14:30 +0200 Subject: [PATCH 7/7] Mention that the fast path is broken without SSE. --- src/libcore/num/dec2flt/algorithm.rs | 14 +++++++++++--- src/libcoretest/num/dec2flt/mod.rs | 7 +++++++ 2 files changed, 18 insertions(+), 3 deletions(-) diff --git a/src/libcore/num/dec2flt/algorithm.rs b/src/libcore/num/dec2flt/algorithm.rs index 97019090b56c7..f166bb9b3eb8d 100644 --- a/src/libcore/num/dec2flt/algorithm.rs +++ b/src/libcore/num/dec2flt/algorithm.rs @@ -21,9 +21,8 @@ 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?) +// We simply store the best approximation for *all* exponents, so the variable "h" and the +// associated conditions can be omitted. This trades performance for a couple kilobytes of space. fn power_of_ten(e: i16) -> Fp { assert!(e >= table::MIN_E); @@ -37,6 +36,15 @@ fn power_of_ten(e: i16) -> Fp { /// /// This is extracted into a separate function so that it can be attempted before constructing /// a bignum. +/// +/// The fast path crucially depends on arithmetic being correctly rounded, so on x86 +/// without SSE or SSE2 it will be **wrong** (as in, off by one ULP occasionally), because the x87 +/// FPU stack will round to 80 bit first before rounding to 64/32 bit. However, as such hardware +/// is extremely rare nowadays and in fact all in-tree target triples assume an SSE2-capable +/// microarchitecture, there is little incentive to deal with that. There's a test that will fail +/// when SSE or SSE2 is disabled, so people building their own non-SSE copy will get a heads up. +/// +/// FIXME: It would nevertheless be nice if we had a good way to detect and deal with x87. pub fn fast_path(integral: &[u8], fractional: &[u8], e: i64) -> Option { 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, diff --git a/src/libcoretest/num/dec2flt/mod.rs b/src/libcoretest/num/dec2flt/mod.rs index bd8cfc74f0c63..b7ef956055e29 100644 --- a/src/libcoretest/num/dec2flt/mod.rs +++ b/src/libcoretest/num/dec2flt/mod.rs @@ -90,6 +90,13 @@ fn zero() { test_literal!(1e-500); } +#[test] +fn fast_path_correct() { + // This number triggers the fast path and is handled incorrectly when compiling on + // x86 without SSE2 (i.e., using the x87 FPU stack). + test_literal!(1.448997445238699); +} + #[test] fn lonely_dot() { assert_eq!(to_f64("."), Ok(0.0));