I was browsing the C++ programs on "The Computer Language Shootout"
http://shootout.alioth.debian.org/ and for the heck of it, I thought
I'd try my hand at seeing if I couldn't shorten the entry for the
N-Body benchmark...
....Without too much concern for speed (admittedly the main emphasis of
the site), I whittled the C++ example down from 127 lines to 67 lines,
which would place it first in the lines-of-code ranking. (The Shootout
doesn't count lines consisting of only whitespace or solitary curly
braces). I'd be interested in seeing if anyone had additional ideas
for compressing it further without resorting to things like stripping
out whitespace and cramming everything onto a couple of unnaturally
long lines. Also keep in mind the constraints of the Shootout: you
only have access to the standard library, no Boost, etc. I've
included the code below, but you can see a syntax highlighted version
at ( http://sleepingsquirrel.org/cpp/body4.cpp.html ) along with an
alternative version which lays out the data structures in a different
manner ( http://sleepingsquirrel.org/cpp/body5.cpp.html ). Comments on
style would also be appreciated. (BTW, many other C++ programs on the
Shootout site look like they could use some TLC.)
Thanks,
Greg Buchholz
#include <iostream> /* The Great Computer Language Shootout */
#include <iomanip> /* http://shootout.alioth.debian.org/ */
#include <numeric> /* originally by:Christoph Bauer & David McCombs*/
#include <valarray>
#include <vector>
#define solar_mass (4 * M_PI * M_PI)
#define days_per_year 365.24
#define NBODIES 5
#define DIM 3
using namespace std;
typedef valarray<double > V;
template <class T> T mag(const valarray<T>& v){ return
sqrt((v*v).sum( )); }
class Planet {
public: V pos;
V vel;
double mass;
Planet () : pos(DIM), vel(DIM) {};
Planet (double m, double x, double y, double z,
double vx, double vy, double vz) :
pos(DIM),vel(DI M)
{
mass = m; pos[0] = x; pos[1] = y; pos[2] = z;
vel[0] = vx; vel[1] = vy; vel[2] = vz;
};
};
V momentum(V& m, const Planet& p){ return m + p.mass * p.vel; }
void offset_momentum (vector<Planet> & b)
{
V z(0.0,DIM);
b[0].vel = (-1/solar_mass) * accumulate(b.be gin()+1,b.end() , z,
momentum);
}
double mv_sqrd(double acc, const Planet& p)
{
return acc + p.mass * (p.vel * p.vel).sum();
}
double energy(const vector<Planet>& b)
{
double kinetic, potential = 0;
kinetic = 0.5 * accumulate(b.be gin(),b.end(), 0.0, mv_sqrd);
for(int i=0; i<NBODIES; i++)
for(int j=i+1; j<NBODIES; j++)
potential += b[i].mass * b[j].mass / mag((V)(b[i].pos -
b[j].pos));
return kinetic - potential;
}
void advance(vector< Planet>& b, double dt)
{
for(int i=0; i<NBODIES; i++)
{
for(int j=i+1; j<NBODIES; j++)
{
V d = b[i].pos - b[j].pos;
double dist = mag(d);
double mag = dt / (dist * dist * dist);
b[i].vel -= b[j].mass * mag * d;
b[j].vel += b[i].mass * mag * d;
}
b[i].pos += dt * b[i].vel;
}
}
Planet dpy(Planet p)
{
p.vel *= days_per_year;
return p;
}
int main(int argc, char *argv[])
{
Planet sun(solar_mass, 0, 0, 0, 0, 0, 0);
Planet jupiter( 9.5479193842432 6609e-04 * solar_mass,
4.8414314424647 2090e+00,-1.1603200440274 2839e+00,-1.0362204447112 3109e-01,
1.6600766427440 3694e-03,
7.6990111841974 0425e-03,-6.9046001697206 3023e-05);
Planet saturn( 2.8588598066613 0812e-04 * solar_mass,
8.3433667182445 7987e+00,
4.1247985641243 0479e+00,-4.0352341711432 1381e-01,
-2.7674251072686 2411e-03, 4.9985280123491 7238e-03,
2.3041729757376 3929e-05);
Planet uranus( 4.3662440433515 6298e-05 * solar_mass,
1.2894369562139 1310e+01,-1.5111151401698 6312e+01,-2.2330757889265 5734e-01,
2.9646013756476 1618e-03,
2.3784717395948 0950e-03,-2.9658956854023 7556e-05);
Planet neptune( 5.1513890204661 1451e-05 * solar_mass,
1.5379697114850 9165e+01,-2.5919314609987 9641e+01,
1.7925877295037 1181e-01,
2.6806777249038 9322e-03,
1.6282417003824 2295e-03,-9.5159225451971 5870e-05);
Planet tmp[] = {sun, jupiter, saturn, uranus, neptune};
vector<Planet> bodies(tmp,tmp+ NBODIES);
transform(bodie s.begin(),bodie s.end(),bodies. begin(),dpy);
offset_momentum (bodies);
cout << setprecision(9) << energy(bodies) << endl;
for(long i=atoi(argv[1]); i>0; i--) advance(bodies, 0.01);
cout << energy(bodies) << endl;
return 0;
}
http://shootout.alioth.debian.org/ and for the heck of it, I thought
I'd try my hand at seeing if I couldn't shorten the entry for the
N-Body benchmark...
....Without too much concern for speed (admittedly the main emphasis of
the site), I whittled the C++ example down from 127 lines to 67 lines,
which would place it first in the lines-of-code ranking. (The Shootout
doesn't count lines consisting of only whitespace or solitary curly
braces). I'd be interested in seeing if anyone had additional ideas
for compressing it further without resorting to things like stripping
out whitespace and cramming everything onto a couple of unnaturally
long lines. Also keep in mind the constraints of the Shootout: you
only have access to the standard library, no Boost, etc. I've
included the code below, but you can see a syntax highlighted version
at ( http://sleepingsquirrel.org/cpp/body4.cpp.html ) along with an
alternative version which lays out the data structures in a different
manner ( http://sleepingsquirrel.org/cpp/body5.cpp.html ). Comments on
style would also be appreciated. (BTW, many other C++ programs on the
Shootout site look like they could use some TLC.)
Thanks,
Greg Buchholz
#include <iostream> /* The Great Computer Language Shootout */
#include <iomanip> /* http://shootout.alioth.debian.org/ */
#include <numeric> /* originally by:Christoph Bauer & David McCombs*/
#include <valarray>
#include <vector>
#define solar_mass (4 * M_PI * M_PI)
#define days_per_year 365.24
#define NBODIES 5
#define DIM 3
using namespace std;
typedef valarray<double > V;
template <class T> T mag(const valarray<T>& v){ return
sqrt((v*v).sum( )); }
class Planet {
public: V pos;
V vel;
double mass;
Planet () : pos(DIM), vel(DIM) {};
Planet (double m, double x, double y, double z,
double vx, double vy, double vz) :
pos(DIM),vel(DI M)
{
mass = m; pos[0] = x; pos[1] = y; pos[2] = z;
vel[0] = vx; vel[1] = vy; vel[2] = vz;
};
};
V momentum(V& m, const Planet& p){ return m + p.mass * p.vel; }
void offset_momentum (vector<Planet> & b)
{
V z(0.0,DIM);
b[0].vel = (-1/solar_mass) * accumulate(b.be gin()+1,b.end() , z,
momentum);
}
double mv_sqrd(double acc, const Planet& p)
{
return acc + p.mass * (p.vel * p.vel).sum();
}
double energy(const vector<Planet>& b)
{
double kinetic, potential = 0;
kinetic = 0.5 * accumulate(b.be gin(),b.end(), 0.0, mv_sqrd);
for(int i=0; i<NBODIES; i++)
for(int j=i+1; j<NBODIES; j++)
potential += b[i].mass * b[j].mass / mag((V)(b[i].pos -
b[j].pos));
return kinetic - potential;
}
void advance(vector< Planet>& b, double dt)
{
for(int i=0; i<NBODIES; i++)
{
for(int j=i+1; j<NBODIES; j++)
{
V d = b[i].pos - b[j].pos;
double dist = mag(d);
double mag = dt / (dist * dist * dist);
b[i].vel -= b[j].mass * mag * d;
b[j].vel += b[i].mass * mag * d;
}
b[i].pos += dt * b[i].vel;
}
}
Planet dpy(Planet p)
{
p.vel *= days_per_year;
return p;
}
int main(int argc, char *argv[])
{
Planet sun(solar_mass, 0, 0, 0, 0, 0, 0);
Planet jupiter( 9.5479193842432 6609e-04 * solar_mass,
4.8414314424647 2090e+00,-1.1603200440274 2839e+00,-1.0362204447112 3109e-01,
1.6600766427440 3694e-03,
7.6990111841974 0425e-03,-6.9046001697206 3023e-05);
Planet saturn( 2.8588598066613 0812e-04 * solar_mass,
8.3433667182445 7987e+00,
4.1247985641243 0479e+00,-4.0352341711432 1381e-01,
-2.7674251072686 2411e-03, 4.9985280123491 7238e-03,
2.3041729757376 3929e-05);
Planet uranus( 4.3662440433515 6298e-05 * solar_mass,
1.2894369562139 1310e+01,-1.5111151401698 6312e+01,-2.2330757889265 5734e-01,
2.9646013756476 1618e-03,
2.3784717395948 0950e-03,-2.9658956854023 7556e-05);
Planet neptune( 5.1513890204661 1451e-05 * solar_mass,
1.5379697114850 9165e+01,-2.5919314609987 9641e+01,
1.7925877295037 1181e-01,
2.6806777249038 9322e-03,
1.6282417003824 2295e-03,-9.5159225451971 5870e-05);
Planet tmp[] = {sun, jupiter, saturn, uranus, neptune};
vector<Planet> bodies(tmp,tmp+ NBODIES);
transform(bodie s.begin(),bodie s.end(),bodies. begin(),dpy);
offset_momentum (bodies);
cout << setprecision(9) << energy(bodies) << endl;
for(long i=atoi(argv[1]); i>0; i--) advance(bodies, 0.01);
cout << energy(bodies) << endl;
return 0;
}
Comment