How to account for the antisymmetric nature of the gravitational force in a simulation of N gravitational bodies

I am attempting a simulation of the gravitational N-body problem in Python. My code for the acceleration on the ith body is of the same structure as

def acc(r,m,eps):
a = np.zeros((len(r),3))
for i in range(len(r)):
    for j in range(len(r)):
        ra2 = ((r[i,:]-r[j,:])**2).sum()
        if (i != j):
            a[i,:] += -(r[i,:]-r[j,:])*m[j]/(ra2**1.5) 
return a # return acceleration

which is found here http://wiki.tomabel.org/index.php?title=Gravitational_N-body_Problem

However in this format, are we not going to be performing unnecessary calculations as the force on particle i due to particle j, is just going to be the negative of the force on particle j due to particle i? How would we take this into account in order for the program to run faster? I was thinking of somehow taking an N x N array, filling half of it, and then taking the transpose, but am unsure as to how to do this, or if there is a better way.

Thanks very much

Solution

I would suggest something like this:

def acc(r,m,eps):
    a = np.zeros((len(r),3))
    for i in range(len(r))[:-1]:
        for j in range(len(r))[i+1:]:
            ra2 = ((r[i,:]-r[j,:])**2).sum()
            f= -(r[i,:]-r[j,:])*m[j]/(ra2**1.5) 
            a[i,:] += f
            a[j,:] += -f
    return a # return acceleration

In this way you update both i and j acceleration at the same time so you can always assume j > i which enables you to avoid double calculation.