Take the 2-minute tour ×

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;
    sum += pow(p[i]-q[i],2);

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.

share|improve this question
add comment

2 Answers

up vote 1 down vote accepted

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!

share|improve this answer
In addition to the vectorizing, you could also try performing the calculation on multiple cores, if that fits your use case. –  user1118321 Aug 24 '12 at 4:35
Unfortunately the Euclidean Norm functions provided in BLAS dont actually do what I want. The Euclidean Norm, or Euclidean Length, is equivalent to calculating the distance between a point in N dimensional space and the origin, not two different points. –  user1357607 Aug 24 '12 at 5:17
@user1357607 Yes, you have to subtract the two vectors with a *axpy first. The code would be something like *axpy(N, -1, p, 1, q, 1); answer = *nrm2(N, p, 1);. –  nevsan Aug 24 '12 at 5:20
Thanks. At the risk of being annoying, is there also a function for multiplying two vectors together, element by element? –  user1357607 Aug 24 '12 at 5:29
add comment

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.

share|improve this answer
add comment

Your Answer


