The Computer Language
22.05 Benchmarks Game

n-body C++ g++ #8 program

source code

/* The Computer Language Benchmarks Game
   https://salsa.debian.org/benchmarksgame-team/benchmarksgame/

   An implementation pretty much from scratch, with inspiration from the Rust
   version, which used the idea of saving some of the ingredients of the
   compution in an array instead of recomputing them.
   
   contributed by cvergu
*/

#include <cstdio>
#include <cmath>

constexpr double SOLAR_MASS = 4 * M_PI * M_PI;
constexpr double DAYS_PER_YEAR = 365.24;
constexpr std::size_t BODIES_COUNT = 5;

struct alignas(32) vector3d {
    double x, y, z;

    constexpr double norm() const noexcept
    {
        return x*x + y*y + z*z;
    }

#pragma GCC optimize ("no-math-errno")
    double magnitude(double dt) const noexcept
    {
        double sum = this->norm();
        return dt / (sum * sqrt(sum));
    }
};

constexpr vector3d operator+(vector3d v1, vector3d v2)
{
    return vector3d {v1.x + v2.x, v1.y + v2.y, v1.z + v2.z};
}


constexpr vector3d operator-(vector3d v1, vector3d v2)
{
    return vector3d {v1.x - v2.x, v1.y - v2.y, v1.z - v2.z};
}

constexpr vector3d& operator+=(vector3d& v1, vector3d v2)
{
    v1.x += v2.x;
    v1.y += v2.y;
    v1.z += v2.z;

    return v1;
}

constexpr vector3d& operator-=(vector3d& v1, vector3d v2)
{
    v1.x -= v2.x;
    v1.y -= v2.y;
    v1.z -= v2.z;

    return v1;
}

constexpr vector3d& operator*=(vector3d& v, double mag)
{
    v.x *= mag;
    v.y *= mag;
    v.z *= mag;

    return v;
}

constexpr vector3d operator*(vector3d v, double mag)
{
    return vector3d {v.x * mag, v.y * mag, v.z * mag};
}

constexpr vector3d operator/(vector3d v, double mag)
{
    return vector3d {v.x / mag, v.y / mag, v.z / mag};
}

struct body {
    vector3d position;
    vector3d velocity;
    double mass;
};

void advance(body state[BODIES_COUNT], double dt)
{
    /*
     * We precompute the quantity (r_i - r_j)
     */
    vector3d rij[BODIES_COUNT][BODIES_COUNT];

    for (std::size_t i = 0; i < BODIES_COUNT; ++i)
        for (std::size_t j = i + 1; j < BODIES_COUNT; ++j)
            rij[i][j] = state[i].position - state[j].position;

    double magnitudes[BODIES_COUNT][BODIES_COUNT];

    for (std::size_t i = 0; i < BODIES_COUNT; ++i)
        for (std::size_t j = i + 1; j < BODIES_COUNT; ++j)
            magnitudes[i][j] = rij[i][j].magnitude(dt);

    /*
     * Compute the new speed using v_i = a_i dt, where
     * a_i = \sum_{j \neq i} m_j (r_i - r_j)/|r_i - r_j|^3
     */
    for (std::size_t i = 0; i < BODIES_COUNT; ++i) {
        for (std::size_t j = i + 1; j < BODIES_COUNT; ++j) {
            vector3d dist = rij[i][j];
            double mag = magnitudes[i][j];
            state[i].velocity -= dist * (state[j].mass * mag);
            state[j].velocity += dist * (state[i].mass * mag);
        }
    }

    /*
     * Compute the new positions
     */
    for (std::size_t i = 0; i < BODIES_COUNT; ++i)
        state[i].position += state[i].velocity * dt;
}

void offset_momentum(body state[BODIES_COUNT])
{
    vector3d& sun_velocity = state[0].velocity;

    for (std::size_t i = 1; i < BODIES_COUNT; ++i)
        sun_velocity -= state[i].velocity * state[i].mass / SOLAR_MASS;
}

double energy(const body state[BODIES_COUNT])
{
    double energy = 0;

    for (std::size_t i = 0; i < BODIES_COUNT; ++i) {
        const body& body1 = state[i];
        energy += 0.5 * body1.mass * body1.velocity.norm();
        #pragma clang loop vectorize(enable)
        for (std::size_t j = i + 1; j < BODIES_COUNT; ++j) {
            const body& body2 = state[j];
            vector3d r12 = body1.position - body2.position;
            energy -= body1.mass * body2.mass / sqrt(r12.norm());
        }
    }

    return energy;
}


body state[] = {
        // Sun
        {.position = {0, 0, 0},
         .velocity = {0, 0, 0},
         .mass = SOLAR_MASS},
        // Jupiter
        {.position = {4.84143144246472090e+00,
                      -1.16032004402742839e+00,
                      -1.03622044471123109e-01},
         .velocity = {1.66007664274403694e-03 * DAYS_PER_YEAR,
                      7.69901118419740425e-03 * DAYS_PER_YEAR,
                     -6.90460016972063023e-05 * DAYS_PER_YEAR},
         .mass = 9.54791938424326609e-04 * SOLAR_MASS},
        // Saturn
        {.position = {8.34336671824457987e+00,
                      4.12479856412430479e+00,
                      -4.03523417114321381e-01},
         .velocity = {-2.76742510726862411e-03 * DAYS_PER_YEAR,
                      4.99852801234917238e-03 * DAYS_PER_YEAR,
                      2.30417297573763929e-05 * DAYS_PER_YEAR},
         .mass = 2.85885980666130812e-04 * SOLAR_MASS},
        // Uranus
        {.position = {1.28943695621391310e+01,
                      -1.51111514016986312e+01,
                      -2.23307578892655734e-01},
         .velocity = {2.96460137564761618e-03 * DAYS_PER_YEAR,
                      2.37847173959480950e-03 * DAYS_PER_YEAR,
                      -2.96589568540237556e-05 * DAYS_PER_YEAR},
         .mass = 4.36624404335156298e-05 * SOLAR_MASS},
        // Neptune
        {.position = {1.53796971148509165e+01,
                      -2.59193146099879641e+01,
                      1.79258772950371181e-01},
         .velocity = {2.68067772490389322e-03 * DAYS_PER_YEAR,
                      1.62824170038242295e-03 * DAYS_PER_YEAR,
                      -9.51592254519715870e-05 * DAYS_PER_YEAR},
         .mass = 5.15138902046611451e-05 * SOLAR_MASS}
};


int main(int argc, char** argv) {
    if (argc < 1) return EXIT_FAILURE;

    int n = atoi(argv[1]);

    offset_momentum(state);

    printf("%.9f\n", energy(state));

    for (int i = 0; i < n; ++i)
        advance(state, 0.01);
    printf("%.9f\n", energy(state));
}
    

notes, command-line, and program output

NOTES:
64-bit Ubuntu quad core
11.2.0-19ubuntu1


Mon, 02 May 2022 20:52:52 GMT

MAKE:
/usr/bin/g++ -c -pipe -O3 -fomit-frame-pointer -march=ivybridge  -std=c++2a nbody.gpp-8.c++ -o nbody.gpp-8.c++.o &&  \
        /usr/bin/g++ nbody.gpp-8.c++.o -o nbody.gpp-8.gpp_run  
rm nbody.gpp-8.c++

2.72s to complete and log all make actions

COMMAND LINE:
./nbody.gpp-8.gpp_run 50000000

PROGRAM OUTPUT:
-0.169075164
-0.169059907