Difference between std::fabs(a * b) and std::fabs(a) * std::fabs(b)

Question

I am working on some numerical code and I was looking at the compiler output. One particular case struck me as odd:

In real numbers, it holds that abs(a) * abs(b) = abs(a * b). I would expect the same to hold in floating point numbers. However, the optimization is performed neither by clang nor by g++ and I wonder whether I am missing some subtle difference there. Both compilers do however realize that abs(abs(a) * abs(b)) = abs(a) * abs(b).

Here is the relevant piece of code:

#include<cmath>

double fabsprod1(double a, double b) {
    return std::fabs(a*b);
}
double fabsprod2(double a, double b) {
    return std::fabs(a) * std::fabs(b);
}
double fabsprod3(double a, double b) {
    return std::fabs(std::fabs(a) * std::fabs(b));
}

And here is the confusing compiler output in godbolt with gcc-10.1 (current stable version as of writing this) and -O3: https://godbolt.org/z/ZEFPgF

Notably, even with -Ofast, which as far as I understand is more lenient with the transformations that are allowed, this optimization is not performed.

As was pointed out by @Scheff in the comments, double and float are not real numbers. But I also fail to see where corner cases with float types, such as getting Infinity or NaN as argument, could produce different outputs.

Answer 1

I believe I have found a counter example. I post this as a separate answer, because I don't think that this is at all analogous to the case for integers.

In the cases I considered, I missed that it is possible to change the rounding mode for floating point arithmetic. Problematically, GCC seems to to ignore that when he (I guess) optimizes "known" quantities at compile time. Consider the following code:

#include <iostream>
#include <cmath>
#include <cfenv>

double fabsprod1(double a, double b) {
    return std::fabs(a*b);
}
double fabsprod2(double a, double b) {
    return std::fabs(a) * std::fabs(b);
}

int main() {
        std::fesetround(FE_DOWNWARD);
        double a  = 0.1;
        double b = -3;
        std::cout << std::hexfloat;
        std::cout << "fabsprod1(" << a << "," << b << "): " << fabsprod1(a,b) << "\n";
        std::cout << "fabsprod2(" << a << "," << b << "): " << fabsprod2(a,b) << "\n";
#ifdef CIN
        std::cin >> b;
#endif
}

Output differs, depending on whether I compile with

g++ -DCIN -O1 -march=native main2.cpp && ./a.out

or

g++ -O1 -march=native main2.cpp && ./a.out

Notably, it only takes O1 (what I would consider completely reliable) to change the output in a way that does not seem reasonable to me.

With -DCIN the output is

fabsprod1(0x1.999999999999ap-4,-0x1.8p+1): 0x1.3333333333334p-2
fabsprod2(0x1.999999999999ap-4,-0x1.8p+1): 0x1.3333333333333p-2

without -DCIN the output is

fabsprod1(0x1.999999999999ap-4,-0x1.8p+1): 0x1.3333333333334p-2
fabsprod2(0x1.999999999999ap-4,-0x1.8p+1): 0x1.3333333333334p-2

Edit: Peter Cordes (thank you for the comment) pointed out, that this surprising result was due to my failure in telling GCC to respect the change of rounding mode. By building with the following command, the expected results are achieved:

g++ -O1 -frounding-math -march=native main2.cpp && ./a.out

(works with O2 and O3 as well on my machine).