Floating Point, Simply

Intro

At this point, I’ve been a heavy user of IEEE 754 floats for well over 20 years (!) on CPUs and GPUs through many hardware generations, so I thought it might be interesting to jot down how I work with and think about floating point values in general.

This isn’t a novel topic by any means, and much of the literature that exists dates back decades. Some of the oldest books on my shelf pertaining to numerical computation contain ideas that are still relevant today, Hamming’s Numerical Methods for Scientists and Engineers from 1986 being just one such example.

And despite all that writing, I think most people would agree that floating point can often be surprisingly tricky. While superficially easy to use, there are tons of seemingly invisible footguns haunting every operation.

My goal here isn’t to expand on every footgun or corner case ever. Instead, this post is about general rules that I’ve more or less stumbled on to live with floating point without too much pain.

Numerical robustness

Here’s the first tip. When coding anything to do with floating point, keep numerical robustness top of mind.

To me, something is numerically robust if small perturbations in the input result in small perturbations in the output. Simple definition, but honestly a very powerful mode of thinking.

Case Study: Quadratic Formula

Let’s take the quadratic formula as a quick example:

This is a formula that I consider not numerically robust. Despite being a well-known and widely used formula, it has a lot of properties that should basically trigger you immediately if you’re the type that works with floating point a good deal.

The main things that should jump out to you are:

• The number of roots can change suddenly when the sign of the discriminant changes.
• The function blows up when the value of a crosses 0 in either direction.

There are other precision and overflow issues to consider here, but those two issues are, in my mind, the main ones. In both cases, a small perturbation in the input can result in a massive change in the outcome! Whenever you notice this sort of thing, you should be prepared to code carefully to ward away nasty surprises later.

A “proper” quadratic solver probably needs to start off by handling all the special cases where a, b, and/or c are zero. From there though, if you were to implement such a solver, you have a number of choices. Numerical robustness isn’t a binary state of being (i.e. something isn’t numerically robust or not); it’s more of a spectrum. Often times, achieving “perfect numerical robustness” can be quite expensive to achieve, so you have to make an engineering decision: “for my use case, how much am I willing to pay for robustness?”

Going back to the quadratic solver, let’s consider some options. Some of the other things to notice are:

• If b*b is large relative to 4*a*c, then sqrt(b^2-4*a*c) will be approximately equal to b, and the numerator will have a loss of precision when the positive root is used.
• If b*b is approximately equal to 4*a*c, we have a different problem where the discriminant itself will experience catastophic cancellation.
• It’s possible that b*b - 4*a*c might overflow to some massive value from which the expression never recovers (even if the roots themselves live in perfectly valid floating point ranges). Similarly, an underflow to a subnormal is possible.

Dragons at every corner indeed. Let’s consider each case and decide if it’s something we need to care about.

b*b >> 4*a*c

In this case, unfortunately, there are so many choices of b, a, and c for which this condition is true so there’s no way around addressing this problem somehow.

The most straightforward solution is actually an algebraic one. In pretty much every expression that contains the addition between a value and a square root, we can multiply the numerator and the denominator by the conjugate of the square root to obtain a different expression with the root sign reversed.

And now we’re left with a an equivalent form that exchanges the sign of the square root. Because algebra has no “numerical precision” issues to speak of, we now have the freedom to choose the first form or this conjugate so-called citardauq form depending on the sign of b.

I’m showing the whole calculation here in gory detail because this type of manipulation extends beyond the quadratic formula to plenty of other cases also.

b*b roughly equal to 4*a*c

In this case, we can say a bit more about the shape of the parabola. Here, the discriminant is approaching 0, so we’re in a realm where the inflection point of the parabola is nearly touching the x-axis.

Unlike the first case, there isn’t any algebraic magic we can use to help us. The signs of both b*b and -4ac are fixed. What we could consider is:

• Computing the determinant at higher precision with double values.
• Use fma for a more accurate result.
• Both.

