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 &lt;vector&gt;
#include &lt;algorithm&gt;
#include &lt;stdexcept&gt;

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

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

template &lt;typename T&gt;
T squared_distance(const std::vector&lt;T&gt;&amp; a, const std::vector&lt;T&gt;&amp; 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&lt;T&gt;(), [](T x,T y){return (y-x)*(y-x);});

}


Here is some driver code to demonstrate usage.

// points.cpp
#include &lt;iostream&gt;
#include &lt;vector&gt;
#include "euclid.h"

template &lt;typename T&gt;
std::ostream&amp; operator&lt;&lt;(std::ostream&amp; out, const std::vector&lt;T&gt; &amp;v)
{
    out &lt;&lt; "{ ";
    for (auto p : v)
        out &lt;&lt; p &lt;&lt; ' ';
    return out &lt;&lt; "}";
}

#define say(x) std::cout &lt;&lt; "" #x " = " &lt;&lt; (x) &lt;&lt; std::endl
int main()
{
    std::vector&lt;double&gt; 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 -6 -8 }
squared_distance(origin,a) = 50
squared_distance(origin,b) = 14


I'm interested in general ideas for improvement or criticism of the existing code, but I also have some specific questions:


can the squared_distance calculation be made more efficient?
is there a nice way to reduce code duplication in the + and - operators?
I've tested with int, float, double and std::complex.  What else might be useful?
is there any point to accommodating unsigned types?
should I do anything about possible numeric overflow or underflow?
is there any point to having an #ifndef header guard in this file?
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 &lt;typename T&gt;
void check_size(std::vector&lt;T&gt; const &amp;a, std::vector&lt;T&gt; const &amp;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 &lt;typename T&gt;
std::vector&lt;T&gt; operator+(std::vector&lt;T&gt; a, const std::vector&lt;T&gt;&amp; b)
{
    check_size(a, b);
    std::transform(a.begin(), a.end(), 
                   b.begin(), 
                   a.begin(), 
                   std::plus&lt;T&gt;());
    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&lt;double&gt; 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&lt;&lt; overload (minutely modified to take a std::valarray parameter instead of a std::vector).

Considering your more specific questions:

specific questions:


  
  can the squared_distance calculation be made more efficient?
  


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.