Bad optimization of std::fabs()?

Question

Recently i was working with an application that had code similar to:

for (auto x = 0; x < width - 1 - left; ++x)
{
    // store / reset points
    temp = hPoint = 0;
    for(int channel = 0; channel < audioData.size(); channel++)
    {
        if (peakmode) /* fir rms of window size */
        {
            for (int z = 0; z < sizeFactor; z++)
            {
                temp += audioData[channel][x * sizeFactor + z + offset];
            }
            hPoint += temp / sizeFactor;
        }
        else /* highest sample in window */
        {
            for (int z = 0; z < sizeFactor; z++)
            {
                temp = audioData[channel][x * sizeFactor + z + offset];
                if (std::fabs(temp) > std::fabs(hPoint))
                hPoint = temp;
            }
        }
        .. some other code
    }
    ... some more code
}

This is inside a graphical render loop, called some 50-100 times / sec with buffers up to 192kHz in multiple channels. So it's a lot of data running through the innermost loops, and profiling showed this was a hotspot.

It occurred to me that one could cast the float to an integer and erase the sign bit, and cast it back using only temporaries. It looked something like this:

if ((const float &&)(*((int *)&temp) & ~0x80000000) > (const float &&)(*((int *)&hPoint) & ~0x80000000))
    hPoint = temp;

This gave a 12x reduction in render time, while still producing the same, valid output. Note that everything in the audiodata is sanitized beforehand to not include nans/infs/denormals, and only have a range of [-1, 1].

Are there any corner cases where this optimization will give wrong results - or, why is the standard library function not implemented like this? I presume it has to do with handling of non-normal values?

e: the layout of the floating point model is conforming to ieee, and sizeof(float) == sizeof(int) == 4

Answer 1

Well, you set the floating-point mode to IEEE conforming. Typically, with switches like --fast-math the compiler can ignore IEEE corner cases like NaN, INF and denormals. If the compiler also uses intrinsics, it can probably emit the same code.

BTW, if you're going to assume IEEE format, there's no need for the cast back to float prior to the comparison. The IEEE format is nifty: for all positive finite values, a<b if and only if reinterpret_cast<int_type>(a) < reinterpret_cast<int_type>(b)

Answer 2

It occurred to me that one could cast the float to an integer and erase the sign bit, and cast it back using only temporaries.

No, you can't, because this violates the strict aliasing rule.

Are there any corner cases where this optimization will give wrong results

Technically, this code results in undefined behavior, so it always gives wrong "results". Not in the sense that the result of the absolute value will always be unexpected or incorrect, but in the sense that you can't possibly reason about what a program does if it has undefined behavior.

or, why is the standard library function not implemented like this?

Your suspicion is justified, handling denormals and other exceptional values is tricky, the stdlib function also needs to take those into account, and the other reason is still the undefined behavior.

One (non-)solution if you care about performance:

Instead of casting and pointers, you can use a union. Unfortunately, that only works in C, not in C++, though. That won't result in UB, but it's still not portable (although it will likely work with most, if not all, platforms with IEEE-754).

union {
    float f;
    unsigned u;
} pun = { .f = -3.14 };

pun.u &= ~0x80000000;

printf("abs(-pi) = %f\n", pun.f);

But, granted, this may or may not be faster than calling fabs(). Only one thing is sure: it won't be always correct.

Answer 3

You would expect fabs() to be implemented in hardware. There was an 8087 instruction for it in 1980 after all. You're not going to beat the hardware.

Answer 4

How the standard library function implements it is .... implementation dependent. So you may find different implementation of the standard library with different performance.

I imagine that you could have problems in platforms where int is not 32 bits. You 'd better use int32_t (cstdint>)

For my knowledge, was std::abs previously inlined ? Or the optimisation you observed is mainly due to suppression of the function call ?

Answer 5

Some observations on how refactoring may improve performance:

• as mentioned, x * sizeFactor + offset can be factored out of the inner loops

• peakmode is actually a switch changing the function's behaviour - make two functions rather than test the switch mid-loop. This has 2 benefits:

  1. easier to maintain
  2. fewer local variables and code paths to get in the way of optimisation.
     
     

• The division of temp by sizeFactor can be deferred until outside the channel loop in the peakmode version.

• abs(hPoint) can be pre-computed whenever hPoint is updated

• if audioData is a vector of vectors you may get some performance benefit by taking a reference to audioData[channel] at the start of the body of the channel loop, reducing the array indexing within the z loop down to one dimension.

• finally, apply whatever specific optimisations for the calculation of fabs you deem fit. Anything you do here will hurt portability so it's a last resort.

Answer 6

In VS2008, using the following to track the absolute value of hpoint and hIsNeg to remember whether it is positive or negative is about twice as fast as using fabs():

int hIsNeg=0 ;
...
//Inside loop, replacing
//    if (std::fabs(temp) > std::fabs(hPoint))
//        hPoint = temp;
if( temp < 0 )
{
    if( -temp > hpoint )
    {
        hpoint = -temp ;
        hIsNeg = 1 ;
    }
}
else
{
    if( temp > hpoint )
    {
        hpoint = temp ;
        hIsNeg = 0 ;
    }
}
...
//After loop
if( hIsNeg )
    hpoint = -hpoint ;