Decimal to floating point conversion

[rkruppe]

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 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+ 3s. 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 FIXMEs 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. 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

[Gankro]

CC @lifthrasiir

[eefriedman]

The decimal exponent is too big. This limitation is shared with the old code, though it accepts any exponent that fits into an i64 while this version stops at i64::MAX / 10, since it needs to do arithmetic with the exponent.

This seems like something which could be handled trivially: if the exponent overflows, either the correct result is inf or 0, or the string has a stretch of zeros way too long to fit into memory.

[rkruppe]

That's a good point. Exponent and number of decimal digits used to be a smaller integer type and I didn't stop to reconsider when I switched to 64 bit integers. However, exploiting this insight requires restructuring some code: The exponent is parsed in the parsing module, which isn't expected to produce any floats. I'll sleep over the least hacky way to integrate this.

A related aspect is that the integral and fractional parts really can't be larger than 1.8 exabyte, which is the limit now. I'll remove those guards when I do the other change.

[eefriedman]

However, the third bullet point would have rejected some sensible inputs, such as 1.7976931348623157e-324 (the smallest positive subnormal number, represented with 17 fractional digits; note that the same float is usually expressed as 5e-324).

Umm, 1.7976931348623157e-324 ? Are you sure you didn't mean some other number?

[eefriedman]

It would probably be a good idea to deprecate f32/f64::from_str_radix. (We don't want to keep them around given that they don't actually work correctly.)

[eefriedman]

This should explicitly document that it only handles floats with a positive sign bit.

[rkruppe]

Done.

[eefriedman]

This might be a silly question... but does fast_path() work correctly on Linux x86-32 with SSE disabled? (The Clinger paper calls out a similar scenario in section 9.)

[rkruppe]

Umm, 1.7976931348623157e-324 ? Are you sure you didn't mean some other number?

Gah! The decimal part belongs to another boundary case. The correct number is 4.9406564584124654e-324 Edited.

It would probably be a good idea to deprecate f32/f64::from_str_radix. (We don't want to keep them around given that they don't actually work correctly.)

Good to hear. I'd like to collect some more support before I go ahead and do it, and I'm not sure if it needs to go into this PR.

This might be a silly question... but does fast_path() work correctly on Linux x86-32 with SSE disabled? (The Clinger paper calls out a similar scenario in section 9.)

Not silly at all! Unfortunately the VM on which I'd normally test this currently has... issues. If anyone could compile a short program doing a float multiplication with a 32 bit rust and -C target-feature="-sse" and post the assembly, that'd be great. x87 can round correctly by setting the control word correctly, but GCC for example never does this. I'd have to check if LLVM allows doing so and if rustc instructs it to.

[rkruppe]

Lifted the unnecessary restrictions @eefriedman pointed out. I immediately squashed that into the last commit, I hope that's okay.

[eefriedman]

Do you need to trim zeros here?

[rkruppe]

Indeed, thanks! That would have been sensible even when rejecting large exponents.

[eefriedman]

Might want to explicitly note why the math involving e in this file can't overflow; it's a bit subtle.

[rkruppe]

Done (in the long comment at the start of the file, and also in parse.rs).

[eefriedman]

You can "round" numbers with too many digits: once you have too many digits, the actual values of the digits cease to matter. Your bigints need to be a bit larger to make that work, but it's the same order of magnitude; glibc uses the formula "1 + ((DBL_MANT_DIG - DBL_MIN_EXP + 2) * 10) / 3" (==3587) for the number of necessary bits.

[rkruppe]

Interesting. Do you have a paper or other write-up explaining this in detail? It is not at all obvious to me how this works, or what the corresponding algorithm looks like.

[rkruppe]

Thinking about it some more, I find it hard to believe that this is always strictly correct. Consider the exact decimal representation (this always exists and is terminating) of some float with an even mantissa. Add half an ULP. This still rounds to the same float as we started with because the mantissa was even. Now add 2^-10000 or some other sufficiently small power of two: Suddenly we need to round up, and it's not visible without inspecting thousands of additional bits! The tiny power of two does not affect the first couple thousand decimal digits either. Is there some magic trick to correctly evaluate the decimal digits you're not putting into the bigint?

[eefriedman]

There's an overview of glibc's algorithm at http://www.exploringbinary.com/how-glibc-strtod-works/ .

