5 ms·

I met William Kahan about 15 years ago and asked him about this. His response: "If your program says double then yes, of course it's a bug if the compiler does

I met William Kahan about 15 years ago and asked him about this. His response: "If your program says double then yes, of course it's a bug if the compiler does long-double precision arithmetic instead. Isn't that obvious? Why are you even asking this?" If the author of the IEEE 754 standard says that the way you do floating-point arithmetic is obviously a bug, you're probably doing it wrong.

acqq 10y ago
And I believe he's right, if he was asked specifically what that program should do. If all the statements in the program are really just: #include <stdio.h> void test(double x, double y) { const double y2 = x + 1.0; if (y != y2) printf("error\n"); } void main() { const double x = .012; const double y = x + 1.0; test(x, y); } Then the comparison has to be done on doubles, and "but it's extended precision internally" mentioned in some answers is the wrong answer. All the variables in the program are stated to be doubles. If I'd write the x86 and x87 asm instructions I could easily write them so that they work as expected, even without SSE2. However if you compile with some "special" optimizations that are documented to relax some specific code generation rules, then it can be OK that these specific optimizations behave differently, but then: you asked for it by choosing that documented optimization. In the answers there, Richard Henderson wrote 2001-01-16 06:30:04 UTC: "You are seeing the results of excess precision in the FPU. Either change the rounding precision in the FPCR, or work around the problem with -ffloat-store." He's also right, if the -ffloat-store is documented to be turned off by -O and if it's documented that that flag also affects the precomputed calculations that are done in the compile time. If the second detail is not documented, the results could be unexpected, but again depends of what you expect from the optimizer. The properly written compiler and optimizer should perform the calculations during the code generation correct to the definition and not cut the corners. If the generated code was generated to actually perform the != during the run time, then the compile-time calculations are irrelevant, but the sane default should be something like -ffloat-store when not using SSE2 even with -O, if -O has the semantic "the optimizations that should not produce really unexpected effects." It sounds complicated, but it has its reasons. P.S. answer to nsajko's post below: I mention -ffloat-store not ffast-math. Read the description of ffloat-store please: https://gcc.gnu.org/onlinedocs/gcc-3.2/gcc/Optimize-Options.html https://gcc.gnu.org/onlinedocs/gcc-3.2/gcc/Optimize-Options.... and also "Note on x86 and m68080 floating-point math" here https://gcc.gnu.org/wiki/FloatingPointMath https://gcc.gnu.org/wiki/FloatingPointMath
- nsajko 10y ago
  Nonstandard compiler options (ffast-math) are not relevant to this bug.
- caf 10y ago
  The relevant option is actually -fexcess-precision=standard, which is implied by -std=c99.
- kccqzy 10y ago
  With -O2 gcc actually compiles the program to xorl %eax, %eax ret That is, it optimized away everything.
  - acqq 10y ago
    Yes, that's the correct optimization. That code obviously can't perform printf("error\n") to manifest the bug, and obviously both additions were made compile-time producing the results that are equal.
