16

Consider this piece of C code:

#include <math.h>
#include <stdbool.h>
#include <stdio.h>

bool foo(int a, int b, int c, int d) {
    double P = atan2(a, b);
    double Q = atan2(c, d);
    return P < Q;
}

bool bar(int a, int b, int c, int d) {
    return atan2(a, b) < atan2(c, d);
}

int main() {
    if (foo(2, 1, 2, 1)) puts("true"); else puts("false");
    if (bar(2, 1, 2, 1)) puts("true"); else puts("false");
    return 0;
}

When I compile it with gcc -lm -m64, it prints

false
false

as expected. However, when I compile it with gcc -lm -m32, the output is

false
true

Why does such a trivial difference as storing the result in a variable change the behavior of the code? Is it because the floating point registers are 80-bit? How should I write my code such that weirdness like this is avoided?

(I originally used the variant in bar in a C++ program as a compare function for std::sort. This triggered undefined behavior because comp(a, a) returned true - how should I compare vectors by their angle without UB?)

Comment summary so far

It seems there is consensus that this is because of the differences in where the results are rounded to the range of the double. The C/C++ standards allow this, even when the results differ as a consequence.

A workaround of comparing the input as a * d < b * c, which would indeed help in this example. However, I seek to avoid this category of issues in the future, and not just fix the bug and move on. The issue could still manifest if the inputs themselves were doubles.

I accept that the comparison might be inaccurate due to the nature of floating-point arithmetic, however if the angles are this close, I still want to achieve the strict weak ordering required by std::sort, such that no UB happens.

dbush
  • 205,898
  • 23
  • 218
  • 273
