>We fuzzed that function and found out that it does not give consistent results between Intel and AMD CPUs when SSE2 instructions are used.
I'm sad this wasn't investigated further. Was one of the implementations not standards compliant?
I'm sad this wasn't investigated further. Was one of the implementations not standards compliant?
In practice, some alternative implementations of these instructions treat the underspecification as license to do just enough to clear the threshold of being technically compliant with the standard and no more. This can show up as material differences in output in the least significant bits in some cases.
For floating point intensive algorithms, these small discrepancies and edge cases will occasionally bubble to the surface as material differences in high-level code behavior. My introduction to this class of bug many years ago was an application that was counting things for tax purposes based on geospatial relationships. On AMD, it produced a count that was off by one.
So far it looks like a CPU bug which produces NaN on some input.
The problem is most implementations include a few extra instructions or modes that are not defined in the IEEE spec. The Fast Reciprocal Estimate instruction that caused the issue in this example is one such instruction.
It's not defined in the IEEE spec. It's only defined in the x86 spec (and corresponding similar instructions are defined for other architectures). As its name implies, it's not defined to give an exact result. It's designed to be fast and it's designed to be a rough estimate to some level of precision.
You should not use such instructions if you care about any kind of accuracy or cross platform reproducibly. Stick to the instructions that are actually defined in the IEEE 754 spec. And stick to one of the recommended rounding modes, make sure all platforms use the same rounding mode.
This seems like an overly strong statement, as one reading of it is that "Fast Reciprocal Estimate" is just a bad rng function. I guess you mean that it's up to the user to provide their own error estimation?
One possibility is that there's an accidental equality test on the path of a runtime (AMD) value vs. a value that is computed at compile time (presumably on an Intel CPU).
The IEEE754 floating point representation gives you an easy way to roughly approximately convert between the log2 of it. The exponent gives you exactly the integer log2 of the number, and the mantissa gives you a fractional linear term that you can drop onto the integer part to make it closer to the actual log2 version. This lets you do some fun things fairly easily.
x -> log2(x) -> -log2(x) -> 1/x
x -> log2(x) -> -0.5 log2(x) -> 1/sqrt(x)
x -> log2(x) -> 0.5 log2(x) -> sqrt(x)
The lossiness of the conversion limits your precision, but there's a lot of times you don't give a shit. So the instructions are still valuable even if they're only approximately correct. For the reciprocal case, it will do something like this:
I'm not sure where the nondeterminism comes in between AMD and Intel. As you can see, with the reciprocal, there's no magic constant like there is with inverse square root. Maybe they're fudging something, or have a different form of Newton's method, I dunno.
https://pdf.sciencedirectassets.com/272990/1-s2.0-S157106610...
An initial thought I had that was that their table lookup is to find a magic constant that might nudge the result of the final value in the right direction. (instead of the magic constant 0x7F000000 that my code boils down to.) But that doesn't seem to be what my Kaby Lake is doing.
This is gonna bug me all night, I know it.
Technically they are. Practically many CPUs, especially older one, have non-trivial latency cost for passing a value between FP and integer ALUs, like couple of cycles each direction.
Ever wondered why there're 3 sets of bitwise instructions, e.g. pandn, andnps, andnpd which do exactly same thing with the bits in these registers? Now you know.
Over the complete range of floats, AMD is more precise on average, 0.000078 versus 0.000095 relative error. However, Intel has 0.000300 maximum relative error, AMD 0.000315.
Both are well within the spec. The documentation says “maximum relative error for this approximation is less than 1.5*2^-12”, in human language that would be 3.6621E-4.
Source code that compares them by creating 16GB binary files with the complete range of floats: https://gist.github.com/Const-me/a6d36f70a3a77de00c61cf4f6c1...
I’m not sure about that. I think it might be. Despite both are within specs, the results are numerically different.
Modern DirectXMath doesn‘t use these approximated instruction for inverting matrices, it’s open source, there: https://github.com/microsoft/DirectXMath/blob/83634c742a85d1...
True, but it's quite common that floating point ops throwaway bits (i.e. beyond the epsilon) will be numerically different between vendors.
> Modern DirectXMath doesn‘t use these approximated instruction for inverting matrices, it’s open source, there: https://github.com/microsoft/DirectXMath/blob/83634c742a85d1....
Good find, neither does Wine, actually: https://doxygen.reactos.org/dc/dd8/d3dx9math_8h.html#a3870c6...
Would be interesting to see if this bug also happens on Linux.
Divps instruction (32-bit float divide, the precise version) was relatively slow back then, e.g. on AMD K8 it had 18-30 cycles both latency and throughput, 21 cycles on AMD Jaguar, Core 2 Duo did in 6-18 cycles.
Fortunately, they fixed the CPUs. Skylake has 11 cycles latency of that instruction, Ryzen 10 cycles.
It’s similar for square root. Modern CPUs compute non-approximated 32-bit sqrt(x) in 9-12 cycles (Ryzen being the fastest at 9-10 but Skylake ain’t bad either at 12 cycles), old CPUs were spending like 20-40 cycles on that.
This is a common sentiment, and it is perhaps a helpful way to look at things if you don't want to dig into the details (even if it's not quite right).
But it's worth understanding the magnitude of errors. rcpps will be 3-4 orders of magnitude "more wrong" compared to the typical operation (if you view the "epsilon" of most operations to be the error after rounding). Or said another way: it would take the cumulative error from many thousands of adds and multiplies to produce the same error as one rcpps operation.
I think these implementation-specific shenanigans only applied to legacy x87 instructions which used 80-bit registers. SSE and AVX instructions operate on 32 or 64-bit floats. The math there, including rounding behavior, is well specified in these IEEE standards.
I’ve made a small test app: https://github.com/Const-me/SimdPrecisionTest/blob/master/ma...
It abuses AES instructions to generate long pseudo-random sequence of bits, then re-interprets these bits as floats, does some FP math on these numbers, and saves the output.
I’ve tested addition, multiplication, FMA, and float to integer conversions. All 4 output files, 1 GB / each, are bitwise equal between AMD desktop and Intel laptop.
Nice work! I'd be extremely curious to see if this still holds on intrinsics like `_mm512_exp2a23_round_ps` or `_mm512_rsqrt14_ps` (I'd wager it probably won't).
In this case (games), they probably made the right tradeoff despite the bug. But in general? I don't want to rehash the argument, I'm just really, really glad I'm not the one flipping that switch (or subject to the hordes of HNers who think I'm a monster for not flipping it).
https://bugs.llvm.org/show_bug.cgi?id=21385
On your graph, what’s red and what’s blue and what is on x axis?
X axis is source value, I've extracted a small subset of the data to graph in Excel, there: https://gist.github.com/Const-me/a6d36f70a3a77de00c61cf4f6c1...
Also, I've just recalled Y is absolute error not relative; that's why it follows shape of 1/x.
Thanks. If you still have the data (or if you can regenerate it) it would actually be possible to make a few small graphs that would cover the whole set:
Basically, the idea is to have x points as much as there are different exponents, which in 32-bit floats is at most 256. Then on the y axis one draws the number of bits of maximum distance between the really correct value of the mantissa and the calculated mantissa in the whole interval of one exponent. Those would allow comparing the implementations of Intel and AMD separately, and these graphs are what I'd be very interested to see. So the idea is to find the maximums in an interval, and there are then a limited number of the intervals. Only if such graphs match between AMD and Intel it would be interesting to compare inside of the intervals, but differences on that level would be these that I'd expect would be the obvious ones making the most problems, where the results like in the article (full black instead of shadow) wouldn't be surprising.
For that don't have to juggle with huge files, it's only 256 values per each CPU in that pass that are to be compared.
The results are not guaranteed to have same exponent, e.g. 1 / 0.499999 can be either 1.999999 or 2.0000001, both are correct within the precision, but have different exponents.
1) Use the binary representation of the numbers! To do so, cast the resulting float to the unsigned integer, then use bit masks and shifts to extract the exponent and mantissa. Note that the leading 1 is not explicit but implicit in the IEEE format unless it's a denormal number (so make it always explicit during the extraction).
2) use the exponent of the correct result as the "interval" reference.
If the exponents are the same, do a subtraction of the smaller from the bigger mantissa, that's the absolute "distance" between the two numbers -- the goal is to find what is the biggest absolute distance in which interval.
3) If one of the 2 values that are compared has different exponent, they can be converted to the same by a bit shift. Shift the mantissa of the one with the bigger exponent left accordingly. Again do the subtraction and use the result as the absolute distance. The goal is to figure out the biggest absolute distance in each interval (maintaining a maximum for each interval).
In short, think binary, not decimal, and measure using these values. Binary are only values that matter, decimal representation doesn't necessarily represent the exact values of bits.
Examples:
float 1.0 == unsigned 0x3f800000 here exponent is 127 == 2^0 and mantissa 0 with implicit 1 at the start i.e. explicit: 0x800000
float 0.999999940395355224609375 = unsigned 0x3f7fffff here exponent is 126 == 2^-1 and mantissa explicit: 0xffffff
The absolute distance between these two numbers is 1 (adding one to the lowest bit of mantissa of the smaller number would result the higher number 0xffffff + 1 = 0x1000000, the later is the mantissa adjusted to the same exponent of the smaller number (0x800000 << 1) ). If the "correct" number was 0x3f800000 and even if the shift was needed to calculate the absolute distance, the interval is still 0 (i.e. 0 is the x axis value, as its exponent was 2^0 i.e. 127, and the value to be plotted is on y is 1 until a bigger distance occurs).
For more examples of the format you can play here:
https://www.h-schmidt.net/FloatConverter/IEEE754.html
Also note that a few exponents are special, meaning infinity or NaN. Whenever the "correct" answer is not a NaN or infinity but the "incorrect" is, that should be treated specially, if it actually happens.
Total, exact, less, greater columns have total count of floats in a bucket. Sum of the “Total” gives 2^32, the total count of unique floats.
Computing max bit that’s different is too slow for the use case, neither SSE nor AVX have vector version of BSR instruction. Instead, I’m re-interpreting floats as integers and computing difference of the integers. maxLess, maxGreater, and maxAbs columns have that maximum error, measured as count of float values of the error. The value 4989 means the mantissa had like 12-13 lowest bits incorrect.
Source code is there: https://github.com/Const-me/SimdPrecisionTest/blob/master/rc... Not particularly readable because I’ve used AVX2 and OpenMP, however this way it takes less than a second on desktop, and maybe 1.5 seconds on a laptop to process all of these floats.
I’m guessing that the 3dnow! path was still being chosen based on a check like processor brand == AMD and processor generation >= x, assuming that 3dnow! would never be removed. This is something that has already been discovered being done in other games.