12 ms·

"(even the same computation on the same system can produce different results for floats)" wait...what? Is this real?

fhars 7y ago
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.
- Iv 7y ago
  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.
  - shaklee3 7y ago
    I've never heard of that happening on a GPU. Can you give an example?
    - llukas 7y ago
      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.
      - jononor 7y ago
        And performance will suffer as a consequence.
    - theamk 7y ago
      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)
      - FabHK 7y ago
        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-algorithm/7163 https://discourse.julialang.org/t/accurate-summation-algorit... EDIT: [2] https://discourse.julialang.org/t/array-ordering-and-naive-summation/1929 https://discourse.julialang.org/t/array-ordering-and-naive-s...
      - dragontamer 7y ago
        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.
  - dragontamer 7y ago
    > 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:_Shorter_span,_more_parallel 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.
- dnautics 7y ago
  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.
  - zeckalpha 7y ago
    Are any C undefined behaviors explicitly undefined spec?
- w0utert 7y ago
  >> 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?
  - Tuna-Fish 7y ago
    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.
dragontamer 7y ago
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.
- deathanatos 7y ago
  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.
  - pure-awesome 7y ago
    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.
- martinpw 7y ago
  > 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.
xhgdvjky 7y ago
this is not allowed by the standard. it is false in practice too
- myaH 7y ago
  The correct comment at the bottom again, while everyone else misinterprets IEEE.

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.

Iv 7y ago
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.
- shaklee3 7y ago
  I've never heard of that happening on a GPU. Can you give an example?
  - llukas 7y ago
    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.
    - jononor 7y ago
      And performance will suffer as a consequence.
  - theamk 7y ago
    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)
    - FabHK 7y ago
      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-algorithm/7163 https://discourse.julialang.org/t/accurate-summation-algorit... EDIT: [2] https://discourse.julialang.org/t/array-ordering-and-naive-summation/1929 https://discourse.julialang.org/t/array-ordering-and-naive-s...
    - dragontamer 7y ago
      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.
- dragontamer 7y ago
  > 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:_Shorter_span,_more_parallel 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.
dnautics 7y ago
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.
- zeckalpha 7y ago
  Are any C undefined behaviors explicitly undefined spec?
w0utert 7y ago
>> 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?
- Tuna-Fish 7y ago
  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.

dragontamer 7y ago

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.

deathanatos 7y ago
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.
- pure-awesome 7y ago
  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.
martinpw 7y ago
> 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.

xhgdvjky 7y ago

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

myaH 7y ago
The correct comment at the bottom again, while everyone else misinterprets IEEE.