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

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

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

c++

numeric

floating-accuracy

complex-numbers