The FMA solution is technically “free” on machines that support it, although you may lose determinism since the exact precision afforded by FMA instructions isn’t stipulated by any IEEE specification.

This is an area where you, the implementer should decide what the behavior ought to be. If you’re using the solver for something that demands the utmost precision, then pay for it. If this solver is used in some code for a particle system, or inverse kinematics, it’s probably not worth addressing. I wouldn’t suggest using an off-the-shelf solution though (at least, not without consulting the implementation), since the implementation may either be outright wrong (surprisingly common) or the implementation might make the wrong tradeoffs.

Overflow and underflow

The last cases to consider: it’s possible that the expression b*b - 4*a*c explodes to some massive value and all precision is crushed into “infinity.” Conversely, it’s possible that the expression underflows and we lose precision by crushing the value into the smallest epsilon value representable.

If this happens, we technically can use algebra to a certain degree. The equation we’re trying to solve is:

Notice that we can multiply both sides of the equation by any constant and the values of x that satisfy the equation are algebraically equivalent. This means that we can scale a, b, and c by whatever we want, and if overflow or underflow is happening, we could choose a value to scale all three values to sidestep the issue.

Or that’s the strategy anyways. In practice, this can be quite tricky to implement so I’ll just refer you to this paper which documents the strategy, citing earlier work from Floating-point Computation by Sterbenz (1974).

As this issue is exceedingly rare, I have not actually found a case myself where I felt the need to handle overflow/underflow within the quadratic solver itself. Still though, the general trick to scale all inputs by some constant is the useful takeaway. Operating in a scaled or shifted space is a general strategy, something we do often in rendering (e.g. rendering view models or rendering large worlds).

Strategies for improving robustness

The quadratic formula is a great case study, not because I expect everyone is solving the quadratic equation often, but because it covers most of the points to think about when approaching floating point issues in general.

The main points to consider are:

• What are the edge cases where the results might change drastically depending on the inputs? These cases probably warrant dedicated handling.
• In cases where catastrophic cancellation could be happening, can I algebraically manipulate the expression to avoid it? If not, before trying to compute the expression at higher precision, does my code behave “well enough” in cases where the cancellation is catastophic to justify the hassle?
• Instead of trying to achieve maximal precision in the general case, can I scale or shift the inputs to the algorithm to stay in “healthy” floating point ranges where overflow, underflow, or even cancellation issues aren’t a problem?
• Does the computation I’m trying to do actually justify the use of double precision instead?

Fast math? More like reckless math!

The -ffast-math flag has got to be one of the worst named compiler flags in history. If I could rename it, I’d call it -freckless-math instead, so I’ll refer to it as such, henceforth.

By requesting reckless-math optimizations from the compiler, you are giving it a license to kill.

Ok maybe not that bad, but you are letting the compiler:

• …assume that +0 and -0 are interchangeable (aka that signed zeroes don’t exist).
• …assume that floating point operations can’t ever produce NaN values or infinities.
• …treat floating-point math as real numbers, meaning that addition and multiplication are assumed to be associative, commutative, and obey other algebraic identities.

By my count, none of these assumptions are correct, so reckless-math is zero for three.

• Signed zeroes are often quite useful because an algorithm might take the sign of an expression to use “in some way,” and the sign carries meaning even if the value in question is zero.
• Even if NaN values aren’t something we typically produce intentionally, plenty of numerical algorithms handle infinities perfectly well. In fact, many algorithms may even return infinities intentionally. In other cases, producing a NaN may be valid if the calling code expects NaN as a potential failure case output.
• Pretending floating-point numbers are real numbers when they are decidedly not is the worst offense of reckless-math.

As we saw in the last section, one of the most powerful tools we have at our disposal is algebra. In math-land, real numbers are in fact, real numbers. There are no precision limitations, no bits to saturate, and no issues with using algebraic identities. In other words, on paper, numbers behave like numbers!

Using reckless-math essentially negates the most powerful tool at our disposal. Any effort expended to carefully rearrange or transform an expression to exhibit desired properties with floating point values may be undone at the compiler’s discretion.

