For an exercise in university, we had to implement a [Leapfrog integrator](http://en.wikipedia.org/wiki/Leapfrog_integration) with exact Newtonian forces in Python. The course is over and our solutions were all more than good enough, but I wondered if/how one could improve the performance of the force calculation even more.

The bottleneck is to calculate all forces (aka accelerations):

![a<sub>i</sub> = Σ<sub>j≠i</sub> Gm<sub>j</sub> / |r<sub>1</sub>-r<sub>2</sub>|<sup>3</sup> * (r<sub>1</sub>-r<sub>2</sub>)](https://i.stack.imgur.com/HhdJE.png)

for a large (1000 and bigger) number of particles N (i,j<N). 

Here r1 and r2 are the 3-dimensional vectors of the position of the particles stored in an ndarray of shape (N,3) and Gm is the particle mass times the gravitational constant which I saved in an ndarray of shape (N).

The fastest version I found so far is the following:

```
def a(self):
    sep = self.r[np.newaxis, :] - self.r[:, np.newaxis]
    # Calculate the distances between all particles with cdist
    # this is much faster than by hand
    dists = cdist(self.r, self.r)
    scale =dists*dists*dists
    # set diagonal elements of dist to something != 0, to avoid division by 0
    np.fill_diagonal(scale,1)
    Fsum = (sep/scale.reshape(self.particlenr,self.particlenr,1))*self.Gm[:,None]
    return np.add.reduce(Fsum,axis=1)

```

But it bugs me that this is probably not the fastest version. The first line seems to be too slow when comparing to cdist which is doing essentially the same calculation. Also, this solution does not use the symmetry of switching r1 and r2 in the problem and calculates all elements twice.

Do you know any performance improvements (without changing the force-calculation to some approximation or changing the programming language)?