Как 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); 
}