How does C++ compute the absolute value of a complex number, preventing overflow?

Question

The C++ header <complex> provides abs(z) and norm(z).

The norm of a complex number z=x+iy is norm(z):=x^2+y^2.

The absolute value of z is abs(z):=sqrt(norm(z)).

However, the following example shows that abs(z) must be implemented differently, since it does not overflow although norm(z) does. At least, it does not overflow under g++ 6.2.1.

Is this non-overflow guaranteed by the standard? How is it achieved?

#include <iostream>
#include <complex>
typedef std::complex<double> complex_t;

int main()
{
    complex_t z = { 3e200, 4e200 };
    double a = abs(z);
    double n = norm(z);

    std::cout << a << " -> " << std::isinf(a) << "\n";
    std::cout << n << " -> " << std::isinf(n) << "\n";

    return 0;
}

Output:

5e+200 -> 0
inf -> 1

Answer 1

The std::complex::abs is equivalent to std::hypot function, which is indeed guaranteed to avoid overflow and underflow at intermediate stages of the computation.

Wikipedia page on Hypot function gives some insight on the implementation.

I'll quote the pseudocode just in case:

  // hypot for (x, y) != (0, 0)
double hypot(double x,double y)
{
    double t;
    x = abs(x);
    y = abs(y);
    t = min(x,y);
    x = max(x,y);
    t = t/x;
    return x*sqrt(1+t*t);
}