Just a small bit of (possibly naive) input. I don't do game programming, but what I am feeling is that your fundamental bottleneck is the gravity-due-to-gravity calculation. Instead of iterating over each object X and then finding the gravitational effect from each object Y and adding it, you can take each pair X,Y and find the force between them. That should cut the number of gravity calculations from O(n^2). Then you will be doing a lot of addition (O(n^2)), but this is normally less expensive.

Also at this point you can implement rules such as "if the gravitational force will be less than \epsilon because these bodies are too small, set the force to zero". It may be advantageous to have this structure for other purposes too (including collision detection).