26 comments

[ 2.7 ms ] story [ 43.2 ms ] thread
I've never heard anyone say any of these supposed myths, except for the first one, sort of, but nobody means what the first one pretends it means, so this whole post feels like a big strawman to me.
Same, but I still learned a fair amount anyways. I say no harm, no foul.
i am laughing and crying at the same time...

a couple of decades ago, I was still a lowly applied mathematics (and aspiring programmer) student when i proposed to a professor to change the compiler options from O2 to O3 to gain performance (reducing run time by a significant percentage like 50% or more...) only to be disregarded by saying "no, compiler optimizations break numerical accuracy".

I like these. They are push back against the sort of… first correction that people make when encountering floating point weirdness. That is, the first mistake we make is to treat floats as reals, and then we observe some odd rounding behavior. The second mistake we make is to treat the rounding events as random. A nice thing about IEEE floats is that the rounding behavior is well defined.

Often it doesn't matter, like you ask for a gemm and you get whatever order of operations blas, AVX-whatever, and OpenMP conspire to give you, so it is more-or-less random.

But if it does matter, the ability to define it is there.

My favorite thing about floating point numbers: you can divide by zero. The result of x/0.0 is +/- inf (or NaN if x is zero). There's a helpful table in "weird floats" [0] that covers all the cases for division and a bunch of other arithmetic instructions.

This is especially useful when writing branchless or SIMD code. Adding a branch for checking against zero can have bad performance implications and it isn't even necessary in many cases.

Especially in graphics code I often see a zero check before a division and a fallback for the zero case. This often practically means "wait until numerical precision artifacts arise and then do something else". Often you could just choose the better of the two options you have instead of checking for zero.

Case in point: choosing the axes of your shadow map projection matrix. You have two options (world x axis or z axis), choose the better one (larger angle with viewing direction). Don't wait until the division goes to inf and then fall back to the other.

[0] https://www.cs.uaf.edu/2011/fall/cs301/lecture/11_09_weird_f...

This is also my favorite thing about floating point numbers. Unfortunately languages like Python try to be smart and prevent me from doing it. Compare:

    >>> 1.0/0.0
    ZeroDivisionError
    >>> np.float64(1)/np.float64(0)
    inf
I'm so used to writing such zero division in other languages like C/C++ that this Python quirk still trips me up.
> you can divide by zero

It is implementation-dependent. It is not obligatory for implementation to respect IEEE 754.

You can exploit the exactness of (specific) floating point operations in test data by using sums of powers of 2. Polynomials with such coefficients produce exact results so long as the overall powers are within ~53 powers of 2 (don't quote me exactly on that, I generally don't push the range very high!). You can find exact polynomial solutions to linear PDEs with such powers using high enough order finite difference methods for example.

However, the story about non-determinism is no myth. The intel processors have a separate math coprocessor that supports 80bit floats (https://en.wikipedia.org/wiki/Extended_precision#x86_extende...). Moving a float from a register in this coprocessor to memory truncates the float. Repeated math can be done inside this coprocessor to achieve higher precision so hot loops generally don't move floats outside of these registers. Non-determinism occurs in programs running on intel with floats when threads are interrupted and the math coprocessor flushed. The non-determinism isn't intrinsic to the floating point arithmetic but to the non-determinism of when this truncation may occur. This is more relevant for fields where chaotic dynamics occur. So the same program with the same inputs can produce different results.

NaN is an error. If you take the square root of a negative number you get a NaN. This is just a type error, use complex numbers to overcome this one. But then you get 0. / 0. and that's a NaN or Inf - Inf and a whole slew of other things that produce out of bounds results. Whether it is expected or not is another story, but it does mean that you are unable to represent the value with a float and that is a type error.

NaN is not necessarily an error. It might be fine. It depends on what you're doing with it.

If NaN is invalid input for the next step, then sure why not treat it as an error? But that's a design decision not an imperative that everybody must follow. (I picture Mel Brooks' 15 commandments, #11-15, that fell and broke. This is not like that.)

When writing code, it's useful to follow simpler rules than strictly needed to make it easier to read for other coders (including you in the future). For example, you might add some unnecessary brackets to expressions to avoid remembering all the operator precedence rules.

The article is right that sometimes you can can use simple equality on floating point values. But if you have a (almost) blanket rule to use fuzzy comparison, or (even better) just avoid any sort if equality-like comparison altogether, then your code might be simpler to understand and you might be less likely to make a mistake about when you really can safely do that.

It's still sensible to understand that it's possible though. Just that you should try to avoid writing tricky code that depends on it.

Every time I use exact floating point equality I simply write a comment explaining why I did this. That's what comments are for.
Note: the improved loop in the "1. They are not exact" can easily hang. If count > 2^24, then the ulp of f is 2, and adding 1.0f leaves f unchanged. What's wild is that a few lines later he notes how above 2^24 numbers start “jumping” every 2. Ok, but FP are always "jumping", by their ulp, regardless of how big or small they are.
People need to understand rounding better, especially the topic of when rounding can happen and when it can't for the basic operations.

Updating with a concrete example: The Fortran standard defines the real MOD and MODULO intrinsic functions as being equivalent to the straightforward sequence of a division, conversion to integer, multiplication, and subtraction. This formula can round, obviously. But MOD can be (and thus should be) implemented exactly by other means, and most Fortran compilers do so instead. This leaves Fortran implementors in a bit of a pickle -- conform to the standard, or produce good results?

Definitely. People think they'll get out of knowing how rounding works by using arbitrary precision arithmetic, but arguably it's even more important there (you run out of precision/memory at some point; what do you think happens then?). You can use floats for money if you do the rounding right.
The core thing to know about floating point numbers in comp langs is they aren't floating point numbers.

They are approximations of floating point numbers, and even if your current approximation of a floating point number to represent a value seems to be accurate as an integer vale...

There is no guarantee if you take two of those floating point number approximations that appear completely accurate that the resulting operation between them will also be completely accurate.

> They are not exact

It's not exactly a myth, as the article mentions, they're only exact for certain ranges of values.

> NaN and INF are indication of an error

This is somewhat semantic, but dividing by zero typically does create a hardware exception. However, it's all handled behind the scenes, and you get "Inf" as the result.

You can make it so dividing by zero is explicitly an error, see the "ftrapping-math" flag.

Here's one that's not a myth: IEEE-754 floats are the only "primitive types" that allow "a == a" to not be true.

I.e., two floats that are _identical_ to each other (even when it's _the same_ variable, on the same memory address) can be not _equal_ to each other, specifically if it's NaN. This is dictated by IEEE-754, and this is true for all programming languages I know, and to this day, this makes zero sense to me, but apparently this is useful for some reason.

>>> A beginner would use them, trusting they are infinitely capable and precise, which can lead to problems. An intermediate programmer knows that they have some limitations, and so by using some good practices the problems can be avoided.

Amusingly, pocket calculators became affordable while I was in high school, and anybody who was interested in math learned the foibles of floating point almost instantly. Now it's an advanced topic. A difference may have been that our calculators had very few significant digits, so the errors introduced by successive calculations became noticeable more quickly. Also, you had to think about what you were doing because, even if you trusted the calculator, you didn't trust your fingers.

In the context of double precision the article says

> the largest integer value that can be represented exactly is 2^53

— I am confused as to why it not 2^52, given that there are 52 bits of mantissa, so relative accuracy is 2^-52, which translates to absolute accuracy larger than 1 after 2^52. Compare this to the table there saying "Next value after 1 = 1 + 2^-52".