C++如何计算复数的绝对值，防止溢出？

Question

C++ 标头<complex> 提供abs(z) 和norm(z)。

复数z=x+iy的范数是norm(z):=x^2+y^2。

z 的绝对值是abs(z):=sqrt(norm(z))。

但是，下面的示例显示abs(z) 必须以不同的方式实现，因为尽管norm(z) 会溢出，但它不会溢出。至少在g++ 6.2.1下不会溢出。

标准是否保证了这种不溢出？它是如何实现的？

#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;
}

输出：

5e+200 -> 0
inf -> 1

Answer 1

std::complex::abs 等价于std::hypot 函数，确实可以保证在计算的中间阶段避免上溢和下溢。

Wikipedia page on Hypot function 提供了一些关于实现的见解。

我将引用伪代码以防万一：

  // 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);
}