std::abs(std::complex) too slow

Question

Why running std::abs over a big complex array is about 8 times slower than using sqrt and norm?

#include <ctime>
#include <cmath>
#include <vector>
#include <complex>
#include <iostream>
using namespace std;

int main()
{
    typedef complex<double> compd;

    vector<compd> arr(2e7);
    for (compd& c : arr)
    {
        c.real(rand());
        c.imag(rand());
    }

    double sm = 0;
    clock_t tm = clock();
    for (const compd& c : arr)
    {
        sm += abs(c);
    }
    cout << sm << ' ' << clock() - tm << endl; // 5.01554e+011 - 1640 ms

    sm = 0;
    tm = clock();
    for (const compd& c : arr)
    {
        sm += sqrt(norm(c));
    }
    cout << sm << ' ' << clock() - tm << endl; // 5.01554e+011 - 154

    sm = 0;
    tm = clock();
    for (const compd& c : arr)
    {
        sm += hypot(c.real(), c.imag());
    }
    cout << sm << ' ' << clock() - tm << endl; // 5.01554e+011 - 221
}

Answer 1

I believe the two are not to be taken as identical in the strict sense.

From cppreference on std::abs(std::complex):

Errors and special cases are handled as if the function is implemented as std::hypot(std::real(z), std::imag(z))

Also from cppreference on std::norm(std::complex):

The norm calculated by this function is also known as field norm or absolute square.

The Euclidean norm of a complex number is provided by std::abs, which is more costly to compute. In some situations, it may be replaced by std::norm, for example, if abs(z1) > abs(z2) then norm(z1) > norm(z2).

In short, there are cases where a different result is obtained from each function. Some of these may be found in std::hypot. There the notes also mention the following:

std::hypot(x, y) is equivalent to std::abs(std::complex<double>(x,y))

In general the accuracy of the result may be different (due to the usual floating point mess), and it seems the functions were designed in such a way to be as accurate as possible.

Answer 2

The main reason is that abs handles underflow and overflow during intermediate computations.

So, if norm under/overflows, your formula returns an incorrect/inaccurate result, while abs will return the correct one (so, for example, if your input numbers are in the range of 10²⁰⁰, then the result should be around 10²⁰⁰ as well. But your formula will give you inf, or a floating point exception, because the intermediate norm is around 10⁴⁰⁰, which is out of range. Note, I've supposed IEEE-754 64-bit floating point here).

Another reason is that abs may give a little bit more precise result.

If you don't need to handle these cases, because your input numbers are "well-behaved" (and don't need the possible more precise result), feel free to use your formula.