I want to find optimizations to my code for the 3BP, and more specifically computing accelerations. I'm using a data-driven approach, so I have a bodies structure that contains arrays for mass, position, velocity, and acceleration, all stored as Vector2d's of the library Eigen. So something like array<Vector2d, 3> position. (Mass is obviously a double array though).

What I have now is as follows:

void updateAccelerations() {
    //The rHat/r^2 part of the forces vectors, which is mutual between pairs of bodies
    Vector2d mutualVector;

    for (Vector2d &acc : bodies.accList)
    //i, j just loop through unordered pairs of {0,1,2}
    int j;
    for (int i = 0; i < NUM; i++) {
        j = (i+1)%NUM;
        mutualVector = directedInverseSquare(i, j);
        bodies.accList[i] += bodies.massList[j] * G * mutualVector;
        bodies.accList[j] += - bodies.massList[i] * G * mutualVector;

Vector2d directedInverseSquare(int i, int j) {

    Vector2d diff = bodies.posList[j] - bodies.posList[i];
    Vector2d rHat = diff.normalized();
    double rSquared = diff.squaredNorm();

    return rHat / rSquared;

Some things I've tried are pushing the setZero command to an existing loop in my code, and multiplying the mutualVector by G instead of multiplying two times when computing the accelerations. To my surprise, both changes slowed my code down by a slight margin.

An idea that I had, that is quite nice mathematically but probably not practical, is to recognize that:

$$G\cdot \begin{pmatrix} m_1 \\\ m_2 \\\ m_3 \end{pmatrix} \times \begin{pmatrix} \vec{R}_{23} \\\ \vec{R}_{31} \\\ \vec{R}_{12} \end{pmatrix} = \begin{pmatrix} \vec{a}_1 \\\ \vec{a}_2 \\\ \vec{a}_3 \end{pmatrix}$$

Where $\vec{R}_{ij} = \frac{\hat{r}_{ij}}{\lVert\vec{r}_{ij}\rVert ^2}$

This is probably useless, especially since Eigen doesn't allow a cross product between a vector of doubles and of Vectors. I thought I'd mention it though.

If anyone can think of an improvement to the algorithm (or some optimization to the code implementation/other things), I'd be interested in hearing them.

Edit: Someone wanted me to post my Bodies struct:

using namespace Eigen;

#define NUM 3

typedef vector<tuple<double, Vector2d, Vector2d>> initialData;
typedef array<Vector2d, NUM> vData;
typedef array<double, NUM> Data;

struct Bodies {

    vData posList;
    vData velList;
    vData accList;
    Data massList;

    Bodies(initialData init) {

        if (init.size() != NUM)
            throw invalid_argument("Not Correct Body Amount");

        int i = 0;
        for (auto const& [mass, pos, vel] : init) {
            posList[i] = pos;
            velList[i] = vel;
            accList[i] = Vector2d(0,0);
            massList[i] = mass;

    //functions for calculating energy conservation (irrelevant)

Keep in mind that the initialization only happens once, so I really don't care for its efficiency.

  • $\begingroup$ could you provide your bodies struct? $\endgroup$
    – MPIchael
    Commented Jan 17 at 16:02
  • $\begingroup$ @MPIchael Added it as an edit $\endgroup$
    – Remeraze
    Commented Jan 17 at 16:32
  • $\begingroup$ What integration method are you using? Boost::odeint, some other general library, or a home-cooked fixed-step method? $\endgroup$ Commented Jan 17 at 18:58
  • 1
    $\begingroup$ With more than 2 bodies, there is no lower bound for the distance between the bodies. This means that one can get very high speeds with very small turning radii, requiring a very small step size, and no guarantee that it is small enough. As the symplectic methods only work as advertised for fixed step size, you will get also very many points on the boring slow, low curvature segments, which is wasteful. Methods with variable step size may not preserve energy, but give similar accuracy with less effort. $\endgroup$ Commented Jan 17 at 22:45
  • 1
    $\begingroup$ Adaptive time steppers are implemented in all of the major ODE solver packages, but are quite difficult to implement yourself. You'll be well served not rolling your own ODE integrator. $\endgroup$ Commented Jan 18 at 19:24

1 Answer 1


Welcome to Scicomp!

For all following points, always measure and compare.

  • Don't fall for premature optimization! It is the root of all evil. Is the time you spend on optimizing your code worth your time compared to your other constraints (project time / available cpu resources / deadlines)? It is always possible to spend month of work to get your code to execute 10% faster...
  • Have you turned on compiler optimizations? (gcc yourprog.cc -O3 -ffastmath ...)
  • It is possible your available hardware supplies optimizations which the compiler can use. provide the compiler with the information for which target you are compiling: -march=native
  • constify your code (tell the compiler all about which values are const after initialization, which arguments will not change during method invocation and which values are compile-time-constant (constexpr))
  • try to add inline keywords to functions as suggestions to the compiler
  • can you eliminate the many multiplications of G? You could just multiply G at the very end, or combine it with the mass in a new constant outside of the loop.
  • you might save time when squaring your norm (are you squaring two times, or are you calculating ^4 in one step and leave the details to the compiler)
  • It might be possible to use SIMD instructions of your cpu. Essentially you may execute arithmetic operations on small arrays of numbers in the same time it usually takes for one single value.
  • If you are into that kind of thing, you can inspect and compare your bytecode in one of the available online compilers like: https://godbolt.org/ for the different variants you are trying.
  • It might be worth a try to unroll your for loop for the three bodies. That makes the code slightly less consise but more predictable (less jumps in reading memory into caches). You may do so manually, or tell the compiler to do it: for one loop: "attribute((optimize("unroll-loops")))"" or you can do it globally for your code: " -funroll-loops"


  • If you want to run statistics on the trajectories of the three bodies, you can consider to simulate them all at once, i.e., increase the number of the TBP systems you solve significantly. If you propagate 10k of them in one run, you should be able to use all the cores of your cpu to their full potential and still wind up with the same statistical properties. (A more primite way to do that would be to start multiple processes.) I would presume your calculation is actually cpu bound and not memory bound. It should be possible to add another loop for the number of systems you want to simulate and surround it with a "#pragma omp parallel for". A neat advantage of that approach is that you can calculate and save your measurements as you go. That might sound like a small achievement but you can inspect and stop your simulation as soon as you see some error or abnormality in the results.(https://de.wikipedia.org/wiki/OpenMP)
  • This is a long shot but try if the compiler optimizations which aim to reduce the executable size help you ("-Os" iirc). The rationale is that you can sometimes reduce the instruction set to a size that fits into the instruction set cache of the cpu. It sounds counterintuitive but I had collegues who reported speedups better than (-O3), which tends to increase the size of the executable.


  • Have you run your program with a sampling profiler yet? The idea is to run a long running simulation while sampling the code statistically as to which line or method it is currently invokating. Its a bit of tedious grunt work to get the profilers running, but once you have results you can directly inspect where the cpu is spending its time and any inefficiency will pop up at the top of the list. I used gprof in the past with good results, but there might be more modern tools.

Also please share your resulting speedups:-)

  • $\begingroup$ I've indeed added -O3 -march=native, but -ffast-math both didn't give me big performance gains, and is less accurate than the default (accuracy is also very important here, including conservation of energy and such). As for const, I read somewhere on SO that it doesn't matter for performance since there's a way to bypass it in code anyway (maybe the mutable keyword? can't remember). Is that false? I've actually tried to combine G with mass, and it weirdly gave me worse preformence! but That was a while ago, maybe I'll try again. As for the norm, I'm not sure what you mean. rephrase? $\endgroup$
    – Remeraze
    Commented Jan 17 at 16:38
  • $\begingroup$ @Remeraze The normalized and squaredNorm functions both compute the square sum of diff. This duplication can be reduced by combining the operations to diff/pow(diff.squaredNorm(),1.5) $\endgroup$ Commented Jan 17 at 17:44
  • $\begingroup$ @LutzLehmann I see. I implemented what you suggested and compared to the old one (ran a loop of a billion calls to both functions), and It's about 25 times slower . I wonder what kind of magic optimizations Eigen does so that using it is so much faster than regular pow. Also, I just want to mention that you've previously helped me on another post on SO, and it was very helpful, so thank you. $\endgroup$
    – Remeraze
    Commented Jan 17 at 18:06
  • $\begingroup$ @Remeraze : Speculation: It is possible that the compiler/optimizer inlines the code of the vector methods and recognizes the duplication of the square sum calculation, so that in the end sqrt(r2)*r2 is computed, with only once computing the square sum r2. This in general will be faster than pow(r2,1.5). In modern libraries, the power function pow(x,y) is a heavy overhead of special cases and argument reduction over the simple but at times insufficiently accurate exp(y*log(x)). $\endgroup$ Commented Jan 17 at 18:54
  • $\begingroup$ Interestingly enough, when I use sqrt(rSquared)*rSquared instead of pow, it's almost as fast (only about 10% slower) than Eigen. $\endgroup$
    – Remeraze
    Commented Jan 17 at 21:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.