pcwalton 10y ago
> If the author of the IEEE 754 standard says that the way you do floating-point arithmetic is obviously a bug, you're probably doing it wrong. That's interesting, because it seems to contradict his attitude toward FLT_EVAL_METHOD. Kahan is a strong supporter of FLT_EVAL_METHOD==2—in fact, he's quite possibly the single person who's most responsible for its existence—and FLT_EVAL_METHOD==2 has been resoundly rejected by compilers on newer architectures due to being totally unintuitive. (For those who don't know, FLT_EVAL_METHOD==2 means that "x∗y+z" can have a different result from "double tmp = x∗y; tmp + z", even if x, y, and z are doubles. It's a very similar situation to the bug, just with expression temporaries vs. explicit storage locations instead of register-memory spills.) In particular, Kahan wrote a widely-shared, hyperbolic essay called "Why Java Floating Point Hurts Everyone Everywhere" (which is still sometimes shared around today as a way to attack Java, despite being obsolete) castigating Java for having FLT_EVAL_METHOD==0. Java added the little-used "strictfp" in response, but history proved Java's original behavior right—with the advent of SSE on x86-64, GCC switched to FLT_EVAL_METHOD==0 by default, matching Java. Besides, Kahan's x87 FPU design is arguably responsible for this bug to begin with. The 80-bit FP register stack is designed makes sense when you consider ancient C compilers that didn't keep variables in registers unless explicitly told to do so. In that case, the register stack allows for a simple implementation of FLT_EVAL_METHOD==2. But if assignment to a variable doesn't result in a spill (as in any non-toy compiler since 1990?) then this becomes awkward, because register spills no longer coincide with those evaluation rules. This is just a design mistake in the x87 FPU, and it's a large part of the reason why we have scalar FP math in SSE.
- acqq 10y ago
  I've did a fast "reread" of "Why Java Floating Point Hurts Everyone Everywhere" and I can't see FLT_EVAL_METHOD mentioned anywhere. What I see is that Kahan argues that Java should let him use 80 bits if he would want to use it, and that the float arithmetic without specific notation for otherwise should be done at least "as in K&R C" always as doubles. Basically his argument is "don't limit what I can program to only your strange subset of what's already in the standard and in the hardware, and also provide the defaults that when the programmer doesn't specify the explicit behavior the behavior is safer numerically." Which in the case of the sample program can't be "one sum is different than another sum" when both times it's the double that is the explicit result of the sum operations. It's not "should it spill or not" there the compiler has to respect the size which the programmer specified. P.S. answer to pcwalton's post below: no it's certainly not "pretty much" "hey, you lose precision when you evaluate in fewer bits." It's: let me do the full IEEE standard in Java. It's: if the programmer doesn't specify what the sizes of the calculations in the expressions are, do these calculations with the biggest sizes (but predictably!), and if he specifies, respect what he specified(1). In the sample program we discuss, doubles as the results are explicit. These doubles can be then compared even in 80 bits but they will be still the same. It's not a "spill" it's an explicit cast that should happen. 1) "For example, tan(x) should delivers a double result no matter whether x is float or double, but tan((float)x) should deliver a float result provided a suitable tan method is available. Thus, a programmer who gives no thought to the question gets the safer default." This is obviously not "hey, you lose precision when you evaluate in fewer bits" he wants to be able to choose.
dbcurtis 10y ago
I have also met Kahan. He gives new meaning to the word "pedantic". The reality is that numerical anslysts want the computer do do exactly what they tell them to, in precisely the order described. The rest of the human race wants approximately the right answer reasonably quickly. (Floating point is a bad choice if you want exactly the right answer). And the people who make purchase decisions read benchmark reports they don't understand.
- smallnamespace 10y ago
  > The rest of the human race wants approximately the right answer reasonably quickly But numerical analysts really require being able to specify the semantics of calculations precisely in order to write libraries that mere mortals can use to get the approximate right answer. A quick perusal of any numerical analysis textbook should convince you that order of evaluation, order of truncation, overflow modes, and the like are very relevant for being able to write a library that has the right order of error. Something like Kahan summation will have vastly different results depend on whether intermediate calculations are done at a higher precision or whether the compiler is allowed to reorder expressions. In fact, an overly aggressive optimizing compiler will make Kahan summation simply fail [1]. If all you're doing is adding a few numbers together, by all means, you don't need or want to think about precise semantics. But if you've ever used BLAS or Numpy or any other math package, you very much depend on languages and compilers being pedantic as possible in spelling out how floating point gets evaluated. [1] https://en.wikipedia.org/wiki/Kahan_summation_algorithm https://en.wikipedia.org/wiki/Kahan_summation_algorithm
  - alfalfasprout 10y ago
    Actually, BLAS implementations usually do guarantee error bounds, but the implementation is almost always tailored for speed over precision. That said, you can always opt for extended precision at the cost of performance.
  - dbcurtis 10y ago
    >But numerical analysts really require being able to specify the semantics of calculations precisely in order to write libraries that mere mortals can use to get the approximate right answer. >A quick perusal of any numerical analysis textbook should convince you that order of evaluation, order of truncation, overflow modes, and the like are very relevant for being able to write a library that has the right order of error. I see your point. I've done more than peruse numerical analysis textbooks, and as an old-school CPU logic designer, I've forgotten more about implementing floating point than most programmers know to start with. But your point is not in conflict with mine. Numerical analysis requires precise operation ordering. Most people not only don't care, they don't have sufficient understanding of the issues that you could even convince them to care.
BuuQu9hu 10y ago
Kahan has another great paper, https://people.eecs.berkeley.edu/~wkahan/Mindless.pdf https://people.eecs.berkeley.edu/~wkahan/Mindless.pdf , which discusses the futility of attempting to rely on floating-point numbers.
EpicEng 10y ago
Interesting. This reminds me of a bug I had to track down when rewriting a bunch of Java image analysis code in C# (yeah I would have preferred neither language for this sort of thing, but that's a different story.) Java has a flag which allows you to specify how 32-bit floating point arithmetic happens (STRICT_FP or something.) C# does not, and was doing extended precision math on floats (80-bit register) and then truncating the result. This led to small differences in feature extraction output.

