6 ms·

> The SIMD code does come with one asterisk, though: because floating-point addition is not associative, and it performs the summation in a different order, it

> The SIMD code does come with one asterisk, though: because floating-point addition is not associative, and it performs the summation in a different order, it may not get the same result as straight-line code. In retrospect, this is likely why the compiler doesn't generate SIMD instructions to compute the sum! > Does this matter for your use case? Only you can know! Anything is possible of course, and the typical way of writing a loop to sum up an array is literally telling the computer to go element by element and accumulate the values. But it seems pretty unlikely that doing, say, four accumulations in parallel with simd and then summing them up at the end is any more wrong than going element by element. Summing floats should by default be taken to have error bounds, and any answer in those bounds is valid. If you know something special about the floats you are inputting, the language should have some means to explicitly encode that. It shouldn’t be the most basic loop, that’s the default case, so it should give the best performance.

neonsunset 2y ago
A lot of code relies on FP operations being deterministic, often within confines of specific ISA. Applying SIMD to loops with FP could have been a default but it breaks a lot of existing code, makes the output frequently straight-up non-deterministic and as a result is something that a programmer has to explicitly choose to use. Moreover, many programmers may not be even aware of these, so when a `float Sum(float[] values)` starts to return something different, they may not have any way to know that it is the vectorization that makes it do so. Which is why, for example, .NET's standard library will use SIMD for `integers.Sum()` but not for `floats.Sum()`.
- crote 2y ago
  > A lot of code relies on FP operations being deterministic, often within confines of specific ISA How much of this is true in practice? I vaguely recall reading about a hard-to-diagnose bug where a seemingly unrelated code change meant that an intermediate value was no longer stored in the x87 extended-precision register but in memory instead, leading to a different result. Wouldn't you run into stuff like that all the time?
  - mindcandy 2y ago
    > A lot of code relies on FP operations being deterministic A lot, but still quite the minority. People do run into situations where the same code with the same data produces slightly different results all the time because of changes in compiler, compiler settings, cpu arch, etc… They just don’t notice. But, then they run into a situation where it matters and it’s a huge PITA to nail down.
    - kardos 2y ago
      Yes this is exactly it. The gold standard is bit-4-bit matching when shifting code between systems. If you move your code to a new system and get different results --> is it due to different summing order, or something else wrong? By ruling out the former you're more likely to get the bit-4-bit match and thus avoid a lot of debugging effort -- or spend it solving a real bug.
  - Sesse__ 2y ago
    This is one of many reasons why we generally don't use the x87 registers anymore.
throwaway260124 2y ago
Sorting the floats reduces the error. I think using multiple accumulators reduces the accuracy then. And sorted data is not uncommon. I think there is always a correct answer and the compiler shouldn’t make changes at least by default that are wrong. But ways for the programmer to express his intent more clearly are always welcome.
- Chabsff 2y ago
  Not necessarily. If the array is large and the values are similar, using accumulators would yield a more accurate result than a single one.
- celrod 2y ago
  Multiple accumulators increases accuracy. See pairwise summation, for example. SIMD sums are going to typically be much more accurate than a naive sum.
  - Chabsff 2y ago
    Not necessarily either. It's not particularly hard to create a vector where in-order addition is the most accurate way to sum its terms. All you need is a sequence where the next term is close to the sum of all prior ones. There just isn't a one-size-fit-all solution to be had here.
    - kardos 2y ago
      > There just isn't a one-size-fit-all solution to be had here. But there is: https://news.ycombinator.com/item?id=40867842 https://news.ycombinator.com/item?id=40867842
      - bee_rider 2y ago
        That’s a one size fits some solution. Not more than 2x as slower than native floats is pretty slow for cases where you don’t need the extra precision. If might be the case that it is the best solution if you do need the extra precision. But that’s just one case, not all.
