0
\$\begingroup\$

I would like to call std::hypot with many arguments (not just 2 or 3).

Mandatory requirements:

  • fast
  • run at compile time when possible

Non-mandatory requirement:

  • avoid overflow

#include <cmath>
#include <iostream>


template <typename T>
constexpr T sum2(const T v) {
  return v * v;
}

template<typename T1, typename... Ts>
constexpr T1 sum2(const T1 v1, const Ts... rest)
{
  return v1 * v1 + sum2(rest...);
}

template <typename T1, typename... Ts>
constexpr T1 qsum(const T1 v1, const Ts... v)
{
  return sqrt(sum2(v1, v...));
}


int main()
{
  double a = 10.;
  std::cout << qsum(1., a) << std::endl;
}

One bad thing is when I ask qsum(1, 1.) since the first argument is int and so the output. I would like that in this case the inputs are cast to the safer common type.

\$\endgroup\$
  • 1
    \$\begingroup\$ Just as a note: POSIX specifies that underflow may only occur when both arguments are subnormal and the correct result is also subnormal **(this forbids naive implementations)** Also note the standard one only works on values of the same type (so standard promotion is forced before the function is called. So your template should only have one T. \$\endgroup\$ – Martin York May 24 '16 at 16:06
  • \$\begingroup\$ Lots more info here: http://en.cppreference.com/w/cpp/numeric/math/hypot \$\endgroup\$ – Martin York May 24 '16 at 16:06
2
\$\begingroup\$

First of all, let me just say that qsum is a really bad name for a generic hypot function. I would have called it... hypot! In fact, I'm going to do that in this answer, because it'll make the code a lot easier to read.

Well, I'll call it my_hypot, because we'll be including <cmath>, and I don't want any confusion over whether we're accidentally referring to the standard library's 2-argument hypot.


You have an anti-pattern in your code: you're writing things like T sum2(const T v) instead of T sum2(T v). The former is exactly equivalent to the latter as far as the caller of the function is concerned; but for the implementor of the function, the latter is strictly better:

  • It's less typing, and more importantly, less reading for the maintainer/reviewer.
  • It doesn't inhibit move semantics.
  • For these reasons, it's idiomatic C++; and therefore you should follow the idiom because breaking idioms for no reason is distracting to the maintainer/reviewer.

I think you could simplify your code by "thinking in C++17" and then working backwards to lower the resulting code back down to C++11.

In C++17, with fold-expressions, your function would be just

template<typename... Ts> constexpr auto my_hypot(Ts... vs) {
    return sqrt(( (vs * vs) + ... ));
}

In C++14, though, we don't have fold-expressions, so we're going to have to replace the construct (pattern + ...) with some kind of function call; say, sum(pattern...). So our C++14 code should look like

template<typename... Ts>
constexpr auto my_hypot(Ts... vs) {
    return sqrt(sum( (vs * vs)... ));
}

In C++11, we need to use the "trailing decltype" trick to get the compiler to accept that code, since C++11 can't deduce return types (except in lambdas).

template<typename... Ts>
constexpr auto my_hypot(Ts... vs) -> decltype(sqrt(sum( (vs * vs)... ))) {
    return sqrt(sum( (vs * vs)... ));
}

And there we go! Other than writing the sum function, this is a very straightforward translation of the code you want to write, into code that should be acceptable to a C++11 compiler.


However, I think your code as it stands is a bit sketchy, because of the overflow issue. my_hypot(65536, 65536, 1.0) should yield 92681.9, not undefined behavior! To mitigate that issue, I agree with you that casting all of the parameters to common_type would be a good idea. Let's implement that.

#include <cmath>
#include <iostream>
#include <type_traits>

template<typename CT> constexpr CT sum() { return {}; }

template<typename CT, typename T, typename... Ts>
constexpr CT sum(T v, Ts... vs)
{
    return v + sum<CT>(vs...);
}

template <typename... Ts>
constexpr auto my_hypot(Ts... vs) -> typename std::common_type<Ts...>::type
{
    using CT = typename std::common_type<Ts...>::type;
    return sqrt(sum<CT>( (CT(vs) * CT(vs))... ));
}

int main()
{
  std::cout << my_hypot(10., 10) << std::endl;
  std::cout << my_hypot(10, 10, 10.) << std::endl;
  std::cout << my_hypot(10, 10) << std::endl;
}

This is all the same straightforward code as before, with the possible exception of sum<CT>(). Here we're making an overload set of sum<CT>() (with no arguments) and sum<CT>(T, Ts...) (with at least one argument). By construction, both T and all the elements of Ts... will always be CT; we could static_assert that, if we wanted to. sum<CT>() returns a value-initialized CT (i.e., the "zero" of whatever type CT happens to be); I used {} instead of 0 just in case CT is some wacky bignum type that doesn't have an implicit conversion from int. Depending on your assumptions, 0 might be the safer choice.


The really nice thing about this code is that because we constructed it by starting with C++17 and moving backward, it's easy to take this code and migrate it forward to C++14 (by removing the trailing return types) and then to C++17 (by replacing the implementation of sum with a one-liner, or inlining it altogether) as your compiler catches up. Looking at your original code, I think it would be a little bit harder to disentangle the C++11isms from the "business logic" of how to compute hypot. For example, you've got the subexpression v * v repeated in two different places, for no good reason that I can see. You could have factored that out... or rather, not factored it in in the first place! :) Don't repeat yourself.

