Looking for Runge-Kutta 8th order in C/C++

If you're doing celestial mechanics over long time scales, using a classical Runge-Kutta integrator will not preserve energy. In that case, using a symplectic integrator would probably be better. Boost.odeint also implements a 4th-order symplectic Runge-Kutta scheme that would work better for long time intervals. GSL does not implement any symplectic methods, as far as I can tell.

I would like to add that while what Geoff Oxberry suggests for long-term integration (using symplectic integrators) is true, in some cases it won't work. More specifically, if you have dissipative forces, your system does not preserve the energy anymore, and therefore you cannot resort to symplectic integrators in that case. The person asking the question was talking about low Earth orbits, and such orbits exhibit a big amount of atmospheric drag, that is a dissipative force that precludes using such symplectic integrators.

In that specific case (and for cases where you cannot use /do not have access to/do not want to use symplectic integrators), I would recommend the use of Bulirsch-Stoer integrator if you need precision and efficiency over long timescales. It works well by experience, and is also recommended by the Numerical Recipes (Press et al., 2007).

summarizing some points:

1. If it's a long-term integration of a non-dissapative model, a symplectic integrator is what you're looking for.
2. Otherwise, since it's an equation of motion, Runge-Kutta Nystrom methods will be more efficient than a transformation to a first order system. There are high order RKN methods due to DP. There are some implementations, like here in Julia they are documented and here's a MATLAB one.
3. High order Runge-Kutta methods are only needed if you want a high accuracy solution. If it's lower tolerances then a 5th order RK will likely be faster (for the same error). The best thing to do if you need to solve this often is to test a bunch of different methods. In this set of benchmarks on 3-body problems we see that (for the same error) high order RK methods are only a marginal improvement in speed, though as error -> 0 you can see that the improvement already goes to >5x against Dormand-Prince 45 (DP5) when you're looking at 4 digits of accuracy (tolerances are wayyyy lower for this though. Tolerances are only a ballpark in any problem). As you pull the tolerances even lower the improvement from a high order RK method grows, but you may need to start using higher precision numbers.
4. Dormand-Prince order 7/8 algorithm has a different 8th order tableau than the DP853 method of Hairer's dop853 and DifferentialEquations.jl's DP8 (which are the same). The latter 853 method cannot be implemented in standard tableau version of an Runge-Kutta method since it's error estimator is non-standard. But this method is much more efficient and I wouldn't recommend even using the older Fehlberg 7/8 or DP 7/8 methods.
5. For high order RK methods, the Verner "Efficient" methods are the gold standard. That shows up in the benchmarks I linked. You can code those into Boost yourself, or use one of the 2 packages which implement them if you want those easier (Mathematica or DifferentialEquations.jl).