Hacker News Reader

Comment by savant_penguin | Hacker News Reader

"(even the same computation on the same system can produce different results for floats)"

wait...what? Is this real?

Yes and no.

1. Yes -- On old x86 systems, the x87 registers were 80-bits and the bottom bits were undefined. You could theoretically have a 1-bit difference on some results depending on what the bottom 26 bits (that were undefined, because you had 64-bit floats most of the time, even though the machine did things 80-bits at a time).

2. No -- Modern x86 systems use SSE registers, which are 64-bits. There are a whole slew of configuration options, but if Windows / Linux does their job, you should have the same rounding-errors across your program (unless you manually change the rounding options, but that's your own fault if you go that route).

3. Kinda yes -- Floating point operations are NON-associative. (A+B)+C does NOT equal A+(B+C). The best example is Python:

>>> (1+1)+2.053 9007199254740994.0 >>> 1+(1+2.053) 9007199254740992.0

253 is the first value which starts to "round off" 1.0 in double-precision. So (1+2.053) == 2.053, because the value was rounded off.

Due to #3, even if you did everything correctly (ie: set the rounding flags so they were consistent), if your arithmetic happened in slightly different orders (very common in network-games, where different players may have their updates in slightly different orders), you end up with a 1-bit error between clients... which completely borks your simulation.

Because #3 is nonintuitive, many people think that you get different results when using floats. But its simply due to non-associativity.

In case you're like me, reading this comment and thinking "how the heck is 1 plus 1 plus 2-ish like 90 gazillion?", it's:

  >>> (1+1)+2.0**53
  9007199254740994.0
  >>> 1+(1+2.0**53)
  9007199254740992.0

The exponent operator (double star) got transformed into italics apparently.

Oh! Thank you! The comment makes a lot more sense now!

The newlines also helped, by the way, because I was wondering what putting the two numbers next to each other was meant to indicate.

> On old x86 systems, the x87 registers were 80-bits and the bottom bits were undefined. You could theoretically have a 1-bit difference on some results

I spent a long time once debugging an issue arising from this. Code was basically:

  double x = y;
  [...]
  if (x>y) fail;

Nothing in between the assignment and the check changed the value of x or y, yet the check triggered the fail condition. Turned out to be that one of the value stayed in the 80 bit x87 register and one was written to RAM as a 64 bit value, then loaded back into an x87 register for the check, resulting in the inequality. Running in a debugger wrote both values out of the registers before the check, making the problem unreproducible.

Nothing in between the assignment and the check changed the value of x or y

What was the check for?

There was code in between that could have changed the value, but it was never triggered in the bug case.

Ah yeah, that makes sense. Relying on FP equality should always make one's eye twitch a little but this is a neat way of hiding it.

Yes, that's a hard bug to catch, but if the mantra of "every floating point comparison must use a context-appropriate epsilon" is followed, then it wouldn't have happened in the first place.

Absolutely, and the fix involved adding a tolerance check. Not my code, I was the new guy at the time, and my job was trying to fix the bugs nobody else wanted to deal with...

I hate that mantra, as it produced more bugs and inaccuracies in the code I have to work with than necessary. Actually you should be able to rely on the equality of a value with its copy no matter what. We are past the FPU days, so I think we can let go of these mindless epsilons sprinkled over the code.

(Even during the FPU days, every ; should have rounded each value to the size of the variables, unless you use a special compiler switch to make FP math faster)

x87 has been functionally superseded (largely by various SSE iterations), but it is still supported for backwards-compatibility.

Modern applications and compilers largely do not use x87, but they can, and x87 results depending on hidden bits persists in modern x86_64 CPUs.

yet it is still used for arithmetic with "long double"

Wouldn't that resovle the issue? long double contains the whole 80 bit state, not much space for undefined bits left.

Yes, it resolves the issue in practice. However, the issue is an implementation detail of the compiler, not of the floating point format itself.

DLang allow to use it explicitly by using Real type, I instead of Double and Float types.

> Modern x86 systems use SSE registers

Are your sure of this? It's my understanding that SSE registers require the use of a special API and standard floating point operations do not use them.

But the last time I worked with them was writing a SIMD vector library 7 years ago.

Compilers now exclusively target the SSE registers (with scalar SSE instructions like mulss), and they did 7 years ago too.

https://gcc.gnu.org/onlinedocs/gcc/x86-Options.html

  -mfpmath=unit
  
      Generate floating-point arithmetic for selected unit unit. The choices for unit are:
  
      ‘387’
  
          Use the standard 387 floating-point coprocessor present on the majority of chips and emulated otherwise. Code compiled with this option runs almost everywhere. The temporary results are computed in 80-bit precision instead of the precision specified by the type, resulting in slightly different results compared to most of other chips. See -ffloat-store for more detailed description.
  
          This is the default choice for non-Darwin x86-32 targets.

I'm pretty certain all 64-bit x86 (so amd64) systems support SSE 2, so FP calcs are done with these registers.

> It's my understanding that SSE registers require the use of a special API and standard floating point operations do not use them.

No, SSE registers are used by most compilers that do floating point these days. They support all the usual IEEE float math operations and a host of bit twiddling operations as well as vector operations. The vector operations do still require using compiler intrinsics in C++, although some autovectorization does occur in gcc, icc, and llvm.

Ah yes. This is what I was remembering. For vector operations you have to use special intrinsics

64-bit x86 CPUs are defined to support SSE2 as a minimum, and all 64-bit operating systems use it by default.

For 32-bit systems, you still theoretically have to check, but I think every halfway modern compiler also uses SSE2 as the default target.

Wrt #3, I've had huge errors (I think on the order of 10e-3) in non linear curve fitting spectroscopy algorithms because of this. One of the physics research fellows in the group looked at me like I was an idiot for not knowing the order of multiplication mattered (I still have no clue why he would think it's a standard thing to know this).

Lots of researchers cargo cult floating point programming. I had a Fortran program a few months ago where the author did "if var > 0.99 and var < 1.01 then" etc., where "var" was an integer. I tried to search back where and how that was ever done, but no scenario made any sense ("var" was a categorical variable and always had been). So I went back to some of the original authors and there too they looked at me like I was an idiot and said "you should never test for equality, always check for distance within a certain epsilon". So I asked about the difference between integers and floating point and then they "didn't have time for technical details in code written 10 years ago." Shrug, they pay me to fix that particular sort of weirdness I guess.

In all likelyhood, the reason that's there isn't because a programmer was an idiot. Almost certainly, that variable was originally a float, then some time later another programmer came in and refactored it into an int, but missed fixing this condition, since it technically "works" for ints as well.

This is a good argument for stricter typechecking than anything else. In Rust, a condition like that wouldn't compile because of the mixed types. That's probably a good call in my opinion: it alerts the programmer that there's something funky going on, and it forces the programmer to be explicit about how the comparison should actually happen.

I do find myself doing the "cargo cult" thing from time to time. I was working on a codebase where Lua was embedded in C++, and there were all these flags that had to work in both systems. Because Lua only supports 64-bit floats as a numeric type, the natural type for these were doubles, and "by convention" they only ever had the values 0.0 or 1.0. It was a bit silly of me, but everytime i saw the line `if (fSomeFlag == 1.0f) ...` in the code, I would wince. It just feels wrong to do that to floats. I found myself regularly writing `if (fSomeFlag > 0.5f) ...` instead.

Of course, I was wrong about this. These weren't values determined as a result of arithmetic, they were always set to whole numbers. There was never going to be a case where the actual value of the flag would be 0.9999998f or whatever. But still: once you've been burned by the "floating point equality comparison", you live in fear of it forever after.

It was the same in the situation I was talking about; all numbers that would always be integers because of their nature (and in the data source they came from, they were integers, there was no way this could ever have been floating point, or that the values ever had any non-0 digits after the comma). And well, they probably weren't 'idiots' in the absolute sense, in fact I'm 100% sure that they are much smarter (and/or just plain better humans in other ways) than I am in many respects.

My point was that many researchers can only cobble together somewhat barely working software by taking snippets from their undergraduate textbooks (or googling/stack overflowing, for those under 40), and have no idea about many/most of the underlying principles. If only they would recognize that, and leave the software development to professionals, instead of treating it as an inconsequential implementation detail that is beneath them. Then again, if they did, they wouldn't have to pay me what they do to fix it up afterwards. So meh?

> . No -- Modern x86 systems use SSE registers, which are 64-bits. There are a whole slew of configuration options, but if Windows / Linux does their job, you should have the same rounding-errors across your program

Knock, knock. Who is there? Floating-point contraction.

Common Lisp over here like:

  * (+ (+ 1 1) (expt 2 256))
  115792089237316195423570985008687907853269984665640564039457584007913129639938

  * (+ 1 (+ 1 (expt 2 256)))
  115792089237316195423570985008687907853269984665640564039457584007913129639938

Many languages support BigFloats and BigInts, but this discussion is about IEEE 754 floats vs Gustavson's proposals.

Python will happily do this too, if you remove the ".0" in the example, the values become integers, instead of floats, and the computation is precise.

The point is to demonstrate a property of IEEE floats, not the language. I'd expect CL has some means of using IEEE floats too, and with the appropriate syntax, the same issues could be demonstrated.

Not really. I suspect they are referring to the original x87 FP functions with their odd 80 bit FP stack and operations, where the exact result depends on when intermediate results are written to RAM and when they are kept in the FP stack, which might depend on intransparent things like compiler optimisations.

They could refer to parallel systems, most likely GPUs, where it is common to have slightly different results even on the same operation, on the same system.

This is due to the fact that several pipelines will do calculations in parallel and depending on which pipeline ends first (which can depend on many factors like current temperature) the sums can happen in a different order, leading to different rounding results.

I've never heard of that happening on a GPU. Can you give an example?

You could have this kind of behavior if you use atomic functions - but you need to explicitly opt-in actually using those and it is not limited to GPU.

And performance will suffer as a consequence.

We have this in our GPU computations, this is very annoying. The errors are usually too small to matter, but they make comparing binary results / reproducible builds impossible. Here is an example I just made:

    $ python3
    >>> s = [10**-x for x in range(16)]
    >>> s
    [1, 0.1, 0.01, 0.001, 0.0001, 1e-05, 1e-06, 1e-07, 1e-08, 1e-09, 1e-10, 1e-11, 1e-12, 1e-13, 1e-14, 1e-15]
    >>> sum(s)
    1.1111111111111112
    >>> sum(reversed(s))
    1.111111111111111

(it uses CPU for simplicity, but GPU has exactly the same math)

Summation of an array is surprisingly non-trivial. Below a discussion about it with many interesting links (including Kahan summation and pairwise summation) if you want to go down some unexpected rabbit holes.

Somewhere I've seen a little routine (by Stefan Karpinski of Julia fame) that starts with a fixed pre-determined array of numbers, and you give it a desired target number, and it permutes the numbers of the array such that when you sum them up (naively, from left to right) the result is your desired target number [2].

https://discourse.julialang.org/t/accurate-summation-algorit...

EDIT: [2] https://discourse.julialang.org/t/array-ordering-and-naive-s...

If you don't care about cancellation error and are only concerned with perfect bitwise / reproducible results, then sort your numbers.

Sorting the numbers is necessary to minimize cancellation error. But absolute minimum cancellation error can only be done sequentially: the parallel version may return different results.

    l = list(s)
    l.sort()
    sum(l)
    l = list(reversed(s))
    l.sort()
    sum(l)

Both return 1.111111111111111 now. Sorting is the "synchronization" step you need to get all your parallel tasks back into a defined order.

> They could refer to parallel systems, most likely GPUs, where it is common to have slightly different results even on the same operation.

Parallel systems have great difficulty in defining an order of operations. Floating-point numbers are non-associative, so order matters for bitwise compatibility.

With that being said, you CAN define an order, and it CAN be done in parallel. The parallel-prefix sum is a fully defined order of operations for example, and if you sort all numbers (from smallest to largest), you'll get a more accurate result.

https://en.wikipedia.org/wiki/Prefix_sum#Algorithm_1:_Shorte...

The + operations always happen in the same order, even on a GPU. So you should get the same results every time with floating-point math. But it is non-intuitive for many people to work with non-associative floating-point math.

-----------

Note that writing a sequential algorithm with fully defined order is still grossly difficult! IMO, if anyone wants a bit-accurate simulation, use 64-bit integers instead.

If you absolutely must use floats, then define a proper order of all operations, and include sorting the numbers (!!) to force a fully defined order as much as possible. You may need to sort your pool of numbers many, many times per calculation. But that's what you have to do to get a defined ordering.

it doesn't happen much in practice, but the official IEEE spec is silent about it, so you can toss it into whatever bin contains 'undefined' behaviour of C.

Are any C undefined behaviors explicitly undefined spec?

>> where the exact result depends on when intermediate results are written to RAM and when they are kept in the FP stack, which might depend on intransparent things like compiler optimisations.

One thing I've never really been sure about is how context-switching affects this. I've always assumed this could also be a source of non-deterministic floating-point results: when the OS decides to switch threads, it has to save the CPU & FPU state, and restore it when switching back. Since context switches are (practically speaking) non-deterministic on a pre-emptive OS, this would imply truncating 80-bit to 64-bit floating point values could also happen on non-deterministic points in time.

Does anyone know if this is actually true? Or do the full 80-bits get saved/restored on context switches?

They are saved.

The FPU state is saved using the F(N)SAVE x86 instruction, which stores 108 bytes to memory, which contains all visible FPU state.

this is not allowed by the standard. it is false in practice too

Reply on news.ycombinator.com