\$\endgroup\$
  • \$\begingroup\$ qsum is "sum in quadrature" or "euclidean distance". hypotenuse is just for triangles... \$\endgroup\$ – Ruggero Turra May 25 '16 at 8:53
  • \$\begingroup\$ Removing the trailing return type in my_hypot might change the output of this function. In fact, your implementation is wrong. Suppose you call my_hypot(1, 1). Then, CT is int and hence, in presence of your trailing return type, the output will be implicitly casted to int; which is not the desired behavior. So, you need to change the trailing return type to be the type of the result of the std::sqrt invocation. \$\endgroup\$ – 0xbadf00d Apr 1 '17 at 19:18
  • \$\begingroup\$ And let me note that your code still doesn't handle possible overflow issues as std::hypot does. Casting to std::common_type_t<Ts...> won't help in any practical application where usually all Ts are the same. \$\endgroup\$ – 0xbadf00d Apr 1 '17 at 19:29
1
\$\begingroup\$

I haven't tried to render it in constexpr form, but I would consider an extended version of the Moler and Morrison algorithm for computing a hypotenuse robustly. Section 2 of the paper shows the basic algorithm for 2 dimensional computation, and section 3 shows how to extend the 2D version to work with arbitrary dimensions.

For what it's worth, I posted a C implementation of the 2D algorithm to Stack Overflow some years ago.

I would take particular note of the fact that you can typically use a fixed number of iterations of the 2D algorithm, which I believe will probably make the code somewhat easier to render in a form that can be executed at compile time if needed (though I'll admit I'm uncertain about that).

The primary advantage of this algorithm is being much more immune to overflow (and underflow). If the actual distance won't fit in your destination type, then overflow becomes almost unavoidable, of course, but that's about the only time this will overflow (or underflow).

This is also relatively fast, at least as a general rule. In my experience, a single square root (alone) is often slower than computing the hypotenuse as a whole with this (though that testing has only been with the 2D version, not the multi-dimensional version). On the other hand, if you have a lot of dimensions, you may lose that advantage. A naive method does a lot of multiplications and additions, then one square root at the end. This does something in the same general speed range as a square root for each term instead.

By contrast, more obvious methods involve squaring the inputs, then summing the squares. This can result in overflow if the inputs are large or underflow if they're small. If we have a large number of dimensions, the summation should probably be done with something like Kahan's summation algorithm to maintain precision. That typically slows the summation by a factor of 2 (or so) compared to naive summation.

\$\endgroup\$
  • \$\begingroup\$ Upvoted for interestingness, but I'm not sure this counts as a code review. ;) \$\endgroup\$ – Quuxplusone May 26 '16 at 18:20
0
\$\begingroup\$

I tend to agree with Loki's comments if you wish to duplicate std::hypot. But, finding a "safer common type" is way more fun than having a common type for all arguments right from the start.

For this, you could use a trailing return type with decltype to make each expression-returning template function return its expression's promoted type, rather than whatever type the function's first argument has.

The example below is kind of ugly because each function's expression appears in the function body and in the decltype, so there is room for improvement, but it demonstrates the basic idea.

This prints 10.0499 on both lines, where the original would print 10.0499, then 10, because the second qsum call's first argument is an int.

#include <cmath>
#include <iostream>

template <typename T>
constexpr T sum2(const T v) {
  return v * v;
}

template<typename T1, typename... Ts>
constexpr auto sum2(const T1 v1, const Ts... rest) -> decltype(v1 * v1 + sum2(rest...))
{
  return v1 * v1 + sum2(rest...);
}

template <typename T1, typename... Ts>
constexpr auto qsum(const T1 v1, const Ts... v) -> decltype(sqrt(sum2(v1, v...)))
{
  return sqrt(sum2(v1, v...));
}

int main()
{
  std::cout << qsum(1.0, 10) << std::endl;
  std::cout << qsum(10, 1.0) << std::endl;
}
\$\endgroup\$
  • \$\begingroup\$ is the sqrt in decltype really executed? \$\endgroup\$ – Ruggero Turra May 25 '16 at 6:51
  • \$\begingroup\$ @RuggeroTurra no, the argument to decltype is never evaluated. However, despite having no value (because it's not evaluated), it certainly has a type (which is all decltype cares about). \$\endgroup\$ – Quuxplusone May 25 '16 at 8:01
  • \$\begingroup\$ In C++14, you could reasonably omit the trailing return types and just use auto return type deduction. The -> decltype(...) stuff matters for SFINAE purposes, but I'd argue that you could reasonably ignore SFINAE if you wanted to. \$\endgroup\$ – Quuxplusone May 25 '16 at 8:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.