acqq 10y ago

And I believe he's right, if he was asked specifically what that program should do. If all the statements in the program are really just: #include <stdio.h> void test(double x, double y) { const double y2 = x + 1.0; if (y != y2) printf("error\n"); } void main() { const double x = .012; const double y = x + 1.0; test(x, y); } Then the comparison has to be done on doubles, and "but it's extended precision internally" mentioned in some answers is the wrong answer. All the variables in the program are stated to be doubles. If I'd write the x86 and x87 asm instructions I could easily write them so that they work as expected, even without SSE2. However if you compile with some "special" optimizations that are documented to relax some specific code generation rules, then it can be OK that these specific optimizations behave differently, but then: you asked for it by choosing that documented optimization. In the answers there, Richard Henderson wrote 2001-01-16 06:30:04 UTC: "You are seeing the results of excess precision in the FPU. Either change the rounding precision in the FPCR, or work around the problem with -ffloat-store." He's also right, if the -ffloat-store is documented to be turned off by -O and if it's documented that that flag also affects the precomputed calculations that are done in the compile time. If the second detail is not documented, the results could be unexpected, but again depends of what you expect from the optimizer. The properly written compiler and optimizer should perform the calculations during the code generation correct to the definition and not cut the corners. If the generated code was generated to actually perform the != during the run time, then the compile-time calculations are irrelevant, but the sane default should be something like -ffloat-store when not using SSE2 even with -O, if -O has the semantic "the optimizations that should not produce really unexpected effects." It sounds complicated, but it has its reasons. P.S. answer to nsajko's post below: I mention -ffloat-store not ffast-math. Read the description of ffloat-store please: https://gcc.gnu.org/onlinedocs/gcc-3.2/gcc/Optimize-Options.html https://gcc.gnu.org/onlinedocs/gcc-3.2/gcc/Optimize-Options.... and also "Note on x86 and m68080 floating-point math" here https://gcc.gnu.org/wiki/FloatingPointMath https://gcc.gnu.org/wiki/FloatingPointMath

nsajko 10y ago
Nonstandard compiler options (ffast-math) are not relevant to this bug.
caf 10y ago
The relevant option is actually -fexcess-precision=standard, which is implied by -std=c99.
kccqzy 10y ago
With -O2 gcc actually compiles the program to xorl %eax, %eax ret That is, it optimized away everything.
- acqq 10y ago
  Yes, that's the correct optimization. That code obviously can't perform printf("error\n") to manifest the bug, and obviously both additions were made compile-time producing the results that are equal.

pcwalton 10y ago

> If the author of the IEEE 754 standard says that the way you do floating-point arithmetic is obviously a bug, you're probably doing it wrong. That's interesting, because it seems to contradict his attitude toward FLT_EVAL_METHOD. Kahan is a strong supporter of FLT_EVAL_METHOD==2—in fact, he's quite possibly the single person who's most responsible for its existence—and FLT_EVAL_METHOD==2 has been resoundly rejected by compilers on newer architectures due to being totally unintuitive. (For those who don't know, FLT_EVAL_METHOD==2 means that "x∗y+z" can have a different result from "double tmp = x∗y; tmp + z", even if x, y, and z are doubles. It's a very similar situation to the bug, just with expression temporaries vs. explicit storage locations instead of register-memory spills.) In particular, Kahan wrote a widely-shared, hyperbolic essay called "Why Java Floating Point Hurts Everyone Everywhere" (which is still sometimes shared around today as a way to attack Java, despite being obsolete) castigating Java for having FLT_EVAL_METHOD==0. Java added the little-used "strictfp" in response, but history proved Java's original behavior right—with the advent of SSE on x86-64, GCC switched to FLT_EVAL_METHOD==0 by default, matching Java. Besides, Kahan's x87 FPU design is arguably responsible for this bug to begin with. The 80-bit FP register stack is designed makes sense when you consider ancient C compilers that didn't keep variables in registers unless explicitly told to do so. In that case, the register stack allows for a simple implementation of FLT_EVAL_METHOD==2. But if assignment to a variable doesn't result in a spill (as in any non-toy compiler since 1990?) then this becomes awkward, because register spills no longer coincide with those evaluation rules. This is just a design mistake in the x87 FPU, and it's a large part of the reason why we have scalar FP math in SSE.

