n-dimensional Euclidean space calculation templates

I have been working with C++11 code that uses std::vector[] to store coordinates. Most often this code uses 2D or 3D but it occurred to me that it may be generally useful for n-dimensional Euclidean distance calculations.

Because these calculations are used within a simulation, they should be as fast as practical without sacrificing generality.

// euclid.h
#include 
#include 
#include 

template 
std::vector operator+(const std::vector& a, const std::vector& b)
{
    if (a.size() != b.size())
        throw std::domain_error("vector addition must have equal length vectors");
    std::vector result;
    result.reserve(a.size());  
    std::transform(a.begin(), a.end(), b.begin(), 
                   std::back_inserter(result), std::plus());
    return result;
}

template 
std::vector operator-(const std::vector& a, const std::vector& b)
{
    if (a.size() != b.size())
        throw std::domain_error("vector subtraction must have equal length vectors");
    std::vector result;
    result.reserve(a.size());
    std::transform(a.begin(), a.end(), b.begin(), 
                   std::back_inserter(result), std::minus());
    return result;
}

template 
T squared_distance(const std::vector& a, const std::vector& b)
{
    if (a.size() != b.size())
        throw std::domain_error("squared_distance requires equal length vectors");
    return std::inner_product(a.begin(), a.end(), b.begin(), T(0), 
        std::plus(), [](T x,T y){return (y-x)*(y-x);});

}

Here is some driver code to demonstrate usage.

// points.cpp
#include 
#include 
#include "euclid.h"

template 
std::ostream& operator &v)
{
    out  origin{0, 0, 0}, a{3, 4, 5}, b{-1, -2, -3},g{0,0,0,0};
    say(origin);
    say(a);
    say(b);
    say(a+b);
    say(a-b);
    say(b-a);
    say(squared_distance(origin,a));
    say(squared_distance(origin,b));
}

The output from the program looks like this:

```
origin = { 0 0 0 }
a = { 3 4 5 }
b = { -1 -2 -3 }
a+b = { 2 2 2 }
a-b = { 4 6 8 }
b-a = { -4

To reduce redundancy, I'd probably move the check for the inputs being the same size into a function template by itself, and just invoke that from the other three:

template 
void check_size(std::vector const &a, std::vector const &b) {
    if (a.size() != b.size())
        throw std::domain_error("vector addition must have equal length vectors");
}

I'd then at least test with passing one of the input parameters by value, using it as the destination for the calculation, and returning it:

template 
std::vector operator+(std::vector a, const std::vector& b)
{
    check_size(a, b);
    std::transform(a.begin(), a.end(), 
                   b.begin(), 
                   a.begin(), 
                   std::plus());
    return a;
}

Although this doesn't always improve speed, it does often enough to be worth testing/profiling to see how well it works in this case.

Another possibility to consider would be using std::valarray instead. It's the forgotten step-child (so to speak) of the C++ standard, but it was designed specifically to support fast numeric calculations. It already defines equivalents of your + and - operators, as well as a sum and * operators, so your code comes out something like this:

std::valarray origin{ 0, 0, 0 }, a{ 3, 4, 5 }, b{ -1, -2, -3 }, g{ 0, 0, 0, 0 };
say(origin);
say(a);
say(b);
say(a + b);
say(a - b);
say(b - a);
say(((a - origin) * (a - origin)).sum());
say(((b - origin) * (b - origin)).sum());

Since std::valarray defines all the operators we're using, the only other code we need is the say macro and the operator


std::valarray and/or std::array might help--but for a relatively small number of dimensions, I wouldn't expect to see much difference.

• is there a nice way to reduce code duplication in the + and - operators?

The code above attempts to answer this, to at least some degree.

• I've tested with int, float, double and std::complex. What else might be useful?

• is there any point to accommodating unsigned types?

I don't see a lot of point in unsigned types for this task. Yes, distances are always positive, but you can pretty easily end up with a negative number from an intermediate calculation, and you don't want that wrapping around to a large number as it would with an unsigned type.

• should I do anything about possible numeric overflow or underflow?

Harder to say. It basically depends on the use to which you're putting the code. If you need the distance, there are ways of computing it that avoid the larger magnitude of the squared distance (even as an intermediate), and are fairly fast to compute as well--but they're still slower than computing the squared distance, in which case the final result is also often the single largest value you deal with, so there's not a lot you can about overflow.

In floating point, subtraction is one of the prime culprits that can/will lead to precision loss. At least in theory, you might prefer to avoid it if possible, but for finding a squared distance I don't know of a lot of alternatives either.

• is there any point to having an #ifndef header guard in this file?

The compiler will (or should, anyway) give errors if you include it twice in the same translation unit, so it wouldn't hurt. Of course, you're unlikely to include it twice in the same translation unit directly, but its getting included via two other headers wouldn't be all that rare an occurrence.

I was going to mention std::array` as well, but while I was writing this, I see @Morwenn has written one about that already.

StackExchange Code Review Q#46410, answer score: 8