Maya
  • 1,490
  • 12
  • 24
  • 3
    [What Every Computer Scientist Should Know About Floating-Point Arithmetic](https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html) – Jesper Juhl Jul 13 '20 at 16:40
  • @JesperJuhl this is definitely a useful reference. Are you referring to any specific section of the paper? – Maya Jul 13 '20 at 16:42
  • 8
    @JesperJuhl not related... – Sebastian Hoffmann Jul 13 '20 at 16:42
  • 3
    @NathanOliver I understand that the results won't be exact, but are they really going to be different between invocations of the same function, with the same parameters? – Maya Jul 13 '20 at 16:43
  • 1
    Funny thing is, when I build the 32-bit variant with optimizations enabled (`-O3`) then it printf `true` for both functions, There's definitely something fishy going on here. – Some programmer dude Jul 13 '20 at 16:44
  • In the `bar()` function, is the compiler allowed to replace the comparison with `(double)c * b < (double)a * d`? – Weather Vane Jul 13 '20 at 16:44
  • 2
    @NieDzejkob Yes. With `atan2(a, b) < atan2(c, d)` both results could be in a larger register, which makes the comparison more precise. With `P < Q;` you "truncate" down to the lower precision so that can change the result. – NathanOliver Jul 13 '20 at 16:45
  • @NieDzejkob Not refering to any specific section. Just thought that it might help you in general. Floating point is not exact and that paper contains lots of insights that are useful to know. Just thought I'd throw that in your direction as it might be helpful, that's all. – Jesper Juhl Jul 13 '20 at 16:46
  • `atan2` exists for both `float` and `double`, maybe this introduces the weirdness. The `foo` snippet will force the double overload whereas the `bar` function might take the float one (especially on 32bit) – Sebastian Hoffmann Jul 13 '20 at 16:46
  • godbolt link: https://godbolt.org/z/b44TPW – Sebastian Hoffmann Jul 13 '20 at 16:47
  • 5
    @SebastianHoffmann this is C, `atan2` for `float`s is called `atan2f`. I specifically switched away from C++ when I was debugging this to rule this out. – Maya Jul 13 '20 at 16:47
  • 4
    Why is this tagged as `C++` if the code is `C`? – PaulMcKenzie Jul 13 '20 at 16:48
  • @PaulMcKenzie as mentioned in the question, the motivation for this code is C++'s std::sort, and the issue occurs in both languages. – Maya Jul 13 '20 at 16:49
  • 8
    Seeing as the only difference in the asm is storing and loading the doubles P and Q (the latter only stores and loads one of them), it seems that this must have to do with the precision difference between the registers and an actual `double`. https://godbolt.org/z/fd79vG – Justin Jul 13 '20 at 16:49
  • C and C++ are two different languages. Whatever answer you get from someone running this as a C++ program may not have any relevance to the C program. – PaulMcKenzie Jul 13 '20 at 16:51
  • @PaulMcKenzie In this case, the way they handle the floating point numbers is identical. – Justin Jul 13 '20 at 16:51
  • @Justin Good to know. Still, how should I go about fixing this, to get a reliable compare function for use in std::sort? – Maya Jul 13 '20 at 16:51
  • 1
    Try `-m32 -mfpmath=sse`. If that generates machine instructions that can't be used in your 32-bit environment, try `-m32 -ffloat-store`. – zwol Jul 13 '20 at 16:52
  • @NieDzejkob No idea... – Justin Jul 13 '20 at 16:52
  • Note that the difference between 32 and 64 bit version is that the former uses FPU instructions/80-bit registers, and the latter uses SIMD/SSE vectorization registers with 64-bit width. – Daniel Langr Jul 13 '20 at 16:52
  • 1
    @NieDzejkob a reliable compare? If you don't need the actual angles, consider a cross-multiply as earlier comment. The integer arithmetic will be *exact* (subject to overflow). – Weather Vane Jul 13 '20 at 16:53
  • 1
    @NieDzejkob -- Don't use floating point calculations for things such as keys, loop variables, etc. Come up with an integer-based solution for these things. – PaulMcKenzie Jul 13 '20 at 16:55
  • @WeatherVane Okay, cross-multiplying like this is a good solution for this particular case, but what if the inputs were doubles in the first place? Could the same truncation issue show up? – Maya Jul 13 '20 at 16:55
  • 2
    Of course they could. – Weather Vane Jul 13 '20 at 16:56
  • @WeatherVane Thank you. What should I do then, to sort an array of 2D vectors by their angle, if the coordinates are floating-point? – Maya Jul 13 '20 at 16:57
  • 1
    Is this a gcc bug maybe? the difference between those two assemblies look like optimizations. is gcc allowed to change behaviour here? – Sebastian Hoffmann Jul 13 '20 at 16:58
  • Possibly relevant: https://stackoverflow.com/q/20869904/1896169 – Justin Jul 13 '20 at 16:59
  • 1
    @SebastianHoffmann I believe the standard allows for extra precision on floating point, even in this case where only one of the floating point numbers has the extra precision. In a case such as that, GCC's behavior (which is enabled even on -O0) would be legal. – Justin Jul 13 '20 at 17:00
  • @Justin But I don't think it allows different behaviors for different optimization settings - unless there is a UB in the code, which is not the case here. – Eugene Sh. Jul 13 '20 at 17:02
  • Note there is an inherent problem in what you consider "greater" since angles are cyclic. `atan2` returns a value in the range -π to π radians. – Weather Vane Jul 13 '20 at 17:06
  • 1
    It could be related to the precision rounding mode of FPU. On X86 it is handled in a non strictly compliant mode with respect to IEEE754, while SSE are compliant. See https://stackoverflow.com/questions/7295861/enabling-strict-floating-point-mode-in-gcc and https://www.linuxtopia.org/online_books/an_introduction_to_gcc/gccintro_70.html – Frankie_C Jul 13 '20 at 17:07
  • @WeatherVane I am aware. I just need to make sure that any "cyclic slice" of the sorted result is a "continuous" subset of the vectors. – Maya Jul 13 '20 at 17:09
  • 1
    @EugeneSh. That's obviously false. Most successful optimizations change the wall time a program uses and a program can modify its behavior by measuring the wall time it uses. There's no UB in this code, it's just that its output is not guaranteed to be any particular thing because the standard doesn't precisely specify what has to happen in this case. – David Schwartz Jul 13 '20 at 17:10
  • @NieDzejkob You're trying to sort based on a parameter that is specified with insufficient precision to support the way you're trying to use it. You can sort by angle, but it's not going to be stable or predictable, only approximate to the accuracy of the guaranteed resolution of the types you are using. – David Schwartz Jul 13 '20 at 17:13
  • @EugeneSh. Generally, it does. For instance, with high-level optimizations, division can be substituted by multiplication if the denominator is known at compile time. Then, the result will be different due to rounding, and, consequently, you get possibly different behavior for different optimization levels. – Daniel Langr Jul 13 '20 at 17:13
  • 2
    @DavidSchwartz If we're being pedantic: A program can't measure the time it takes to execute. It can *ask the system* to please tell it how long it's been executing. I.e. the time it takes for a program to execute is actually an *input* to the program. If you actually set up a system that provides the exact same input to an optimized program (including lying about its execution time) as it does to a non-optimized program then they're supposed to behave the same. We're not taking any input here anyway; the optimized and non-optimized programs *are* supposed to behave identically. – HTNW Jul 13 '20 at 17:16
  • @DanielLangr That's not allowed (except if the divisor is a power of two) for any of the default optimization levels up to `O3`. You have to explicitly tell GCC "Yes, please break my math" with `-funsafe-math-optimizations`/`-ffast-math`/`-Ofast`, which is your own fault if you use it. – HTNW Jul 13 '20 at 17:19
  • @HTNW As I said in the second part of my comment, they're not supposed to behave identically because the code relies on behavior not precisely specified. What precision you get in the result of calling `atan2` is, by your criteria, an input to the program. (You're still wrong about optimization though. For example, no rule prohibits an optimization level from changing the size of an `int`, which code can certainly detect.) – David Schwartz Jul 13 '20 at 17:20
  • @DavidSchwartz I do understand your point, but how does http://port70.net/~nsz/c/c11/n1570.html#5.1.2.3p6 apply here then? *"At program termination, all data written into files shall be identical to the result that execution of the program according to the abstract semantics would have produced."* - is the optimization changing the semantics of `atan2` function? – Eugene Sh. Jul 13 '20 at 17:23
  • 1
    You could always define a struct/class with x, y, and angle. Compute the angle once for each point. Then sort based on the precomputed angles. (Bonus: saves time versus computing atan2 twice for each comparison.) – user3386109 Jul 13 '20 at 17:25
  • 3
    @DavidSchwartz No, the behavior *is* precisely specified. Floating point math *is exact* in the sense that there is always *one right answer*, and the right answer here is that both functions give `false`. The wonkiness here is because GCC *gives up on implementing the specification* under `-m32`, as shown by the link given by @Frankie_C and also [this one](https://gcc.gnu.org/wiki/FloatingPointMath#Note_on_x86_and_m68080_floating-point_math). You're right that sorting by the result of `atan2` is wrong, though. It assigns an infinite triangular "slice" of the plane to each return value. – HTNW Jul 13 '20 at 17:25
  • 1
    @user3386109 Thank you. I think this is a good, generalizable solution. – Maya Jul 13 '20 at 18:01
  • 1
    @EugeneSh.: GCC with `-O0` is a different C implementation than GCC with `-O3`, and GCC with `-m32` is a different C implementation than GCC with `-m64`. So it does not matter that the behavior “changes” when optimization is turned on or other switches are changed. – Eric Postpischil Jul 13 '20 at 18:43
  • 1
    @HTNW: Re “No, the behavior is precisely specified”: Even if “floating point math” is exactly specified (it is actually not; IEEE-754 has some wiggle room, such as when tininess is detected), bindings of IEEE-754 (or other floating-point specifications) to C are not. – Eric Postpischil Jul 13 '20 at 18:44
  • Post result of `printf("%d\n", FLT_EVAL_METHOD);` in those 2 cases for more insight. – chux - Reinstate Monica Jul 13 '20 at 23:34
  • "This triggered undefined behavior because comp(a, a) returned true" --> not UB. – chux - Reinstate Monica Jul 13 '20 at 23:38

0 Answers0