n-dimensional Euclidean space calculation templates
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.