GCCBench.c

#include <stdint.h>
#include <stdalign.h>
#include <immintrin.h>
#include <math.h>
#include <stdio.h>

// intptr_t should be the native integer type on most sane systems.
typedef intptr_t intnative_t;

typedef struct{
    double position[3], velocity[3], mass;
} body;

#define SOLAR_MASS (4*M_PI*M_PI)
#define DAYS_PER_YEAR 365.24
#define BODIES_COUNT 5

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


// Advance all the bodies in the system by one timestep. Calculate the
// interactions between all the bodies, update each body's velocity based on
// those interactions, and update each body's position by the distance it
// travels in a timestep at it's updated velocity.
static void advance(body bodies[]){

    // Figure out how many total different interactions there are between each
    // body and every other body. Some of the calculations for these
    // interactions will be calculated two at a time by using x86 SSE
    // instructions and because of that it will also be useful to have a
    // ROUNDED_INTERACTIONS_COUNT that is equal to the next highest even number
    // which is equal to or greater than INTERACTIONS_COUNT.
    #define INTERACTIONS_COUNT (BODIES_COUNT*(BODIES_COUNT-1)/2)
    #define ROUNDED_INTERACTIONS_COUNT (INTERACTIONS_COUNT+INTERACTIONS_COUNT%2)

    // It's useful to have two arrays to keep track of the position_Deltas
    // and magnitudes of force between the bodies for each interaction. For the
    // position_Deltas array, instead of using a one dimensional array of
    // structures that each contain the X, Y, and Z components for a position
    // delta, a two dimensional array is used instead which consists of three
    // arrays that each contain all of the X, Y, and Z components for all of the
    // position_Deltas. This allows for more efficient loading of this data into
    // SSE registers. Both of these arrays are also set to contain
    // ROUNDED_INTERACTIONS_COUNT elements to simplify one of the following
    // loops and to also keep the second and third arrays in position_Deltas
    // aligned properly.
    static alignas(__m128d) double
      position_Deltas[3][ROUNDED_INTERACTIONS_COUNT],
      magnitudes[ROUNDED_INTERACTIONS_COUNT];

    // Calculate the position_Deltas between the bodies for each interaction.
    for(intnative_t i=0, k=0; i<BODIES_COUNT-1; ++i)
        for(intnative_t j=i+1; j<BODIES_COUNT; ++j, ++k)
            for(intnative_t m=0; m<3; ++m)
                position_Deltas[m][k]=
                  bodies[i].position[m]-bodies[j].position[m];

    // Calculate the magnitudes of force between the bodies for each
    // interaction. This loop processes two interactions at a time which is why
    // ROUNDED_INTERACTIONS_COUNT/2 iterations are done.
    for(intnative_t i=0; i<ROUNDED_INTERACTIONS_COUNT/2; ++i){

        // Load position_Deltas of two bodies into position_Delta.
        __m128d position_Delta[3];
        for(intnative_t m=0; m<3; ++m)
            position_Delta[m]=((__m128d *)position_Deltas[m])[i];

        const __m128d distance_Squared=
          position_Delta[0]*position_Delta[0]+
          position_Delta[1]*position_Delta[1]+
          position_Delta[2]*position_Delta[2];

        // Doing square roots normally using double precision floating point
        // math can be quite time consuming so SSE's much faster single
        // precision reciprocal square root approximation instruction is used as
        // a starting point instead. The precision isn't quite sufficient to get
        // acceptable results so two iterations of the Newton–Raphson method are
        // done to improve precision further.
        __m128d distance_Reciprocal=
          _mm_cvtps_pd(_mm_rsqrt_ps(_mm_cvtpd_ps(distance_Squared)));
        for(intnative_t j=0; j<2; ++j)
            // Normally the last four multiplications in this equation would
            // have to be done sequentially but by placing the last
            // multiplication in parentheses, a compiler can then schedule that
            // multiplication earlier.
            distance_Reciprocal=distance_Reciprocal*1.5-
              0.5*distance_Squared*distance_Reciprocal*
              (distance_Reciprocal*distance_Reciprocal);

        // Calculate the magnitudes of force between the bodies. Typically this
        // calculation would probably be done by using a division by the cube of
        // the distance (or similarly a multiplication by the cube of its
        // reciprocal) but for better performance on modern computers it often
        // will make sense to do part of the calculation using a division by the
        // distance_Squared which was already calculated earlier. Additionally
        // this method is probably a little more accurate due to less rounding
        // as well.
        ((__m128d *)magnitudes)[i]=0.01/distance_Squared*distance_Reciprocal;
    }

    // Use the calculated magnitudes of force to update the velocities for all
    // of the bodies.
    for(intnative_t i=0, k=0; i<BODIES_COUNT-1; ++i)
        for(intnative_t j=i+1; j<BODIES_COUNT; ++j, ++k){
            // Precompute the products of the mass and magnitude since it can be
            // reused a couple times.
            const double
              i_mass_magnitude=bodies[i].mass*magnitudes[k],
              j_mass_magnitude=bodies[j].mass*magnitudes[k];
            for(intnative_t m=0; m<3; ++m){
                bodies[i].velocity[m]-=position_Deltas[m][k]*j_mass_magnitude;
                bodies[j].velocity[m]+=position_Deltas[m][k]*i_mass_magnitude;
            }
        }

    // Use the updated velocities to update the positions for all of the bodies.
    for(intnative_t i=0; i<BODIES_COUNT; ++i)
        for(intnative_t m=0; m<3; ++m)
            bodies[i].position[m]+=0.01*bodies[i].velocity[m];
}


// Calculate the momentum of each body and conserve momentum of the system by
// adding to the Sun's velocity the appropriate opposite velocity needed in
// order to offset that body's momentum.
static void offset_Momentum(body bodies[]){
    for(intnative_t i=0; i<BODIES_COUNT; ++i)
        for(intnative_t m=0; m<3; ++m)
            bodies[0].velocity[m]-=
              bodies[i].velocity[m]*bodies[i].mass/SOLAR_MASS;
}


// Output the total energy of the system.
static void output_Energy(body bodies[]){
    double energy=0;
    for(intnative_t i=0; i<BODIES_COUNT; ++i){

        // Add the kinetic energy for each body.
        energy+=0.5*bodies[i].mass*(
          bodies[i].velocity[0]*bodies[i].velocity[0]+
          bodies[i].velocity[1]*bodies[i].velocity[1]+
          bodies[i].velocity[2]*bodies[i].velocity[2]);

        // Add the potential energy between this body and every other body.
        for(intnative_t j=i+1; j<BODIES_COUNT; ++j){
            double position_Delta[3];
            for(intnative_t m=0; m<3; ++m)
                position_Delta[m]=bodies[i].position[m]-bodies[j].position[m];

            energy-=bodies[i].mass*bodies[j].mass/sqrt(
              position_Delta[0]*position_Delta[0]+
              position_Delta[1]*position_Delta[1]+
              position_Delta[2]*position_Delta[2]);
        }
    }

    // Output the total energy of the system.
    printf("%.9f\n", energy);
}


int main(int argc, char *argv[]){
    offset_Momentum(solar_Bodies);
    output_Energy(solar_Bodies);
    for(intnative_t n=atoi(argv[1]); n--; advance(solar_Bodies));
    output_Energy(solar_Bodies);
}