High Performance Correctly Rounded Math Libraries for 32-Bit Floating Point | Hacker News Reader

High Performance Correctly Rounded Math Libraries for 32-Bit Floating Point

blog.sigplan.org

107 pointsingve73 comments

This work essentially allows one to go beyond the IEEE754 standard which mandates correctly-rounded elementary operations +, -, x and / : they provide implementations of libm functions (cos, sin, log, exp, etc.) that are also correctly rounded.

This is really nice, not necessarily because anyone cares about the correctness of the last bit in a float. But more importantly, because since there is only one "correct" answer, requiring correct rounding means that those functions are then fully specified and deterministic. No more getting slightly different result when changing platform, OS, or even just updating libm (some transcendental functions are still actively being improved, thus changed, in glibc!). This is amazing for reproducibility.

Since the performance is good (even better than existing libm implementations, they claim), this is a true win-win.

This is obviously made possible by the fact that one can easily enumerate all 32-bit floats. It would be amazing to have something similar, even at a big performance cost, for 64-bit doubles (MPFR can do it of course but the perf cost there is truly massive). Things are much harder with 64 bits, unfortunately.

> No more getting slightly different result when changing platform, OS, or even just updating libm

Optimizing compilers are usually used with flags that cause them to break IEEE754.

Still, it would be nice if results from low-optimization runs might be bit-exact.

Yes, indeed it is the case with fused multiply-add operations, as I wrote below in another thread. By default at -O3, unfortunately, gcc and icc both take responsibility for deciding whether or not to use FMA instructions when they see (a x b + c) expressions. As a result we can even see slightly different result across compiles!!

However, as far as I understand, specifying a proper language standard (for example -std=c99, -std=c++14, etc.) restores sane behavior (see -ffast-math, and more specifically -fexcess-precision).

Are there other cases where -O3 breaks IEEE754 on modern setups?

Fast math flags are actually far more complicated in C/C++ compilers than just a simple -ffast-math flag, although I suppose most users may not be acquainted with the degree of control possible. In LLVM, for example, fast math flags are actually a set of 7 independently-toggled flags (no-signed-zero, no-infinity, no-NaN, allow FMAs, allow reassociation, allow approximate reciprocal, and allow approximate functions), that at the C/C++ level can be toggled either by command-line options or on a per-scope basis. ICC actually has two different levels of fast math [1].

The dark underside of IEEE 754 is that full support for it also requires features that most languages don't provide access for: rounding mode and exception handling (aka sticky bits). In C/C++, using these features correctly requires STDC FENV_ACCESS support which is only in ICC and recent clang (neither GCC nor MSVC support it). And there's usually two hardware bits for handling denormals (flush-to-zero and denormals-are-zero) that aren't in IEEE 754, which tend to be on by default when you're compiling for some targets (most notably GPUs).

[1] https://software.intel.com/content/www/us/en/develop/documen... -- note fp-model=fast=[1|2]

What is a "per scope basis"?

You can change it between any pair of curly braces.

E.g.:

  void foo(int N) {
    for (int i = 0; i < N; i++) {
      #pragma STDC FENV_ACCESS ON
      // This applies only to the for loop...
    }
    // But not here, for example!
  }

> Yes, indeed it is the case with fused multiply-add operations, as I wrote below in another thread

This is not always true. Some FMA operations round after the multiply, and some do not.

> Optimizing compilers are usually used with flags that cause them to break IEEE754.

For languages that don't mandate IEEE-754. There are plenty of languages that do, like Java, JavaScript, C#, and WebAssembly.

It's worth noting that usually such optimizations are things like enabling algebraic reassociation and commutation (of primitives like +, -, /, ==, etc) operations, fused-multiply add, and others. It's really not clear to me how much disabling all those optimization costs in the real world.

Things get crazy when compilers want to do major reorganization of loops, optimizing stencils, polyhedral optimizations. I am not an expert here, but I think these rely on properties of float arithmetic that don't hold in all cases, i.e. they rely on UB.