Needless to say, I think users really should think twice about using this flag, and I suspect that plenty of folks enable it because of branding. Who wouldn’t want their math to be fast? Maybe someone should run a campaign to rename the flag monkey-paw-fast-math instead!

Fast math in shaders

Shaders is an interesting case where fast-math is employed ubiquitously and its not actually easy to get away from it.

If you compile a shader with DXC for example, the default optimization level is -O3, and the moment the optimizer kicks in, -ffast-math is assumed. For plenty of shader calculations, admittedly, this is probably “fine” for some definition of fine.

Your options to disable it are to mark individual expressions/variables with the precise qualifier, or to disable fast math globally with -Gis. Unfortunately, neither of these options are very good.

The precise qualifier is actually a compiler bug factory since it’s quite challenging for the compiler to essentially reverse-propagate the annotation through the expression dependency graph.

On the other hand, the -Gis flag works by implicitly declaring all values as precise. This probably would be nice, except lots of driver unfortunately won’t emit fused-multiply-add instructions at all once expressions are marked precise (this is the case with many AMD cards I’ve tested for example).

Because FMA is such an important instruction for both accuracy and performance on modern GPUs, it’s not really something you can afford to lose. As such, the advice right now is just to let the compiler apply reckless-math everywhere and mark expressions/variables as precise only in cases where it’s needed (e.g. computing positions that need to match across shader programs to avoid z-fighting).

For the time being, I’ve filed a couple issues to DXC to hopefully have a more satisfactory resolution to this in the future:

• This issue is a feature request for an fma intrinsic that accepts single-precision floats that just demotes to a mad operation on hardware that doesn’t actually support FMA. The mad operation won’t emit an FMA for a lot of drivers if precise is used, and the existing fma intrinsic promotes arguments to double precision.
• This issue is a separate feature request that allows us to disable individual flags that control how the optimizer treats floating point. In an ideal world, we’d be able to, for example, disable FP associativity, no-signed-zeroes, or other reckless-math behaviors while still preserving operations we might want (FMA contraction being the chief optimization of interest).

Further reading

The Floating Point, Simply title unfortunately wasn’t meant to say that floating point was simple, or that there was some magical panacea to make all your floating point problems go away. Moreso that by foregoing reckless math and applying the techniques in the section on numerical robustness, I can avoid trouble most of the time.

Still though, there are plenty of exceptions to the rule where even more careful treatment of a computation is needed and the useful tricks aren’t sufficient. If what’s written above wasn’t enough, you’re probably in a space where you actually need to analyze sources of error propagation and turn to more sophisticated methods and/or existing literature to resolve your issues.

Instead of going through each case, I wanted to end on a bunch of writings I found useful over the years:

• Lecture Notes on the Status of IEEE Standard 754 for Binary Floating-Point Arithmetic
  • From the inventor of floats himself, William Kahan. You should totally read his writing if you haven’t already as it offers a great combination of technical precision and ranting in ideal proportions. I’ll suggest How Java’s Floating-Point Hurts Everyone Everywhere if you’re looking for a nice rant to unwind to some evening.
• What Every Computer Scientist Should Know About Floating-Point Arithmetic
  • This is a useful document to read to grok what everyone is talking about when they refer to ULPs, which sounds like a utterance of an alien, but is actually a way to define and systematically analyze error in floating point operations.
• Bruce Dawson on Floating Point
  • Bruce has a lot of posts on floating point computation that have helped me over the years. Unlike the other articles, his posts have more of a game dev slant, but also touch on idiosyncrasies of particular operating system or hardware quirks. Some of the material is less relevant in an x64-only world, but there’s still a lot of useful material here.
• Numerical Methods for Scientists and Engineers by Richard Hamming
  • This was one of the first books I’ve read that taught me how to “think about” numbers. It’s not exactly a light read, but the exposition is clear and has a great mixture of mathematical insight and practical considerations.

Discuss this post on Patreon.