# 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 <vector>
#include <algorithm>
#include <stdexcept>

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

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

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

}


Here is some driver code to demonstrate usage.

// points.cpp
#include <iostream>
#include <vector>
#include "euclid.h"

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

#define say(x) std::cout << "" #x " = " << (x) << std::endl
int main()
{
std::vector<double> 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:

1. can the squared_distance calculation be made more efficient?
2. is there a nice way to reduce code duplication in the + and - operators?
3. I've tested with int, float, double and std::complex. What else might be useful?
4. is there any point to accommodating unsigned types?
5. should I do anything about possible numeric overflow or underflow?
6. 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 <typename T>
void check_size(std::vector<T> const &a, std::vector<T> 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 <typename T>
std::vector<T> operator+(std::vector<T> a, const std::vector<T>& b)
{
check_size(a, b);
std::transform(a.begin(), a.end(),
b.begin(),
a.begin(),
std::plus<T>());
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<double> 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<< 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.

-
Sorry about that :/ – Morwenn Apr 6 '14 at 9:52
I had completely forgotten about std::valarray and just now read the relevant bits of the standard about it. Seems tailor-made for what I'm trying to do, so I'll do some performance tests with std::vector, std::array and std::valarray and post results soon. – Edward Apr 6 '14 at 12:01
I think std::valarray was forgotten by its designers first, then by the users. When you start using it you realize you want more, but this more is simply not there :-( – iavr Apr 6 '14 at 18:02

Regarding duplication, you could have a single function for applying a function elementwize:

template< typename T, typename U, typename BinaryFunc >
auto elementwise_apply(const vector<T> &x, const vector<U> &y, BinaryFunc func)
-> vector<decltype(func(x[0], y[0]))>
{
// apply func to all arguments and return a vector
}


and then declare your other binary operations in terms of it

template< typename T>
vector<T> operator+(const vector<T> &x, const vector<T> &y)
{
return elementwise_apply(x, y, plus<T>());
}


and so forth.

I also agree with seeing if you can use array rather than vector -- the more that you can make happen at compile-time rather than run time, the better. Similarly, fewer dynamic allocations and fewer non-inlined functions would generally lead to more efficient code.

I would also suggest that you consider making a point struct that has an array member (or a vector member), rather than using a naked array, so that all of your functions will only applicable to objects meant to actually refer to points.

-

I am working on this topic so I have much to say, please bear with me. I am starting with your question (2), which is the most interesting. At the same time, I am making a series of abstraction and generalization steps. The final result will appear weird at first, but is quite more powerful than your current setup, in different ways. I will then discuss your remaining questions quite shortly. All my suggestions can be seen in a live example.

## Generalizations

Before coming to your question, let's look at a couple of limitations in the current approach

operator+(const std::vector<T>& a, const std::vector<T>& b)


(A) Why should the vectors have the same value type? Better, why have the same kind of containers in the first place?

(B) Why limit the operation to two arguments? there are unary operators out there, and there are n-ary functions as well.

Issue (A) can be dealt with by a few type functions:

template <typename C, typename T>
struct subs_type;

template <template <typename...> class C, typename T, typename... A, typename S>
struct subs_type<C<T, A...>, S> { using type = C<S, A...>; };

template <typename... A>
using container_common_type = std::common_type<typename A::value_type...>;

template <typename A, typename... B>
using container_result =
subs_type<A, typename container_common_type<A, B...>::type>;


In words, we are finding the common type of all value types of the given container types, and we are generating a new container of the same kind (e.f. std::vector with the first container and the common value type. Additional template parameters like allocators etc. are reused.

Issue (B) can be dealt with by a quite generic function object apply:

template <typename F>
struct apply : F
{
using F::operator();

template <typename A, typename... B>
typename container_result<A, B...>::type
operator()(const A& a, const B&... b) const
{
using R = typename container_result<A, B...>::type;

if (_or()(a.size() != b.size()...))
throw std::domain_error("vector operation must have equal length vectors");

R result;
result.reserve(a.size());
auto r = std::back_inserter(result);
transform(*this, r, a.begin(), a.end(), b.begin()...);
return result;
}
};


In words, apply applies function object of type F to input arguments of types A, B.... There are a few interesting generalizations here:

• It requires at least one input argument A, but otherwise the fully variadic definition allows unary, binary, or arbitrary n-ary functions.

• It derives function object F, overloads its operator(), and uses itself as a function when entering transform. This allows recursive application so that the same object can apply to scalars as well as any kind of container. We'll see an example below under your question 3.

However, it keeps the same logic of storing the result, using reserve and back_inserter. I am discussing alternatives later on.

The underlying operations required by apply are the following:

struct _or
{
constexpr bool operator()() const { return false; }

template <typename A, typename... B>
constexpr bool operator()(A&& a, B&&... b) const
{ return std::forward<A>(a) || operator()(std::forward<B>(b)...); }
};

struct _do { template <typename... A> _do(A&&...) { } };

template <typename F, typename R, typename I, typename E, typename... J>
void transform(const F& f, R r, I i, E e, J... j)
{
for (; i != e; _do{++r, ++i, ++j...})
*r = f(*i, *j...);
}


Function object _or is pretty much self-explanatory: it's an n-ary generalization of operator||.

Function transform is maybe the weirdest of all. It is more or less an n-ary generalization of std::transform. r is the output iterator, i,e are the current iterator and end respectively of the first container, and j... are the iterators of the remaining iterators. The function object and the output iterator come first, to allow for variadic input iterators in the sequel. All iterators are passed by value.

In words, transform has an output and a number of input iterators. All iterators advance in parallel. At each iteration, function object f is applied to the dereferenced values of all input iterators *i, *j... and the result is stored to the value of the output iterator *r.

Auxiliary struct _do allows for variadic, sequenced evaluation of its contained sub-expressions. I use it quite frequently, see Do here for more.

## Code duplication

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

Ok, whatever tools the language provides for generic programming, there comes a time when macros become inevitable. Macros? Before rushing to criticize, look at the following and please keep your mind open:

MYLIB_OP_UNARY(plus,  +)
MYLIB_OP_UNARY(minus, -)

MYLIB_OP_BINARY(sub, -)
MYLIB_OP_BINARY(mul, *)
MYLIB_OP_BINARY(div, /)
MYLIB_OP_BINARY(mod, %)

MYLIB_OP_BINARY(eq,  ==)
MYLIB_OP_BINARY(neq, !=)
MYLIB_OP_BINARY(gt,  >)
MYLIB_OP_BINARY(lt,  <)
MYLIB_OP_BINARY(ge,  >=)
MYLIB_OP_BINARY(le,  <=)

MYLIB_OP_UNARY(_not, !)
MYLIB_OP_BINARY(_and, &&)
MYLIB_OP_BINARY(_or,  ||)

MYLIB_OP_UNARY(bit_not, ~)
MYLIB_OP_BINARY(bit_and, &)
MYLIB_OP_BINARY(bit_or,  |)
MYLIB_OP_BINARY(bit_xor, ^)


These macros are defined and then only used once to generate the code for the required operators. You are then free to undefine them to avoid pollution. Be careful to use a prefix in the name of the macros to avoid collisions. I have assumed here you are making a library called mylib.

Macros are needed because each expands to around 20 lines of code (which is the minimal possible, despite abstractions), and I wouldn't like to copy this code 20 times or more. And if you think that operator== is an overkill, someone else will not.

Now here is how MYLIB_OP_BINARY is defined (MYLIB_OP_UNARY is similar):

#define MYLIB_OP_BINARY(NAME, OP)                               \
\
namespace operators {                                           \
struct NAME                                                     \
{                                                               \
template <typename A, typename B>                            \
constexpr auto                                               \
operator()(A&& a, B&& b) const                               \
-> decltype(std::forward<A>(a) OP std::forward<B>(b))        \
{ return std::forward<A>(a) OP std::forward<B>(b); }      \
};                                                              \
}                                                               \
\
template <typename A, typename B>                               \
typename container_result<A, B>::type                           \
operator OP(const A& a, const B& b)                             \
{                                                               \
return apply<operators::NAME>()(a, b);                       \
}                                                               \


First, inside a separate namespace operators we are capturing the meaning of each C++ operator into a function object. This is similar to std::plus etc. but

• applies to any C++ operator
• is non-template, allowing easier use

Then, for each operator, we define the corresponding "vectorized" operator for containers and scalars by just calling apply.

Note that this latter definition, although taking arbitrary types of input arguments, does not interfere with other potential overloads, because the result type

typename container_result<A, B>::type


is only valid if A,B are instances of class templates defining their own value_type (so, usually containers). Otherwise, other overloads are considered only. So we don't need enable_if in this; container_result is enough for substitution failure.

I hope you can now see that this macro is not possible to simplify further by existing language abstraction mechanisms. It would only be convenient to eliminate it if we could make it an one-liner, but then we'd loose all the nice generalizations exposed above.

In my work, I separately capture all C++ operators here and then define the "vectorized" counterparts here.

## Namespaces

Regardless of any other precautions, we are defining a number of very generic operators that extend the language in a way that may not be desirable in combination with (or compatible to) other code or other libraries. As an extra measure of precaution, define everything in a separate namespace. Check the live example for details. Again I have assumed you are building a library called mylib.

## Efficiency (accumulation)

(1) can the squared_distance calculation be made more efficient?

I don't want to elaborate too much here. A simple change would be fine:

return std::inner_product(a.begin(), a.end(), b.begin(), T(0),
std::plus<T>(), [](T x,T y){T z = y-x; return z*z;});


We just introduce a temporary to avoid doing the subtraction twice. But I would suggest to measure because most probably the optimizer is doing this anyway for us.

## Efficiency (element-wise operations)

I have left the issue of storing the result inside apply. Using std::vector::reserve() along with std::back_inserter means that

• We are limited to containers having reserve(), hence most probably to std::vector.
• We resize the container at each iteration. Even if we are not reallocating, we need to at least update the vector size.

There are numerous options here, but in most cases one would have go to a lower level of implementation, or use another kind of container.

For instance, if we could build a vector to allow uninitialized storage, we could give the appropriate size without initializing elements. transform would then only initialize each element instead of copying a new value. This can be implemented with a specialized allocator to be used with std::vector.

But then, why should apply decide what is the returned container type and how it should be initialized, resized, or copied? To me, the ultimate solution is to delay evaluation of the result by using expression templates. Then, expression a - b only returns a temporary object holding references to a, b and knowing the function type (like std::minus), while c = a - b does the actual work, initializing according to c's type.

As I said, I am working on this topic. Here is the sequence view (object behaving like a container) corresponding to apply as given above, and here is its iterator. I am afraid this code will be probably unreadable because it's using much more entities defined elsewhere, and because it has taken many more steps of abstraction and generalization than shown here. Libraries like boost::multi_array, boost::ublas or Eigen are similar in this respect (actually, in using expression templates and in being unreadable internally)!

I am sorry I am not giving any more concrete solution here, but I feel this is a huge topic, and if you start in this direction, you may soon realize you should be using an existing library in the first place.

Anyhow, the good news are that you now have a single generic implementation of all operations so whenever you decide to invest on a more efficient implementation, you'll only have to define it once.

## Other types

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

(4) is there any point to accommodating unsigned types?

You won't need anything for unsigned types. The definitions given here will work fine and conversions will be made according to language rules by std::common_type, even for std::complex.

Other types? Anything! For instance, why not other containers? Consider:

std::vector<std::vector<double> >
u{{0, 0, 0}, {3, 4, 5}},
v{{-1, -2, -3}, {5, 6, 7}};
say(u);
say(v);
say(-u);
say(+v);
say(u+v);
say(u-v);
say(v-u);


This works automatically now, for any level of nested containers, thanks to recursion. It wouldn't work with your code.

One "glitch" is that separation between vectors and scalars is done by checking for member type ::value_type as I said above. Unfortunately this is not quite appropriate because e.g. std::complex has this type but is not a container. I have a simple workaround in my live example by also requiring ::size_type, but all such issues will be resolved by C++ concepts.

One comes soon to realize that adding a scalar to a vector would also be convenient. This results in four combinations of operations. Then how about an n-ary function? This would result in 2^n combinations. I have a generic solution and here is a test sample of my code showing just a few of the things than can be done. But this goes way too far.

Another issue is that we have mostly ignored move semantics here. If we were to write this code in the most generic way, we should use rvalue references and forwarding everywhere. This is not so hard to do, but I have skipped it and used const& parameters almost everywhere.

## Type promotion

(5) should I do anything about possible numeric overflow or underflow?

For element-wise operations, nothing.

For accumulation operations like squared_distance, maybe yes. For instance, adding up 1000 values of type int might only fit into a long, and 1000 values of type float might only fit into a double. It would be convenient to have type promotions automatically done somewhere, but to decide exactly when a promotion is needed requires run time information, otherwise we may be promoting "just in case", which is contrary to the C++ principle

don't pay for what you don't use

I believe the user should be aware of such issues. Lazy evaluation may come to the rescue again: in an expression c = f(a, b), if f(a, b) is just a temporary intermediate object, then the actual work is only done when the type of c is known. This is exactly why std::accumulate and std::inner_product take another argument for the output. I also have a generic solution here, but again this goes too far.

For now, I would suggest to just use the right types for the input vectors (e.g. double) when an operation is expected to overflow the result. This is maybe too much a burden to the user, but we don't live in a perfect world.

## UPDATE: Corrections

I now realize that using common_type is not the right approach for finding the return type. It works for operators like +,* but not e.g. for boolean operators where the result is always bool. Moreover, common_type won't work with nested containers with different value types; it's simply not defined for such cases.

It's not too hard to make the return type depend on the operation F as well, using decltype as in the "capturing" function objects, or std::result_of. However, a further complication when considering nested containers with different value types is that substitution of the value type is not enough. For instance, std::vector<T> is actually

std::vector<T, std::allocator<T> >


so one has to substitute T for a new type, say S, in both the container and the allocator, to yield

std::vector<S, std::allocator<S> >


Assuming that containers follow this structure (in particular, they all have allocators), I have updated the previous code accordingly, but I am only giving here a link to a new live example to avoid making this answer even longer. For something more general without too much boilerplate, the containers should provide for such type substitution themselves; but to my knowledge, unfortunately this is not the case in STL.

We can now apply binary operators between containers of as diverse value types as

std::vector<std::vector<double> > u{{-2, 0, 0}, {3, 4, 5}};
std::vector<std::vector<int> >    v{{-1, 0, -3}, {5, 6, 5}};


as shown in the live example, along with some interesting cases like u < v.

Ok, that's it, I hope I will write no more!

-

I know that this won't answer any of your questions, but instead of using std::vector<T>, you should consider using std::array<T, N> to represent your data. Here is why this design could be more interesting:

• std::array<T, N> uses stack memory while std::vector<T> uses heap memory. I bet that a stack-allocated array will be a little bit more efficient than a heap-allocated vector.
• Mathematically speaking, it does not make any sense to change the size of a vector. Therefore, a fixed size container will also ensure that the size does not change.
• You won't have to bother with reserve anymore.
• In your functions, you make sure that the vector you add or subtract have the same size. If you use std::array<T, N> to store your data, you will get an implicit compile time check instead of an explicit runtime one:

template <typename T, size_t N>
std::array<T, N> operator+(const std::array<T, N>& a, const std::array<T, N>& b)
{
std::array<T, N> result;
for (size_t i = 0 ; i < N ; ++i)
{
result[i] = a[i] + b[i];
}
return result;
}


Moreover, since your loop has a static condition, it may be unrolled at compile time by the compiler.

-
Just a note: efficiency-wise (regardless of changing the container type), you are dropping the need to update the container size at each iteration, but at the same time you are implictly default-initializing all elements before entering the loop. This is exactly why Edward is using reserve and back_inserter, to avoid unnecessary initialization. Otherwise, one could say std::vector<T> result(N); and then write the same for loop as yours. Without measurements I cannot say which approach is better, but for me none is good enough (see discussion in my answer). – iavr Apr 6 '14 at 17:52
Plus, I bet you can't fit 10^7 elements on the stack, so it really depends on the application. The fact that STL algorithms are defined on iterators rather than containers has something to say (not on convenience, but at least on performance and genericity). – iavr Apr 6 '14 at 17:56
@iavr Those were some of my concerns, but since he mainly talked about 2D and 3D vectors, I considered that those weren't major drawbacks. I don't believe that many people perform computations in 10^7 dimensions Euclidean space. – Morwenn Apr 6 '14 at 18:15
@Morwenn You're right. Since modern physics says that all of spacetime can be described in 11 dimensions, more than that seems extravagant. – Edward Apr 6 '14 at 20:26
@Morwenn,Edward Well, algorithms for 2D/3D or even 11D spaces would rather use std:array if the number of dimensions is (compile-time) fixed. On the other hand, a color 1024x768 image has roughly 2,5M elements and I routinely do vector operations on such objects, or even worse. Element-wise operations, norms and squared distances remain the same in both cases. – iavr Apr 6 '14 at 20:50

# Results

Based on the feedback received, (thank you all for taking the time!) I wrote three different versions using std::vector, std::array and std::valarray. I then modified my test code to create a vector of one million randomly created 3D points (using double) and then called the squared_distance template function to populate a second results array with the squared distance between a fixed random point and each of the million other points. I then used the time command under Linux to record the User time taken for each of the three variants. Ten trials were done for each and the summary of the timing is shown below.

          valarray  array   vector
average     0.976   0.904   1.730
variance    0.000   0.002   0.020


As is clear from the results, both std::valarray and std::array are almost two times faster than my original std::vector version with std::array being fastest. That version of the code is shown below:

template <typename T, size_t N>
T squared_distance(const std::array<T,N>& a, const std::array<T,N>& b)
{
T accum(0);
for (int i=0; i<N; ++i)
{
T delta = b[i]-a[i];
accum += delta*delta;
}
return accum;
}

-