acqq 10y ago
I've did a fast "reread" of "Why Java Floating Point Hurts Everyone Everywhere" and I can't see FLT_EVAL_METHOD mentioned anywhere. What I see is that Kahan argues that Java should let him use 80 bits if he would want to use it, and that the float arithmetic without specific notation for otherwise should be done at least "as in K&R C" always as doubles. Basically his argument is "don't limit what I can program to only your strange subset of what's already in the standard and in the hardware, and also provide the defaults that when the programmer doesn't specify the explicit behavior the behavior is safer numerically." Which in the case of the sample program can't be "one sum is different than another sum" when both times it's the double that is the explicit result of the sum operations. It's not "should it spill or not" there the compiler has to respect the size which the programmer specified. P.S. answer to pcwalton's post below: no it's certainly not "pretty much" "hey, you lose precision when you evaluate in fewer bits." It's: let me do the full IEEE standard in Java. It's: if the programmer doesn't specify what the sizes of the calculations in the expressions are, do these calculations with the biggest sizes (but predictably!), and if he specifies, respect what he specified(1). In the sample program we discuss, doubles as the results are explicit. These doubles can be then compared even in 80 bits but they will be still the same. It's not a "spill" it's an explicit cast that should happen. 1) "For example, tan(x) should delivers a double result no matter whether x is float or double, but tan((float)x) should deliver a float result provided a suitable tan method is available. Thus, a programmer who gives no thought to the question gets the safer default." This is obviously not "hey, you lose precision when you evaluate in fewer bits" he wants to be able to choose.

dbcurtis 10y ago

I have also met Kahan. He gives new meaning to the word "pedantic". The reality is that numerical anslysts want the computer do do exactly what they tell them to, in precisely the order described. The rest of the human race wants approximately the right answer reasonably quickly. (Floating point is a bad choice if you want exactly the right answer). And the people who make purchase decisions read benchmark reports they don't understand.

smallnamespace 10y ago
> The rest of the human race wants approximately the right answer reasonably quickly But numerical analysts really require being able to specify the semantics of calculations precisely in order to write libraries that mere mortals can use to get the approximate right answer. A quick perusal of any numerical analysis textbook should convince you that order of evaluation, order of truncation, overflow modes, and the like are very relevant for being able to write a library that has the right order of error. Something like Kahan summation will have vastly different results depend on whether intermediate calculations are done at a higher precision or whether the compiler is allowed to reorder expressions. In fact, an overly aggressive optimizing compiler will make Kahan summation simply fail [1]. If all you're doing is adding a few numbers together, by all means, you don't need or want to think about precise semantics. But if you've ever used BLAS or Numpy or any other math package, you very much depend on languages and compilers being pedantic as possible in spelling out how floating point gets evaluated. [1] https://en.wikipedia.org/wiki/Kahan_summation_algorithm https://en.wikipedia.org/wiki/Kahan_summation_algorithm
- alfalfasprout 10y ago
  Actually, BLAS implementations usually do guarantee error bounds, but the implementation is almost always tailored for speed over precision. That said, you can always opt for extended precision at the cost of performance.
- dbcurtis 10y ago
  >But numerical analysts really require being able to specify the semantics of calculations precisely in order to write libraries that mere mortals can use to get the approximate right answer. >A quick perusal of any numerical analysis textbook should convince you that order of evaluation, order of truncation, overflow modes, and the like are very relevant for being able to write a library that has the right order of error. I see your point. I've done more than peruse numerical analysis textbooks, and as an old-school CPU logic designer, I've forgotten more about implementing floating point than most programmers know to start with. But your point is not in conflict with mine. Numerical analysis requires precise operation ordering. Most people not only don't care, they don't have sufficient understanding of the issues that you could even convince them to care.

BuuQu9hu 10y ago

Kahan has another great paper, https://people.eecs.berkeley.edu/~wkahan/Mindless.pdf https://people.eecs.berkeley.edu/~wkahan/Mindless.pdf , which discusses the futility of attempting to rely on floating-point numbers.

EpicEng 10y ago

Interesting. This reminds me of a bug I had to track down when rewriting a bunch of Java image analysis code in C# (yeah I would have preferred neither language for this sort of thing, but that's a different story.) Java has a flag which allows you to specify how 32-bit floating point arithmetic happens (STRICT_FP or something.) C# does not, and was doing extended precision math on floats (80-bit register) and then truncating the result. This led to small differences in feature extraction output.