import std.stdio: writefln;
import std.math: sqrt;
import std.traits: Unqual;
import std.typetuple: allSatisfy;
import core.simd: double2, float4;
import core.runtime: Runtime;

extern(C) pure nothrow int atoi(in char*); // For LDC.

enum isDouble(T) = is(Unqual!T == double);

version (LDC) {
    pragma(LDC_intrinsic, "llvm.x86.sse2.cvtpd2ps") float4 cvtpd2ps(double2) pure nothrow;
    pragma(LDC_intrinsic, "llvm.x86.sse2.cvtps2pd") double2 cvtps2pd(float4) pure nothrow;
    pragma(LDC_intrinsic, "llvm.x86.sse.rsqrt.ps") float4 rsqrtps(float4) pure nothrow;
} else /*DMD*/ {
    void16 cvtpd2ps(void16 v) pure nothrow { return __simd(XMM.CVTPD2PS, v); }
    void16 cvtps2pd(void16 v) pure nothrow { return __simd(XMM.CVTPS2PD, v); }
    void16 rsqrtps(void16 v) pure nothrow { return __simd(XMM.RSQRTPS, v); }
}

enum double PI = 3.141592653589793;
enum double solarMass = 4 * PI * PI;
enum double daysPerYear = 365.24;

struct Body {
    double[3] x;
    double filler_;
    double[3] v;
    double mass;

    this(in double x0, in double x1, in double x2,
         in double v0, in double v1, in double v2,
         in double mass_) pure nothrow {
        this.x = [x0, x1, x2];
        this.filler_ = 0;
        this.v = [v0, v1, v2];
        this.mass = mass_;
    }
}

struct NBodySystem {
private:
    __gshared static Body[5] bodies = [
        // Sun.
        Body(0., 0., 0.,
             0., 0., 0.,
             solarMass),

        // Jupiter.
        Body( 4.84143144246472090e+00,
             -1.16032004402742839e+00,
             -1.03622044471123109e-01,
              1.66007664274403694e-03 * daysPerYear,
              7.69901118419740425e-03 * daysPerYear,
             -6.90460016972063023e-05 * daysPerYear,
              9.54791938424326609e-04 * solarMass),

        // Saturn.
        Body( 8.34336671824457987e+00,
              4.12479856412430479e+00,
             -4.03523417114321381e-01,
             -2.76742510726862411e-03 * daysPerYear,
              4.99852801234917238e-03 * daysPerYear,
              2.30417297573763929e-05 * daysPerYear,
              2.85885980666130812e-04 * solarMass),

        // Uranus.
        Body( 1.28943695621391310e+01,
             -1.51111514016986312e+01,
             -2.23307578892655734e-01,
              2.96460137564761618e-03 * daysPerYear,
              2.37847173959480950e-03 * daysPerYear,
             -2.96589568540237556e-05 * daysPerYear,
              4.36624404335156298e-05 * solarMass),

        // Neptune.
        Body( 1.53796971148509165e+01,
             -2.59193146099879641e+01,
              1.79258772950371181e-01,
              2.68067772490389322e-03 * daysPerYear,
              1.62824170038242295e-03 * daysPerYear,
             -9.51592254519715870e-05 * daysPerYear,
              5.15138902046611451e-05 * solarMass)
    ];

    void offset_momentum() nothrow {
        foreach (const ref b; bodies)
            foreach (immutable i; 0 .. 3)
                bodies[0].v[i] -= b.v[i] * b.mass / solarMass;
    }

public:

    void initialize() {
        offset_momentum;
    }

    void advance(in double dt) nothrow {
        enum N = ((bodies.length - 1) * bodies.length) / 2;
        static double[3][N] r = void;
        static double[N] mag = void;

        {
            size_t i = 0;
            foreach (immutable j, ref const bi; bodies) {
                foreach (const ref bj; bodies[j + 1 .. $]) {
                    foreach (immutable m; 0 .. 3)
                        r[i][m] = bi.x[m] - bj.x[m];
                    i++;
                }
            }
        }

        for (size_t i = 0; i < N; i += 2) {
            double2[3] dx = [[r[i][0], r[i + 1][0]],
                             [r[i][1], r[i + 1][1]],
                             [r[i][2], r[i + 1][2]]];

            /*immutable*/ double2 dsquared = dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2];
            double2 distance = cvtps2pd(rsqrtps(cvtpd2ps(dsquared)));

            foreach (immutable _; 0 .. 2)
                distance = distance * 1.5 - ((0.5 * dsquared) * distance) * (distance * distance);

            /*immutable*/ double2 dt2 = [dt, dt];
            /*immutable*/ double2 dmag = dt2 / dsquared * distance;
            mag[i    ] = dmag.array[0];
            mag[i + 1] = dmag.array[1];
        }

        {
            size_t i = 0;
            foreach (immutable j, ref bi; bodies) {
                foreach (ref bj; bodies[j + 1 .. $]) {
                    foreach (immutable m; 0 .. 3) {
                        bi.v[m] -= r[i][m] * bj.mass*mag[i];
                        bj.v[m] += r[i][m] * bi.mass*mag[i];
                    }
                    i++;
                }
            }
        }

        foreach (ref bdy; bodies)
            foreach (immutable m; 0 .. 3)
                bdy.x[m] += dt * bdy.v[m];
    }

    @property double energy() const nothrow {
        double e = 0.0;

        foreach (immutable i, ref const bi; bodies) {
            e += bi.mass * (bi.v[0] * bi.v[0] +
                            bi.v[1] * bi.v[1] +
                            bi.v[2] * bi.v[2]) / 2.0;

            foreach (const ref bj; bodies[i + 1 .. $]) {
                double[3] dx = void;
                foreach (immutable m; 0 .. 3) {
                    dx[m] = bi.x[m] - bj.x[m];
                    dx[m] *= dx[m];
                }

                immutable double distance = sqrt(dx[0] + dx[1] + dx[2]);
                e -= (bi.mass * bj.mass) / distance;
            }
        }

        return e;
    }
}

void main() {
    immutable uint n = atoi(Runtime.cArgs.argv[1]);
    NBodySystem system;
    system.initialize;

    writefln("%.9f", system.energy);
    foreach (immutable _; 0 .. n)
        system.advance(0.01);
    writefln("%.9f", system.energy);
}