Quickest way to calculate the Euclidean-Distance between two points in N-Dimensional space

Question

I have to calculate the Euclidean Distance between two points in N-Dimensional space, and speed is crucial. I have two C style float arrays representing the two points in N-Dimensional space.

The formula for the distance between them is (^ just means to the power of, not XOR): sqrt(sum((p1-q1)^2 + (p2-q1)^2 + .... (pn-qn)^2))

My current code looks like this:

sum = 0;
for(int i=0;i<N;++i){
    sum += pow(p[i]-q[i],2);
sqrt(sum)

This code is quite slow and I was wondering if there any libraries to speed this up? I imagine someone has written a quick library on performing maths operations on arrays in c that lets me quickly do elementwise operations on arrays.

Edit: In answer to nevsan, I am doing many calculations with a small N, approximately 10 or 20.

Answer 1

Definitely get rid of pow(). A large part of optimization of this depends on how you are using it. Are you doing this once for very large N and it takes too long? Or, more likely, are you doing this many times in a tight loop?

If you are using very large N (>1000 or so), there are highly optimized numerical libraries that can do this. BLAS, for example, has an *nrm2 function that will compute the euclidean norm (dnrm2, snrm2, cnrm2, znrm2, depending on the data type [single, double, complex single, complex double]). GotoBLAS is probably the fastest out there for certain processor architectures. MKL features Intel's hand tuned BLAS implementation, but it's not free. Finally, ATLAS is a self-tuning library implementing BLAS.

If you have a tight loop with small or not quite large N, then you may have to do some hand-tuning to get it faster. You could turn on auto-vectorization with the -O3 or -ftree-vectorize compiler flags. You could also vectorize by hand, but it can be painful to learn how to do this.

You could do loop unrolling (that is, divide N up into chunks of, say, 4 and explicitly write out the computation for 4 successive values inside the for loop body. This has the effect of tricking the compiler to use more registers for immediate computation---and registers are the fastest form of memory that you have to work with. Also, you may be able to take advantage of prefetching (reading a stretch of data with one memory access call).

Another thing to do in this situation is to try overwriting one of your inputs. That is, maybe you could get away with writing the output into p or q. This helps because the positions of p that you compute will still be in the cache when you are ready to write. Caches often won't write the data to memory unless they absolutely have to---one reason is that the cache line is needed and we need to kick the last one out. You use fewer cache lines by writing to one of your inputs.

There are a half million other things to try, but I think I'll stop here. Good luck!

Answer 2

I would never use pow() - my guess without profiling is that this is slowing you down lots.

You need to make a temp, then square that.

double diff = p[i] - q[i];
sum += diff*diff;

the sqrt is a little slow, but the only options here are some approximations. If you have N > than about 10 likely sqrt will not be the bottleneck.

There are also libraries like boost, etc, that could likely speed this up, but first try getting rid of pow(). Remember that diff*diff is one floating point instruction, where pow() is a whole program designed for non integer powers, etc.