When to use `std::hypot(x,y)` over `std::sqrt(x*x + y*y)`
Asked Answered
F

1

41

The documentation of std::hypot says that:

Computes the square root of the sum of the squares of x and y, without undue overflow or underflow at intermediate stages of the computation.

I struggle to conceive a test case where std::hypot should be used over the trivial sqrt(x*x + y*y).

The following test shows that std::hypot is roughly 20x slower than the naive calculation.

#include <iostream>
#include <chrono>
#include <random>
#include <algorithm>

int main(int, char**) {
    std::mt19937_64 mt;
    const auto samples = 10000000;
    std::vector<double> values(2 * samples);
    std::uniform_real_distribution<double> urd(-100.0, 100.0);
    std::generate_n(values.begin(), 2 * samples, [&]() {return urd(mt); });
    std::cout.precision(15);

    {
        double sum = 0;
        auto s = std::chrono::steady_clock::now();
        for (auto i = 0; i < 2 * samples; i += 2) {
            sum += std::hypot(values[i], values[i + 1]);
        }
        auto e = std::chrono::steady_clock::now();
        std::cout << std::fixed <<std::chrono::duration_cast<std::chrono::microseconds>(e - s).count() << "us --- s:" << sum << std::endl;
    }
    {
        double sum = 0;
        auto s = std::chrono::steady_clock::now();
        for (auto i = 0; i < 2 * samples; i += 2) {
            sum += std::sqrt(values[i]* values[i] + values[i + 1]* values[i + 1]);
        }
        auto e = std::chrono::steady_clock::now();
        std::cout << std::fixed << std::chrono::duration_cast<std::chrono::microseconds>(e - s).count() << "us --- s:" << sum << std::endl;
    }
}

So I'm asking for guidance, when must I use std::hypot(x,y) to obtain correct results over the much faster std::sqrt(x*x + y*y).

Clarification: I'm looking for answers that apply when x and y are floating point numbers. I.e. compare:

double h = std::hypot(static_cast<double>(x),static_cast<double>(y));

to:

double xx = static_cast<double>(x);
double yy = static_cast<double>(y);
double h = std::sqrt(xx*xx + yy*yy);
Foretoken answered 7/9, 2015 at 9:55 Comment(5)
I think you should also compare it with std::abs(std::complex<double>(x,y)) as in std::hypot pageHubie
Late, but the cppreference documentation also says as a note (so no guarantee by the standard) that "Implementations usually guarantee precision of less than 1 ulp (units in the last place)." x*x + y*y can lose a couple of bits of precision, if with round to nearest set. This means that std::sqrt(x*x+y*y) can be off by a bit or two. A better algorithm than std::sqrt(x*x+y*y) is needed to get that guarantee. (continued)Fionnula
To make matters worse, suppose you've mucked with the rounding? That will definitely get in the way of achieving that sub-ulp precision. hypot has to set the rounding so as to attain that accuracy and then restore the rounding back to your settings. This setting and resetting of rounding behavior is what makes std:hypot(x,y) considerably slower than std::sqrt(x*x+y*y).Fionnula
I enjoyed this question, but I still wanted to know why the performance disparity. #3765478 has a discussion of this. Specifically, https://mcmap.net/q/161274/-why-hypot-function-is-so-slow explains It for me.Epact
The sqrt function has the property that any relative error present in the input is halved in the result of the square root --- i.e. sqrt(x*(1+e)) ~=~ sqrt(x)*(1+e/2) --- (whereas squaring doubles it), so the square root method is not as bad it as would seem from the above. The extra runtime for hypot is partly made of up choosing among different methods to get extra precision, and steps to avoid over/underflow, but also special tests for inf (e.g. hypot(inf,NaN) -> inf, whereas the other approach gives you NaN).Dymphia
M
43

The answer is in the documentation you quoted

Computes the square root of the sum of the squares of x and y, without undue overflow or underflow at intermediate stages of the computation.

If x*x + y*y overflows, then if you carry out the calculation manually, you'll get the wrong answer. If you use std::hypot, however, it guarantees that the intermediate calculations will not overflow.

You can see an example of this disparity here.

If you are working with numbers which you know will not overflow the relevant representation for your platform, you can happily use the naive version.

Myosotis answered 7/9, 2015 at 10:11 Comment(7)
std::hypot's parameters are floating-pont types, so if you use int here it'll also be promoted to floating-point.Hubie
I might not have been 100% clear in my question. I'm only interested in difference between the two when double precision arguments are used. Because integer overflow is handled by the type promotion, not the hypot function. But something in hypot makes it 20x slower than the naive implementation (even with type promotion). Why is this and when will it save my tushie?Foretoken
As can be seen in the example here there is no disparity when using double precision arguments. Specifically I'm looking for double x,y where a disparity will appear.Foretoken
@EmilyL. just make the value large enough.Myosotis
@Myosotis I still don't see a difference.Foretoken
@Myosotis Now we're getting somewhere :) So when the intermediary result exceeds 1.7e308 then there will be a difference. If you edit your answer to include this and skip the integer overflow discussion. I'll be happy to accept it.Foretoken
@EmilyL. for a given platform, sure. More generally, when it exceeds std::numeric_limits<double>::max().Myosotis

© 2022 - 2024 — McMap. All rights reserved.