source code
/*
* The Computer Language Benchmarks Game
* https://salsa.debian.org/benchmarksgame-team/benchmarksgame/
*
* contributed by Jeffrey Beu
*
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
namespace
{
template< typename ... > struct tuple{};
template<typename H, typename... T>
struct tuple<H,T...>
{
typedef H head;
typedef tuple<T...> tail;
head fst;
tail snd;
tuple(const head& a,const tail& b):fst(a),snd(b){}
tuple() = delete;
};
tuple<> construct_tuple()
{ return tuple<>();
}
template<class H, class...T> tuple<H,T...> construct_tuple(H h, T... a)
{
tuple<H,T...> r(h,construct_tuple(a...));
return r;
}
template<class OP> void FOREACH(OP,tuple<>){}
template<typename OP,class ...T> void FOREACH(OP& op,tuple<T...>& t)
{
op(t.fst);
FOREACH(op,t.snd);
}
template<class OP> void FOR(OP,tuple<>){}
template<typename OP, class ...T> void FOR(OP& op,tuple<T...>& t)
{
op(t);
FOR(op,t.snd);
}
struct body {
double x, y, z;
double vx, vy, vz;
const double mass;
};
const double pi = 3.141592653589793;
const double solar_mass = 4 * pi * pi;
const double days_per_year = 365.24;
const double dt = 0.01;
const body sun =
{
0,0,0,
0,0,0,
solar_mass
};
const body 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
};
const body 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
};
const body 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
};
const body 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
};
struct
{ struct
{ template<class T> struct foreach_block
{ T &t;
template<class T2> void operator()(T2& b2)
{
auto& b = t.fst;
double dx = b.x - b2.x;
double dy = b.y - b2.y;
double dz = b.z - b2.z;
double distance = sqrt(dx * dx + dy * dy + dz * dz);
double mag = dt / (distance * distance * distance);
b.vx -= dx * b2.mass * mag;
b.vy -= dy * b2.mass * mag;
b.vz -= dz * b2.mass * mag;
b2.vx += dx * b.mass * mag;
b2.vy += dy * b.mass * mag;
b2.vz += dz * b.mass * mag;
}
};
template<class T> void operator()(T& bodies)
{
foreach_block<decltype(bodies)> run = {bodies};
FOREACH(run,bodies.snd);
}
}for_block;
struct
{
template<class T> void operator()(T& body)
{
body.x += dt * body.vx;
body.y += dt * body.vy;
body.z += dt * body.vz;
}
}foreach_block;
template<class T> void operator()(T& bodies)
{
FOR(for_block,bodies);
FOREACH(foreach_block,bodies);
}
}advance;
struct
{
struct for_block
{ double &e;
template<class T> struct foreach_block
{ T &t;
double &e;
template<class T2> void operator()(T2& b2)
{
auto& b = t.fst;
double dx = b.x - b2.x;
double dy = b.y - b2.y;
double dz = b.z - b2.z;
double distance = sqrt(dx * dx + dy * dy + dz * dz);
e -= (b.mass * b2.mass) / distance;
}
};
template<class T> void operator()(T& bodies)
{
auto& p = bodies.fst;
e += 0.5 * p.mass * (p.vx * p.vx + p.vy * p.vy + p.vz * p.vz);
foreach_block<decltype(bodies)> run = {bodies,e};
FOREACH(run,bodies.snd);
}
};
template<class T> double operator()(T& bodies)
{
double e = 0.0;
for_block run = {e};
FOR(run,bodies);
return e;
}
}energy;
struct
{
struct foreach_block
{
double &px, &py, &pz;
template<class T> void operator()(T& body)
{
px += body.vx * body.mass;
py += body.vy * body.mass;
pz += body.vz * body.mass;
}
};
template<class T> void operator()(T& bodies)
{
double px = 0.0, py = 0.0, pz = 0.0;
foreach_block run = {px,py,pz};
FOREACH(run,bodies);
bodies.fst.vx = -px / solar_mass;
bodies.fst.vy = -py / solar_mass;
bodies.fst.vz = -pz / solar_mass;
}
}offset;
}
int main(int,char** argv)
{
auto solar_system = construct_tuple(sun,jupiter,saturn,uranus,neptune);
offset(solar_system);
printf ("%.9f\n", energy(solar_system));
int n = atoi(argv[1]);
for (int i = 1; i <= n; i++)
{
advance(solar_system);
}
printf ("%.9f\n", energy(solar_system));
return 0;
}