> It's really not clear to me how much disabling all those optimization costs in the real world.

I assure you that SPECcpu scores are a lot lower, if that's what you care about. And yes, it's extremely helpful to do major reorganization of loops if the program has multi-dimensional loops in the wrong order.

Apple had a library that did that for 80-bit floating point, for both the 6502 (!) and 68000: SANE (https://en.wikipedia.org/wiki/Standard_Apple_Numerics_Enviro...)

And yes, performance wasn’t good, especially compared to the x87 that IBM PCs used. It did claim to round the last bit correctly, though, where early x87s could have large errors (https://randomascii.wordpress.com/2014/10/09/intel-underesti...)

> Things are much harder with 64 bits, unfortunately.

Solved and published over 15 years ago. The underlying methods were published over 20 years ago (V. Lef`evre, J.M. Muller, and A. Tisserand. Towards correctly rounded transcendentals. and Lefvre's associated PhD thesis).

https://hal-ens-lyon.archives-ouvertes.fr/ensl-01529804/file...

The Correctly Rounded Math Library (CRLIBM) is available for several languages.

CRLIBM is great, but it is significantly slower than alternatives since it requires very high internal precision. The cool thing with this paper is that it gives correct rounding without as much of a cost for 32 bit.

That's a bit glib. I don't believe that the bounds are fully proven for all the transcendentals (and I believe that exponentiation is only "proven" over restricted subdomains).

Single precision floating point is a small enough domain that it can be fully enumerated--you can prove that everything rounds correctly and find the problematic cases. 24-bits just isn't that large to modern computers.

Double precision floating point is still too large to fully enumerate. So, there could be an exception lurking in there that goes out a long way.

In reality, it seems that rounding requires about 2n + log2 n bits and never 3n bits (I don't think we've ever found an exception to this). However, if you can prove that, you should submit it as it would be a significant advance in the world of mathematics.

In general, there shouldn't be a rule like 3n bits. For random functions, 1 in 8^n values should require more than 3n bits. Since n bit inputs give 2^n possible values, this implies that for bit length n, there is a 4^-n chance that any given function will require more than 3n bits.

Except that this isn't the case. We know empirically because we have enumerated the single precision space completely.

So, the problem isn't bounded by simple probabilities. There are other properties involved that cause it to diverge from the "naive" probabilistic estimates. What those properties are we don't know yet.

The "rounding" problem is called the "Table Maker's Dilemma" and it descends into some pretty fundamental transcendental number theory.

This is not an "easy" problem. Proof that PI was transcendental is relatively modern (1885).

This doesn't contradict what I said. I said that the odds of a function needing more than 3n digits was 4^-n. For single precision, that means you would expect almost all functions to be OK with 3n bits (1 in 2^46 wouldn't be).

if it was a randomly selected function, which it is not.

In practice the low-order bits of transcendental functions are pretty much indistinguishable from random for _most_ inputs. You do have to be very careful close to especially nice values though, because this falls apart catastrophically there (e.g. logarithms close to 1)--it should be pretty obvious from the Taylor expansions why this happens. The good news is that when this does happen, almost by construction you usually don't end up very close to an exact halfway point for rounding.

> In practice the low-order bits of transcendental functions are pretty much indistinguishable from random for _most_ inputs.

In my experience, low order bits of transcendental functions make terrible rng sources (even discounting their runtime). Your example of the logarithm is great -- for large inputs, you need large changes to flip even low order bits.

If you're just looking for hard-to-round inputs, you can also enumerate the double precision space completely with some cleverness and modest patience.

Elkies [0] computes a piecewise linear approximation to his (number-theoretic) function of interest and does a lattice reduction on each piece to filter out the vast majority of the space. The same approach, works almost everywhere on all of the transcendental functions I've tried, and I believe the CRlibm group uses it.

[0] https://arxiv.org/pdf/math/0005139.pdf

Oooh, that's an interesting take.

I've read that, if one is going to use floats, then one should eschew 32-bit floats, and go with 64-bit floats. 80-bit floats seem only useful for very special applications.

I converted a satellite orbital prediction program that used doubles (64-bit floats) to use long doubles (80-bit floats) but only saw a difference in the 11th decimal place, with lots of sin/cos/atan/sqrt involved in the calculations. With 32-bit floats it was a joke: the accumulated error made the 'predictions' useless.

The basic problem with 32-bit is you have to be quite careful with your arithmetic so you don't lose significance. It's almost as bad as using Q-notation. https://en.wikipedia.org/wiki/Q_(number_format)

The main use for 80-bit floats is probably to confirm that when you switch from 64-bit to 80-bit, nothing changes.

This is a very important use, though!

Errors grow exponentially in most floating-point operations.

The exception is cancellation error, which can suddenly cause the error rate to spike dramatically. But if you sort your numbers before adding/subtracting, you can often negate this overwhelming cancellation error.

The exponential growth of other floating point errors is almost unstoppable however. The only real way to stop it is to design an algorithm that operates with fewer steps.

---------

64-bit has a much smaller initial error (1/2^53), while 32-bit has (1/2^23) initial errors... starting with 30-magnitudes more error than double-precision.

Kahan summation[1] may be better than just doing a sorted sum & doesn't require you have all the numbers up-front to sort (i.e. running sum). Although I think even it benefits from doing it in sorted order if you do have the full history available & the wikipedia article lists a parallelization technique you can leverage in that case. Klein sum might be better but I haven't investigated.

EDIT: Apparently there's a nice Rust crate for doing this stuff: https://crates.io/crates/accurate

[1] https://en.wikipedia.org/wiki/Kahan_summation_algorithm#Alte...

One thing to consider here is XSum. It's a C library that has surprisingly fast exact floating point sums.

Kahan summation is basically just synthesizing a larger floating point type (with 46/106 bits of mantissa for 32/64-bit base type) and then using that to do the sum, though. You still lose precision if your errors exceed that new, larger margin. Eg try adding 34000 + 56 + .078 - 56 - 34000, carrying two digits per float.

It really depends what you're using the floats for. For ML applications? 16 bit floats do the job with some trickery (quantization).

For scientific computing? Generally 64 bit floats give you the needed precision. But if you're chaining a lot of operations that exacerbate floating point error you may need more.

For usage in trading? Again, depends on what computations you're doing. There's no "hard and fast" rule.

32-bit floats are sufficient and desired for performance in many cases, like graphics or audio. 64-bit floats are not enough in some cases, like galaxy-scale coordinates or nanosecond precision over long time periods or some currency calculations. There's no single best answer, you have to know what your requirements are.

> 32-bit floats are sufficient and desired for performance in many cases, like graphics or audio

I was googling this recently, when deciding whether to use floats or doubles in a simple game I was writing. Some answers suggested that (32-bit) floats aren't necessarily any faster in real-world applications on modern hardware. If you or anyone else feels like elaborating on this I'd be interested.

For scalar computations 32 bit floats are not going to be faster except maybe divisions. There are second order effects due the reduced memory usage.

For anything that is vectorizable or otherwise amenable to SIMD optimizations you can get up to 2x speedup.

All the transcendental functions are also an extra 2-4x slower for Float64 since you need higher order approximations, and the compensated arithmetic is slower.

For CPU cycle counts it's not that different, for GPU cycle counts it's something in the range of 50x FLOPS difference for consumer GPUs so avoid it in shaders and whatnot if you can. There is also the matter of SIMD on the CPU, you can generally get twice the SIMD rate by halving the bit width of the data which is quite the gain.

Outside raw computation rate though there is also the matter of cache and memory bandwidth. You'll use 2x with doubles which can really bog things down as even if the CPU can crunch the calculation quickly it could spend a lot of time waiting for the data to crunch to arrive.

All that being said if you know your game is going to be simple then it probably doesn't matter one way or the other even if it does come at a big performance cost.

The big win is memory. Half the memory/cache usage and half the bandwidth. For SIMD you'll want 32 bit floats. But in many cases instead of reaching for SIMD you should go to compute shaders, which are definitely going to be 32 bit. And of course any graphics shaders will be using 32 bit floats.

Games generally don't "accumulate" in floats for extended periods, unlike typical scientific applications where they might be using iterative solvers, simulating processes, or inverting huge matrices.

Instead, games tend to take some input parameters, use those to transform a read-only data set, render a frame, and then throw everything away. The next frame starts from scratch. There's no accumulation.

In my experience, typically only the "player camera angle" is kept around and accumulated into. The precision errors here are negligible. Nobody cares if the viewport is 0.00001 degrees to the left or right of where it "should" be.

Most other parameters that are kept frame-to-frame tend to be "world simulation" stuff. In that case, consistency is so important for network codes that most games use fixed point or integers only. This avoids any kind of discrepancies when using delta compression or multiplayer where people have different CPUs.

It seems to be too late to edit, so I'll post here rather than clog up the thread with individual replies: thanks, I appreciate all the explanations. I will have to look up some details but it seems the basic takeaway is 'if in doubt, just use floats'.

> I've read that, if one is going to use floats, then one should eschew 32-bit floats, and go with 64-bit floats. 80-bit floats seem only useful for very special applications.

This is what Kahan has to say on the matter:

For now the 10-byte Extended format is a tolerable compromise between the value of extra-precise arithmetic and the price of implementing it to run fast; very soon two more bytes of precision will become tolerable, and ultimately a 16-byte format ... That kind of gradual evolution towards wider precision was already in view when IEEE Standard 754 for Floating-Point Arithmetic was framed

(The quote is apparently from 1994 unpublished manuscript, so I can't find the source for it)

More than 25 years later and we are rarely using even 80 bit doubles, nevermind 128 bit ones.

Isn't that because 80-bit floats have been deprecated a decade or two ago in x86, and were never implemented in other architectures to begin with?

Both architectures implemented 80-bit floats - Intel x86 and Motorola 68k.

They were deprecated in x86 because they weren't widely used.

Depends on what you are doing. 8-bit floats are widely used in large neural network applications. (No that is not a typo. I mean 8-bit, not 80-bit.)

Really? I thought that for 8 bits, fixed point formats were much more common.

At that point, it would seem more sensible to just define a list of values that are represented by each possible 8-bit number, and use LUTs for operations. Why even torture this exponent and mantissa concept down to 8 bits when one can have a domain specific set of representable numbers that have exactly the properties that are needed.

I did exactly this - 8-bit floats. The application needed fairly accurate representations around 0; but didn't wildly swing. It was for an IMU that had to have the internal representation converted to 'regular units'. Ultimately, Q.15 won; but it was interesting. In fact, I encourage you to try it with e.g. 4 mantissa 4 exponent, and enumerate all values.

For a game boy (original) project I wanted to represent distance values in range 0.3..32 with a u8, but with higher precision for closer values, but approximately linear relationship for distances that are close to one another, with precision of ~0.05 around d=0 and <0.30 around d=32.

I ended up creating a set of values that follows different mathematical functions in different segments (exponential, linear, quadratic). As a float this wouldn´t have worked, because the precision halves for every new exponent, so the distances between adjacent values are not approximately linear.

I love everything about this project, but I am fairly skeptical of their benchmarking. In particular, they use a _lot_ of polynomial intervals for some of these functions, it's easy to accidentally construct microbenchmarks that miss the cache effects that fall out of that (the full set of coefficients for one of these math functions can easily take up a substantial fraction of your L1; for some use cases that won't matter, for some other scenarios it will be catastrophic).

Historically, high-performance libm implementations have gone to great pains in their internal benchmarking projects to avoid accidentally ignoring these effects and shipping functions that are significantly slower in real use. None of this is to say that they have necessarily fallen into this trap, but historically most "fast" math function implementations that use piecewise polynomials with high segment counts do (and none of it should take away from the accomplishment--even if they're really 2x slower than claimed, they're fast enough for most code to simply use it unconditionally and get the portability benefits).

There are ~30 year old legends in the client/server gamedev community about floating point differing across platforms, but I assumed they were ironed out by now. Wild to find out they haven't been.

Decades ago we were deeply concerned with syncronized simulations; the dominant replication method was lock-step input replication, or similar, and that made it very important that everyone's simulation was identical or bad things would happen.

Nowadays that's not an issue, as there's a slew of different ways of overcoming that and allowing clients to operate asynchronously; sometimes widely asynchronous.

Any good resources (or even search term hints) you could recommend on the more modern/async approaches? I’d love to read about them.

https://gafferongames.com/post/state_synchronization/

Also search for "game networking" + the following:

  - state replication
  - dead reckoning
  - client side prediction
  - networked physics

It's still super important for clientside prediction. If the client and server are getting different math results that's a desync. Now that can matter not at all or a lot depending on the game, but regardless it's still fundamental.

Nowadays, all major platforms follow IEEE754, which specifies exactly the results of elementary operations (+, -, x, /). You must get bit-for-bit reproducible results across different runs, compiles, libraries, OSes and hardware. Any difference is a bug.

However:

1. As the paper emphasizes, math library functions (sin, cos, log, etc.) are not covered by that standardization, at least not mandatorily and not in practice. This is the problem that the paper tackles.

2. In the bad old days of the x87, the FPU carried out computations on eight 80-bit registers. The problem is that whenever those registers spilled out to the stack, precision was truncated down to whatever type the language actually meant (64-bit double or 32-bit float). Because the register-spilling decision is mostly taken by the compiler, it was largely unpredictable. This resulted in code that was not reproducible, even across compiles! Since x87 was deprecated in favor of SSE/AVX, the precision of operations is always that of the type used in the code. Lower than 80-bits but much better in practice because of consistency and reproducibility.

3. The bad old days are making a comeback with FMA instructions (fused multiply-add). This is a single operation (a x b + c) and so it is correctly rounded to the closest representable number. This number can be different from ((a x b) + c) in which two roundings occur. If we give the compiler freedom to decide whether or not to fuse such expressions, we will again see non-reproducibility across compiles. Unfortunately, gcc and icc both do that by default at -O3 (see -ffast-math, and more specifically -fexcess-precision). Specifying a proper language standard (for example -std=c99, -std=c++14, etc.) restores sane behavior.

On the one hand, the ambiguity of FMA is kind of annoying, but on the other, FMA is a godsend for writing accurate arithmatic. One reason for this, is that the error of `ab` is `fma(a, b, -ab)` which makes compensated algorithms much faster to write. Secondly, polynomial evaluations using horner's method converge to .5 ULP (units in last place) with FMA, but to 1.5 ULP without. This is incredibly useful since polynomial approximations are at the heart of pretty much everything else.

I think Julia takes an interesting approach here of never reordering floating point computations unless you explecitly give permission using @fastmath (which acts locally instead of globally like a --fast-math in C/C++). This makes it much easier to use ieee specific tricks and be confident that the compiler won't mess with you, while still making it easy to tell the compiler to do what it wants.

> FMA is a godsend for writing accurate arithmatic.

True, and the small performance boost (on some platforms) is nice too

> never reordering floating point computations unless you explecitly give permission

Yes, requiring explicit code for FMA and arithmetic grouping/reordering in general seems like the sane approach.

> I think Julia takes an interesting approach here of never reordering floating point computations

WebAssembly as well. IEEE-754 is mandated and there are no unsafe compiler optimizations allowed. This is the only sane choice. It's just too hard to reason about programs otherwise. Nothing about a super-intelligent machine completely reorganizing your code in a way that subtly breaks it is good, IMHO.

x87 can store and load full 80-bit register values. The truncation you mention is optional.

> 2. In the bad old days of the x87, the FPU carried out computations on eight 80-bit registers. The problem is that whenever those registers spilled out to the stack, precision was truncated down to whatever type the language actually meant (64-bit double or 32-bit float). Because the register-spilling decision is mostly taken by the compiler, it was largely unpredictable. This resulted in code that was not reproducible, even across compiles! Since x87 was deprecated in favor of SSE/AVX, the precision of operations is always that of the type used in the code. Lower than 80-bits but much better in practice because of consistency and reproducibility.

even today it's still possible to end up with things in the x87 registers, I don't know of a compiler toggle to entirely and unequivocally disable it.

That seems like a fault of the language/compiler, not the 80but float. There are operations for storing the 80bit float in full precision, and I think Pascal even used real=80bit.

Code using "long double" on x86_64 will compile to x87 instructions (at least with GCC on Linux).

This can be triggered unintentionally. For example, string-to-double conversions in Boost use long doubles.

(I have been bitten by that when using valgrind: Valgrind emulates x87 80-bit registers using 64-bit doubles. As a result, some Boost-using code will have a different behavior under valgrind.)

Yup, exactly. You can imagine that a lot of collision code uses tan for example.

I know this is off topic, but the inconsistent use of "i.e." and "e.g." is particularly bad in this article. If you are not sure of the difference in your own writing, use "that is" in stead of "i.e." and use "for example" instead of "e.g."

Notice in this article, in the section entitled "What is the correct result of f(x)?" they say "e.g. ln(x)" in one case and "i.e. ln(x)" three sentences later. Well which is it? It is not just that I am a curmudgeon. Precision matters. When they say, for example, "Everyone uses math libraries (i.e. libm)" do they mean "Everyone uses the libm math library" or, more likely, "Everyone uses math libraries such as libm". Yes, the resolution in this case is fairly obvious, but I would get much more out of the article if I was not distracted by these slight inadvertent ambiguities and misuses.

The nemonic I was taught was "In Explanation" and "Example Given."

I always just think of them as “in effect” and “(for) example”, which is their direct expansion. But that distinction makes sense to me because I use the phrase “in effect” in speech as well sometimes.

It hadn’t occurred to me that this is unusual - I suppose it is? I wonder if it’s a aus / US cultural difference? (I’m Australian) or maybe something I picked up years ago as a grad student?

I just think “egsample” and then i.e. must be the other one.

Its unfortunate that they don't mention the prior work of the crlibm project in this domain.

https://gforge.inria.fr/scm/browser.php?group_id=5929&extra=...

Github mirror: https://github.com/taschini/crlibm

Python bindings: https://pypi.org/project/crlibm/

Julia bindings: https://juliapackages.com/p/crlibm

crlibm is discussed in their paper: https://people.cs.rutgers.edu/~sn349/papers/rlibm32-pldi-202...

Based on their press material, you'd think that nobody else has done this before at any precision. You can see the resulting confusion in this very forum. That kind of behavior is disingenuous at best.

How does this compare with sleef[0], which also provides high performance (as well as vectorization), and implements all the same routines as rlibm (among others) with ≤1 ulp error?

0. https://sleef.org/

It will be a good bit slower. The difference between 1 ULP and .5 ULP is pretty big. Also the method of using a bunch of different polynomials won't vectorize especially well.

I'm a libc author and I want to use this. Then I clicked on the repository and there wasn't any code implementing these libm routines. Instead it just had tools and they told me I need to install a 6gb docker container to use them.

> Then I clicked on the repository and there wasn't any code implementing these libm routines

Whats this then?

https://github.com/rutgers-apl/rlibm-32/tree/main/source/flo...

Thanks!

I love that they have an implementation for posits too.

Reply on news.ycombinator.com