Tightening up Shootout code

[Greg Buchholz]

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...


http://shootout.alioth.debian.org/debian/benchmark.php?test=nbody&lang=all


....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(DIM)
{
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.begin()+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.begin(),b.end(), 0.0, mv_sqrd);

for(int i=0; i<NBODIES; i++)
for(int j=i+1; j<NBODIES; j++)
potential += b.mass * b[j].mass / mag((V)(b.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.pos - b[j].pos;
double dist = mag(d);
double mag = dt / (dist * dist * dist);
b.vel -= b[j].mass * mag * d;
b[j].vel += b.mass * mag * d;
}
b.pos += dt * b.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.54791938424326609e-04 * solar_mass,

4.84143144246472090e+00,-1.16032004402742839e+00,-1.03622044471123109e-01,
1.66007664274403694e-03,
7.69901118419740425e-03,-6.90460016972063023e-05);
Planet saturn( 2.85885980666130812e-04 * solar_mass,
8.34336671824457987e+00,
4.12479856412430479e+00,-4.03523417114321381e-01,
-2.76742510726862411e-03, 4.99852801234917238e-03,
2.30417297573763929e-05);
Planet uranus( 4.36624404335156298e-05 * solar_mass,

1.28943695621391310e+01,-1.51111514016986312e+01,-2.23307578892655734e-01,
2.96460137564761618e-03,
2.37847173959480950e-03,-2.96589568540237556e-05);
Planet neptune( 5.15138902046611451e-05 * solar_mass,
1.53796971148509165e+01,-2.59193146099879641e+01,
1.79258772950371181e-01,
2.68067772490389322e-03,
1.62824170038242295e-03,-9.51592254519715870e-05);

Planet tmp[] = {sun, jupiter, saturn, uranus, neptune};
vector<Planet> bodies(tmp,tmp+NBODIES);
transform(bodies.begin(),bodies.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;
}

 

[vishal.k.s]

Good job. I am not very surprised. The last time I looked at the
source code for some of the C++ programs on that site, I was
able to write shorter and faster C++ equivalents for at least
6 of the listed programs.

I entirely agree with you, the C++ programs need to be fixed,
firstly for performance and secondly for number of lines of code...
as you put it TLC.

 

[Greg Buchholz]

I entirely agree with you, the C++ programs need to be fixed,
firstly for performance and secondly for number of lines of code...
as you put it TLC.

Here's another update. I got almost all of the speed back (~10%
slower) by manually unrolling the implicit inner valarray loops in
advance(), while keeping the LOC constant. Maybe a sufficiently smart
compiler could figure this out on its own. The valarrays are a
constant length of 3 after all. Oh, and it does use a somewhat ugly
looking macro...

http://sleepingsquirrel.org/cpp/body6.cpp.html

....and here's a slow (3x) but shorter version (64 non-trivial lines of
code)...

http://sleepingsquirrel.org/cpp/body7.cpp.html

Greg Buchholz