pjc50 2y ago
> Summing floats should by default be taken to have error bounds, and any answer in those bounds is valid. If you know something special about the floats you are inputting, the language should have some means to explicitly encode that. It shouldn’t be the most basic loop, that’s the default case, so it should give the best performance. This is a lot of "should" for something that basically doesn't happen. The only information you give is the order of arithmetic in the original expression. It is a total nightmare if arithmetic is not stable between builds. Rebuilding the software and running it on the same input should not produce different results. (Long ago I encountered the special Intel version of this: because the FPU used 80-bit registers internally but 64-bit in memory, changing when a register fill/spill happened would change when your answer got rounded, and thereby change the results. You can set a global FPU flag at the start of the program to force rounding on every operation.)
- cwzwarich 2y ago
  > You can set a global FPU flag at the start of the program to force rounding on every operation This doesn’t do quite the same thing. It still uses the wider exponent range of the 80-bit type.
  - deleted 2y ago
    [deleted]
- bee_rider 2y ago
  It does happen in languages and libraries with a higher level of abstraction. MKL for example will do whatever it wants for accumulations (which practically means that it’ll take advantage of SIMD because that’s a big reason why people use the library) unless you specifically request otherwise via there “Conditional Numerical Reproducibility” flag. I think that was the right way to do it. BLAS made the right decision by defining these things in terms of sums and dot products instead of step-by-step instructions. It will always be possible to write programs that run differently on different hardware or with different optimization levels. If somebody is writing code for floating point computations and expects exact bit patterns—it is possible, of course, all the rounding behavior is clearly specified. But usually this is an error.
crote 2y ago
It turns into a serious problem if the floats are of different magnitude. Consider [1e50, -1e50, 1e3, 1e3]. Evaluating it as (((1e50 + -1e50) + 1e3) + 1e3) gives 2e3, evaluating it as ((1e50 + 1e3) + (-1e50 + 1e3)) gives 0. You get similar issues if you're trying to add a large number of small-magnitude numbers to a single larger-magnitude number. (((1e3 + 1e3) + 1e3) ... + 1e50) is substantially different from (((1e50 + 1e3) + 1e3) ... + 1e3).
- jandrese 2y ago
  Adding and subtracting numbers of wildly different magnitudes is always cursed. If you're doing this I'd say the majority of the time you are making an error. Your equation ends up being like "We calculated the diameter of the observable universe, but what if we added the Empire State Building to one side, how does that affect the result?"
  - skybrian 2y ago
    It’s also what we do when incrementing a counter, if it’s a very big counter. Operations like that don’t have a physical meaning, but can happen in things like random number generators where it’s a way of mixing the bits. Or maybe allocating unique ids, if you started from a random number.
    - bee_rider 2y ago
      Of course it is just an example so we shouldn’t nitpick you too much, but that sounds more like a job for an integer.
      - skybrian 2y ago
        Yes, absolutely. Though in JavaScript we often make do with the safe integer range. (I don’t know how fast BigInt is for 64 bit integers.)
        bee_rider 2y ago
RobotToaster 2y ago
> this is likely why the compiler doesn't generate SIMD instructions to compute the sum! Isn't this the type of thing ffast-math is supposed to enable?
- not2b 2y ago
  -ffast-math can be read as -fwrong-math if the accuracy if the answers matters. For some applications accuracy might not matter much. But the differences you get by pretending that floating point math is associative can be enormous, and scientific and engineering code is often carefully designed to take the likely scale of values into account, and rearranging it is not a good idea.
  - Dylan16807 2y ago
    Assuming associativity is not anti-accuracy in general, it just gets in the way of certain clever algorithms.
  - dzaima 2y ago
    For code that's been explicitly manually optimized already, indeed you wouldn't want to turn on -ffast-math; and I'd hope anyone who has actually learned what benefits from manual optimization and how to properly do that would also at some point have read something about not coupling that with -ffast-math. But I really doubt such meaningfully-manually-optimized code is at all frequent, and where it isn't it's a rather simple improvement (potentially with using some subset flags to allow NaN or infinity if desired).
  - rurban 2y ago
    Rearranging, like sorting before, is actually a good idea to minimize errors.
- Dylan16807 2y ago
  Fast math does a lot of risky things in addition to using basic math rules.