In terms of correctness, what I meant by "round" is to treat the existence extra digits like a sticky bit. For example, take the exact decimal representation of half of the smallest denormal double. If you convert that exact number, it gets rounded down to zero. If you append any combination of digits (assuming they aren't all zero), it gets rounded up. The values of the digits don't actually matter.

[rkruppe]

That sounds sensible. However, from my experience over the last few weeks, I predict that actually implementing this algorithm (like all the algorithms that I did implement) will be fraught with many small yet tricky issues, making it quite time-consuming to polish and bring up to production quality. Therefore I suggest that this too is filed under future work. The remaining restrictions of the code in this PR should not be a problem for any practical inputs.

[eefriedman]

That sounds fine.

[arielb1]

The result of dec2flt depends on the order of the input with respect to the numbers half-between floating points. These are numbers of the form (2^53 + 2u + 1) * 2^-s for some 0 <= u < 2^52 for normals, and (2u+1)*(2^-1075) for denormals. These numbers all become integers when multiplied by 10^1075, and therefore have less than 1075 nonzero digits (to be more precise, after multiplying by 10^s, you are left with numbers smaller than 2^54*5^s, so smaller than 10^[log(2,10)*54 + log(5,10)*1075] < 10^768, so with no more than 768 nonzero digits) - beyond that last digit the only thing that matters is the tie-breaker sticky bit.

For the record, one such number of the maximum possible length is halfway between the largest subnormal and the smallest normal,

2.22507385850720113605740979670913197593481954635164564802342610972482222202107694551652952390813508791414915891303962110687008643869459464552765720740782062174337998814106326732925355228688137214901298112245145188984905722230728525513315575501591439747639798341180199932396254828901710708185069063066665599493827577257201576306269066333264756530000924588831643303777979186961204949739037782970490505108060994073026293712895895000358379996720725430436028407889577179615094551674824347103070260914462157228988025818254518032570701886087211312807951223342628836862232150377566662250398253433597456888442390026549819838548794829220689472168983109969836584681402285424333066033985088644580400103493397042756718644338377048603786162277173854562306587467901408672332763671875E-308

That number is (2^53-1)/2^1075, and has 307 zeros and 768 significant digits.

[pnkfelix]

I'm fine with leaving this extension for future work (i don't want this PR to be delayed further than it already has been, due to my tardiness in reviewing it).

[Gankro]

re deprecation: kill it with fire

[pnkfelix]

Okay, done with review! I didn't see anything major, and I trust @rkruppe will be good about addressing any of my nits as @rkruppe sees fit.

So, r=me once it is ready to land. (I'll be out of town for a week, feel free to ping other people in #rust-internals for r=pnkfelix marker once its ready.)

[rkruppe]

So I finally got around to testing whether the fast path is correct on x86 without SSE (thanks, @eefriedman!). Turns out the exact same problem described in the paper (in the context of a Motorola 68881/68882 floating point coprocessor) also applies to more common hardware. 1.448997445238699 is rounded incorrectly, to the float whose shortest decimal representation is 1.4489974452386991 (off by one ULP). Bummer. I'll do something about this before merging, but I don't know yet what.

[rkruppe]

Now that PR #27647 has been merged and the fast path brokenness is at least somewhat addressed (see last commit), I declare this PR officially ready.

The resolution is somewhat unsatisfying and maybe I'll find a better way in the future, but that shouldn't stop the very real improvements (even on platforms where the fast path is technically "broken") from landing.

[eefriedman]

Isn't (target_arch="x86", target_os="ios") targeting the iOS simulator? In that case, the fact that it doesn't default to SSE2 enabled is probably an accident.

[rkruppe]

Yes, it's targeting the iOS simulator. Since @cmr brought that up and the very same @cmr later called the fact that x86 MSVC does not use SSE an accident, I can only assume that disabling SSE on the iOS simulator is not an accident.

Edit: Actually @cmr just confessed to having no idea either.

[rkruppe]

Since there are no more in-tree targets that disable SSE, "deal with it" now means "explain that it's theoretically broken and add a test that will expose said brokenness if some unfortunate soul is forced to compile without SSE".

[eddyb]

@bors r=pnkfelix

[bors]

📌 Commit 15518a9 has been approved by pnkfelix

[bors]

⌛️ Testing commit 15518a9 with merge bb954cf...

[bors]

☀️ Test successful - auto-linux-32-nopt-t, auto-linux-32-opt, auto-linux-64-nopt-t, auto-linux-64-opt, auto-linux-64-x-android-t, auto-mac-32-opt, auto-mac-64-nopt-t, auto-mac-64-opt, auto-win-gnu-32-nopt-t, auto-win-gnu-32-opt, auto-win-gnu-64-nopt-t, auto-win-gnu-64-opt, auto-win-msvc-32-opt, auto-win-msvc-64-opt

[IvanUkhov]

@rkruppe, thank you for this work!

[pnkfelix]

\o/

[donbright]

"absolutely positively completely 100% accurate "

i question the basic assumption that base-10 numbers can be converted into base-2 with 100% accuracy. there will be information lost in the transfer, to me it seems that would result in accuracy of some number less than 100%. and that number should be measurable.