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

user_13059497 0
Viewed 222

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.

compiler-optimization clang++ g++ x86 c++
2 Answers

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).

user_130594970

Comparable problem with intregrals (int) :

#include<cmath>

int fabsprod1(int a, int b) {
    return std::abs(a*b);
}
int fabsprod2(int a, int b) {
    return std::abs(a) * std::abs(b);
}
int fabsprod3(int a, int b) {
    return std::abs(std::abs(a) * std::abs(b));
}

Results in (using your options -O3 -std=c++2a -march=cannonlake):

fabsprod1(int, int):
        mov     eax, edi
        imul    eax, esi
        cdq
        xor     eax, edx
        sub     eax, edx
        ret
fabsprod2(int, int):
        mov     eax, edi
        cdq
        xor     eax, edx
        sub     eax, edx
        mov     edx, esi
        sar     edx, 31
        xor     esi, edx
        sub     esi, edx
        imul    eax, esi
        ret
fabsprod3(int, int):
        mov     eax, edi
        cdq
        xor     eax, edx
        sub     eax, edx
        mov     edx, esi
        sar     edx, 31
        xor     esi, edx
        sub     esi, edx
        imul    eax, esi
        ret

https://godbolt.org/z/tf3nZN

This contradicts your statement about "real/floating point" numbers.

Overall you shouldn't expect the compiler to shortcut math problems for you. That being said some optimizations is made possible. Please provide documentation or a comparable example where you see the optimaztion.

user_22615380
Related