- bee_rider 2y ago
lscharen 2y ago
Surprisingly there are different algorithms for doing something as simple as summing up a list of numbers. The naïve way of adding numbers one-by-one in a loop is an obvious way, but there are more sophisticated methods that give better bounds of the total accumulated error; Kahan summation[1] being one of the better-known ones. Like most things, it can get complicated depending on the specific context. For streaming data, adding numbers one at a time may be the only option. But what if one could use a fixed-size buffer of N numbers? When a new number arrives what should be done? Take a partial sum of some subset of the N numbers in the buffer and add that to a cumulative total? If a subset if chosen, how? Are there provable (improved) error bounds for the subset selection method? Fun stuff. [1] https://en.wikipedia.org/wiki/Kahan_summation_algorithm https://en.wikipedia.org/wiki/Kahan_summation_algorithm
- kardos 2y ago
  As far as I know, xsum [1] more or less solves the problem completely: order invariant and exact, at less than 2x the cost of the naive summation. Further speeding it up may be the only avenue left for improvement. [1] https://gitlab.com/radfordneal/xsum https://gitlab.com/radfordneal/xsum
  - lscharen 2y ago
    I was not aware of this (2015) work -- very nice! A couple of pull-quotes from the paper to summarize: Much work has been done on trying to improve the accuracy of summation. Some methods aim to somewhat improve accuracy at little computational cost, but do not guarantee that the result is the correctly rounded exact sum. Many methods have been developed that instead compute the exact sum of a set of floating-point values, and then correctly round this exact sum to the closest floating-point value. This obviously would be preferable to any non-exact method, if the exact computation could be done sufficiently quickly Exact summation methods fall into two classes — those implemented using standard floating point arithmetic operations available in hardware on most current processors, such as the methods of Zhu and Hayes (2010), and those that instead perform the summation with integer arithmetic, using a “superaccumulator”. I present two new methods for exactly summing a set of floating-point numbers, and then correctly rounding to the nearest floating-point number. ... One method uses a “small” superaccumulator with sixty-seven 64-bit chunks, each with 32-bit overlap with the next chunk, allowing carry propagation to be done infrequently. The small superaccumulator is used alone when summing a small number of terms. For big summations, a “large” superaccumulator is used as well. It consists of 4096 64-bit chunks, one for every possible combination of exponent bits and sign bit, plus counts of when each chunk needs to be transferred to the small superaccumulator. On modern 64-bit processors, exactly summing a large array using this combination of large and small superaccumulators takes less than twice the time of simple, inexact, ordered summation, with a serial implementation

neonsunset 2y ago

A lot of code relies on FP operations being deterministic, often within confines of specific ISA. Applying SIMD to loops with FP could have been a default but it breaks a lot of existing code, makes the output frequently straight-up non-deterministic and as a result is something that a programmer has to explicitly choose to use. Moreover, many programmers may not be even aware of these, so when a `float Sum(float[] values)` starts to return something different, they may not have any way to know that it is the vectorization that makes it do so. Which is why, for example, .NET's standard library will use SIMD for `integers.Sum()` but not for `floats.Sum()`.

crote 2y ago
> A lot of code relies on FP operations being deterministic, often within confines of specific ISA How much of this is true in practice? I vaguely recall reading about a hard-to-diagnose bug where a seemingly unrelated code change meant that an intermediate value was no longer stored in the x87 extended-precision register but in memory instead, leading to a different result. Wouldn't you run into stuff like that all the time?
- mindcandy 2y ago
  > A lot of code relies on FP operations being deterministic A lot, but still quite the minority. People do run into situations where the same code with the same data produces slightly different results all the time because of changes in compiler, compiler settings, cpu arch, etc… They just don’t notice. But, then they run into a situation where it matters and it’s a huge PITA to nail down.
  - kardos 2y ago
    Yes this is exactly it. The gold standard is bit-4-bit matching when shifting code between systems. If you move your code to a new system and get different results --> is it due to different summing order, or something else wrong? By ruling out the former you're more likely to get the bit-4-bit match and thus avoid a lot of debugging effort -- or spend it solving a real bug.
