Tried this, and it doesn't help, the system remains sensitive to the order.

It turns out that using compensated summation (Increment or Decrement) is not correct: these methods are only valid if the addend is smaller that the current value. For accelerations we don't have any guarantees to that effect, and it's clear that reordering the bodies will affect the final result. Experiments with Mathematica show that reordering the bodies results in an error which is of 52 bits at the beginning (fine) but which occasionally jumps to 40 bits. When these errors accumulate the trajectories of the planets end up being affected significantly (remember that the stock system is very unstable).

Using full double precision fixes the problem, but we probably don't want to pay the cost of 10 additions/subtractions instead of 1 in the inner loop of the force computation. For a stable solar system these inaccuracies have no effect whatsoever.

Code pushed to pleroy/1351.

Would either pairwise summation, or Kahan summation help this?

See the reference to "compensated summation" in the comment above, which is another name for "Kahan summation". Kahan summation can yield large relative errors for ill-conditioned sums. That's one of these cases.

And pairwise isn't helping enough?

Pairwise is less precise than Kahan and has the same conditioning properties.

The numerical analysis above was mostly done to check that the sensitivity to order doesn't hide a bug. It doesn't. Done. The problem of accurately predicting the Kerbin system over long periods of time is spectacularly uninteresting as that system is extremely unstable (on a scale of months in the stock game).