mikros/rust
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

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

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

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

# Correctness

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

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

# Limitations

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

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

# Performance

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

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

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

# Binary size

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

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

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

# Future Work

## Directory layout

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

## Performance

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

## Binary size

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

## Other bases

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

Fixes #24557
Fixes #14353

r? @pnkfelix

cc @lifthrasiir @huonw 

[paper]: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.45.4152
[deprecate-radix]: https://internals.rust-lang.org/t/deprecate-f-32-64-from-str-radix/2405
This commit is contained in:
bors 2015-08-13 13:29:38 +00:00
commit bb954cfa75
29 changed files with 3828 additions and 18 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

134
src/etc/dec2flt_table.py Normal file
View File

@ -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 <LICENSE-APACHE or
# http://www.apache.org/licenses/LICENSE-2.0> or the MIT license
# <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your
# option. This file may not be copied, modified, or distributed
# except according to those terms.
"""
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 <LICENSE-APACHE or
// http://www.apache.org/licenses/LICENSE-2.0> or the MIT license
// <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your
// option. This file may not be copied, modified, or distributed
// except according to those terms.
"""
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()

26
src/etc/test-float-parse/_common.rs Normal file
View File

@ -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 <LICENSE-APACHE or
// http://www.apache.org/licenses/LICENSE-2.0> or the MIT license
// <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your
// option. This file may not be copied, modified, or distributed
// except according to those terms.
use std::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();
}

27
src/etc/test-float-parse/few-ones.rs Normal file
View File

@ -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 <LICENSE-APACHE or
// http://www.apache.org/licenses/LICENSE-2.0> or the MIT license
// <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your
// option. This file may not be copied, modified, or distributed
// except according to those terms.
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());
}
}
}
}

21
src/etc/test-float-parse/huge-pow10.rs Normal file
View File

@ -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 <LICENSE-APACHE or
// http://www.apache.org/licenses/LICENSE-2.0> or the MIT license
// <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your
// option. This file may not be copied, modified, or distributed
// except according to those terms.
mod _common;
use _common::validate;
fn main() {
for e in 300..310 {
for i in 0..100000 {
validate(format!("{}e{}", i, e));
}
}
}

39
src/etc/test-float-parse/many-digits.rs Normal file
View File

@ -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 <LICENSE-APACHE or
// http://www.apache.org/licenses/LICENSE-2.0> or the MIT license
// <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your
// option. This file may not be copied, modified, or distributed
// except according to those terms.
#![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<R: Rng>(n: u32, range: &mut Range<u32>, 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
}

32
src/etc/test-float-parse/rand-f64.rs Normal file
View File

@ -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 <LICENSE-APACHE or
// http://www.apache.org/licenses/LICENSE-2.0> or the MIT license
// <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your
// option. This file may not be copied, modified, or distributed
// except according to those terms.
#![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;
}
}
}

399
src/etc/test-float-parse/runtests.py Normal file
View File

@ -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 <LICENSE-APACHE or
# http://www.apache.org/licenses/LICENSE-2.0> or the MIT license
# <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your
# option. This file may not be copied, modified, or distributed
# except according to those terms.
"""
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()

29
src/etc/test-float-parse/short-decimals.rs Normal file
View File

@ -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 <LICENSE-APACHE or
// http://www.apache.org/licenses/LICENSE-2.0> or the MIT license
// <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your
// option. This file may not be copied, modified, or distributed
// except according to those terms.
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));
}
}
}

23
src/etc/test-float-parse/subnorm.rs Normal file
View File

@ -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 <LICENSE-APACHE or
// http://www.apache.org/licenses/LICENSE-2.0> or the MIT license
// <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your
// option. This file may not be copied, modified, or distributed
// except according to those terms.
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));
}
}

21
src/etc/test-float-parse/tiny-pow10.rs Normal file
View File

@ -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 <LICENSE-APACHE or
// http://www.apache.org/licenses/LICENSE-2.0> or the MIT license
// <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your
// option. This file may not be copied, modified, or distributed
// except according to those terms.
mod _common;
use _common::validate;
fn main() {
for e in 301..327 {
for i in 0..100000 {
validate(format!("{}e-{}", i, e));
}
}
}

19
src/etc/test-float-parse/u32-small.rs Normal file
View File

@ -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 <LICENSE-APACHE or
// http://www.apache.org/licenses/LICENSE-2.0> or the MIT license
// <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your
// option. This file may not be copied, modified, or distributed
// except according to those terms.
mod _common;
use _common::validate;
fn main() {
for i in 0..(1 << 19) {
validate(i.to_string());
}
}

28
src/etc/test-float-parse/u64-pow2.rs Normal file
View File

@ -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 <LICENSE-APACHE or
// http://www.apache.org/licenses/LICENSE-2.0> or the MIT license
// <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your
// option. This file may not be copied, modified, or distributed
// except according to those terms.
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());
}
}

361
src/libcore/num/dec2flt/algorithm.rs Normal file
View File

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

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

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

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

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

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

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

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

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

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

File diff suppressed because it is too large Load Diff

153
src/libcore/num/flt2dec/bignum.rs
View File

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

2
src/libcore/num/flt2dec/strategy/dragon.rs
View File

@ -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];

6
src/libcore/num/flt2dec/strategy/grisu.rs
View File

@ -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;

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

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

1
src/libcoretest/lib.rs
View File

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

181
src/libcoretest/num/dec2flt/mod.rs Normal file
View File

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

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

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

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

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

89
src/libcoretest/num/flt2dec/bignum.rs
View File

@ -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));

2
src/libcoretest/num/flt2dec/strategy/dragon.rs
View File

@ -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]

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

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