- Sesse__ 2y ago
  This is one of many reasons why we generally don't use the x87 registers anymore.

throwaway260124 2y ago

Sorting the floats reduces the error. I think using multiple accumulators reduces the accuracy then. And sorted data is not uncommon. I think there is always a correct answer and the compiler shouldn’t make changes at least by default that are wrong. But ways for the programmer to express his intent more clearly are always welcome.

Chabsff 2y ago
Not necessarily. If the array is large and the values are similar, using accumulators would yield a more accurate result than a single one.
celrod 2y ago
Multiple accumulators increases accuracy. See pairwise summation, for example. SIMD sums are going to typically be much more accurate than a naive sum.
- Chabsff 2y ago
  Not necessarily either. It's not particularly hard to create a vector where in-order addition is the most accurate way to sum its terms. All you need is a sequence where the next term is close to the sum of all prior ones. There just isn't a one-size-fit-all solution to be had here.
  - kardos 2y ago
    > There just isn't a one-size-fit-all solution to be had here. But there is: https://news.ycombinator.com/item?id=40867842 https://news.ycombinator.com/item?id=40867842
    - bee_rider 2y ago
      That’s a one size fits some solution. Not more than 2x as slower than native floats is pretty slow for cases where you don’t need the extra precision. If might be the case that it is the best solution if you do need the extra precision. But that’s just one case, not all.

pjc50 2y ago

> Summing floats should by default be taken to have error bounds, and any answer in those bounds is valid. If you know something special about the floats you are inputting, the language should have some means to explicitly encode that. It shouldn’t be the most basic loop, that’s the default case, so it should give the best performance. This is a lot of "should" for something that basically doesn't happen. The only information you give is the order of arithmetic in the original expression. It is a total nightmare if arithmetic is not stable between builds. Rebuilding the software and running it on the same input should not produce different results. (Long ago I encountered the special Intel version of this: because the FPU used 80-bit registers internally but 64-bit in memory, changing when a register fill/spill happened would change when your answer got rounded, and thereby change the results. You can set a global FPU flag at the start of the program to force rounding on every operation.)

cwzwarich 2y ago
> You can set a global FPU flag at the start of the program to force rounding on every operation This doesn’t do quite the same thing. It still uses the wider exponent range of the 80-bit type.
- deleted 2y ago
  [deleted]
bee_rider 2y ago
It does happen in languages and libraries with a higher level of abstraction. MKL for example will do whatever it wants for accumulations (which practically means that it’ll take advantage of SIMD because that’s a big reason why people use the library) unless you specifically request otherwise via there “Conditional Numerical Reproducibility” flag. I think that was the right way to do it. BLAS made the right decision by defining these things in terms of sums and dot products instead of step-by-step instructions. It will always be possible to write programs that run differently on different hardware or with different optimization levels. If somebody is writing code for floating point computations and expects exact bit patterns—it is possible, of course, all the rounding behavior is clearly specified. But usually this is an error.

crote 2y ago

It turns into a serious problem if the floats are of different magnitude. Consider [1e50, -1e50, 1e3, 1e3]. Evaluating it as (((1e50 + -1e50) + 1e3) + 1e3) gives 2e3, evaluating it as ((1e50 + 1e3) + (-1e50 + 1e3)) gives 0. You get similar issues if you're trying to add a large number of small-magnitude numbers to a single larger-magnitude number. (((1e3 + 1e3) + 1e3) ... + 1e50) is substantially different from (((1e50 + 1e3) + 1e3) ... + 1e3).

jandrese 2y ago
Adding and subtracting numbers of wildly different magnitudes is always cursed. If you're doing this I'd say the majority of the time you are making an error. Your equation ends up being like "We calculated the diameter of the observable universe, but what if we added the Empire State Building to one side, how does that affect the result?"
- skybrian 2y ago
  It’s also what we do when incrementing a counter, if it’s a very big counter. Operations like that don’t have a physical meaning, but can happen in things like random number generators where it’s a way of mixing the bits. Or maybe allocating unique ids, if you started from a random number.
  - bee_rider 2y ago
    Of course it is just an example so we shouldn’t nitpick you too much, but that sounds more like a job for an integer.
    - skybrian 2y ago
      Yes, absolutely. Though in JavaScript we often make do with the safe integer range. (I don’t know how fast BigInt is for 64 bit integers.)
      - bee_rider 2y ago

