One fairly simple way to get much better performance is to only consider nearby neighbors. You can find your nearby neighbors by storing lists of particles by their integer coordinates in a hash table, and looking through your list to find your closest neighbors.
See my hashtable
gravity simulation, which immediately scales to over 1000
particles. But there's a problem: the galaxy can't hold
together anymore, even if we artificially inflate the
gravitational constant, because the long-range forces that bind
one side of the galaxy to the other are never evaluated.
The hashtable version gives the wrong answer because it ignores
long-range interactions. There are several neat ways to
approximate these interactions.
A very simple approach is to sum
all the particles into one "megaparticle", which will end up
at the center of the galaxy. We can then compute each
particle's interaction with its local neighbors, and then with the
megaparticle. This doesn't even conserve momentum, because
the influence of a particle on the megaparticle isn't symmetric,
but at least the galaxy holds together!
A more accurate version of this is called Barnes-Hut:
for neighboring particles, we compute full particle-particle
interactions. To interact with a group of distant particles,
we lump the group into a megaparticle. To accurately define
which particles form a "group", we impose a spatial tree across
the simulation, using the coarse levels of the tree to approximate
distant regions, and increasingly fine levels for more nearby
regions.