RobotToaster 2y ago

> this is likely why the compiler doesn't generate SIMD instructions to compute the sum! Isn't this the type of thing ffast-math is supposed to enable?

not2b 2y ago
-ffast-math can be read as -fwrong-math if the accuracy if the answers matters. For some applications accuracy might not matter much. But the differences you get by pretending that floating point math is associative can be enormous, and scientific and engineering code is often carefully designed to take the likely scale of values into account, and rearranging it is not a good idea.
- Dylan16807 2y ago
  Assuming associativity is not anti-accuracy in general, it just gets in the way of certain clever algorithms.
- dzaima 2y ago
  For code that's been explicitly manually optimized already, indeed you wouldn't want to turn on -ffast-math; and I'd hope anyone who has actually learned what benefits from manual optimization and how to properly do that would also at some point have read something about not coupling that with -ffast-math. But I really doubt such meaningfully-manually-optimized code is at all frequent, and where it isn't it's a rather simple improvement (potentially with using some subset flags to allow NaN or infinity if desired).
- rurban 2y ago
  Rearranging, like sorting before, is actually a good idea to minimize errors.
Dylan16807 2y ago
Fast math does a lot of risky things in addition to using basic math rules.
bee_rider 2y ago

lscharen 2y ago

Surprisingly there are different algorithms for doing something as simple as summing up a list of numbers. The naïve way of adding numbers one-by-one in a loop is an obvious way, but there are more sophisticated methods that give better bounds of the total accumulated error; Kahan summation[1] being one of the better-known ones. Like most things, it can get complicated depending on the specific context. For streaming data, adding numbers one at a time may be the only option. But what if one could use a fixed-size buffer of N numbers? When a new number arrives what should be done? Take a partial sum of some subset of the N numbers in the buffer and add that to a cumulative total? If a subset if chosen, how? Are there provable (improved) error bounds for the subset selection method? Fun stuff. [1] https://en.wikipedia.org/wiki/Kahan_summation_algorithm https://en.wikipedia.org/wiki/Kahan_summation_algorithm

kardos 2y ago
As far as I know, xsum [1] more or less solves the problem completely: order invariant and exact, at less than 2x the cost of the naive summation. Further speeding it up may be the only avenue left for improvement. [1] https://gitlab.com/radfordneal/xsum https://gitlab.com/radfordneal/xsum
- lscharen 2y ago
  I was not aware of this (2015) work -- very nice! A couple of pull-quotes from the paper to summarize: Much work has been done on trying to improve the accuracy of summation. Some methods aim to somewhat improve accuracy at little computational cost, but do not guarantee that the result is the correctly rounded exact sum. Many methods have been developed that instead compute the exact sum of a set of floating-point values, and then correctly round this exact sum to the closest floating-point value. This obviously would be preferable to any non-exact method, if the exact computation could be done sufficiently quickly Exact summation methods fall into two classes — those implemented using standard floating point arithmetic operations available in hardware on most current processors, such as the methods of Zhu and Hayes (2010), and those that instead perform the summation with integer arithmetic, using a “superaccumulator”. I present two new methods for exactly summing a set of floating-point numbers, and then correctly rounding to the nearest floating-point number. ... One method uses a “small” superaccumulator with sixty-seven 64-bit chunks, each with 32-bit overlap with the next chunk, allowing carry propagation to be done infrequently. The small superaccumulator is used alone when summing a small number of terms. For big summations, a “large” superaccumulator is used as well. It consists of 4096 64-bit chunks, one for every possible combination of exponent bits and sign bit, plus counts of when each chunk needs to be transferred to the small superaccumulator. On modern 64-bit processors, exactly summing a large array using this combination of large and small superaccumulators takes less than twice the time of simple, inexact, ordered summation, with a serial implementation