Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- # The Computer Language Benchmarks Game
- # http://benchmarksgame.alioth.debian.org/
- #
- # originally by Kevin Carson
- # modified by Tupteq, Fredrik Johansson, and Daniel Nanz
- # modified by Maciej Fijalkowski
- # 2to3
- from libc.math cimport pow, sqrt
- import sys
- PI = 3.14159265358979323
- SOLAR_MASS = 4 * PI * PI
- DAYS_PER_YEAR = 365.24
- cdef struct body_t:
- double x[3]
- double v[3]
- double m
- DEF NBODIES = 5
- BODIES = {
- 'sun': ([0.0, 0.0, 0.0], [0.0, 0.0, 0.0], 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)
- }
- def combinations(l):
- result = []
- for x in range(len(l) - 1):
- ls = l[x+1:]
- for y in ls:
- result.append((l[x],y))
- return result
- cdef void make_cbodies(list bodies, body_t *cbodies, int num_cbodies):
- cdef body_t *cbody
- for i, body in enumerate(bodies):
- if i >= num_cbodies:
- break
- (x, v, m) = body
- cbody = &cbodies[i]
- cbody.x[0], cbody.x[1], cbody.x[2] = x
- cbody.v[0], cbody.v[1], cbody.v[2] = v
- cbodies[i].m = m
- cdef list make_pybodies(body_t *cbodies, int num_cbodies):
- pybodies = []
- for i in range(num_cbodies):
- x = [cbodies[i].x[0], cbodies[i].x[1], cbodies[i].x[2]]
- v = [cbodies[i].v[0], cbodies[i].v[1], cbodies[i].v[2]]
- pybodies.append((x, v, cbodies[i].m))
- return pybodies
- def advance(double dt, int n, bodies):
- # Note that this does not take a `pairs` argument, and *returns* a new
- # `bodies` list. This is slightly different from the original pure-Python
- # version of `advance`.
- cdef:
- int i, ii, jj
- double dx, dy, dz, mag, b1m, b2m, ds
- body_t *body1
- body_t *body2
- body_t cbodies[NBODIES]
- make_cbodies(bodies, cbodies, NBODIES)
- for i in range(n):
- for ii in range(NBODIES-1):
- body1 = &cbodies[ii]
- for jj in range(ii+1, NBODIES):
- body2 = &cbodies[jj]
- dx = body1.x[0] - body2.x[0]
- dy = body1.x[1] - body2.x[1]
- dz = body1.x[2] - body2.x[2]
- ds = dx * dx + dy * dy + dz * dz
- mag = dt / (ds * sqrt(ds))
- b1m = body1.m * mag
- b2m = body2.m * mag
- body1.v[0] -= dx * b2m
- body1.v[1] -= dy * b2m
- body1.v[2] -= dz * b2m
- body2.v[0] += dx * b1m
- body2.v[1] += dy * b1m
- body2.v[2] += dz * b1m
- for ii in range(NBODIES):
- body2 = &cbodies[ii]
- body2.x[0] += dt * body2.v[0]
- body2.x[1] += dt * body2.v[1]
- body2.x[2] += dt * body2.v[2]
- return make_pybodies(cbodies, NBODIES)
- def report_energy(bodies):
- # Note that this takes just a `bodies` list, and computes `pairs`
- # internally, which is a slight modification from the original pure-Python
- # version of `report_energy`.
- e = 0.0
- pairs = combinations(bodies)
- for (((x1, y1, z1), v1, m1),
- ((x2, y2, z2), v2, m2)) in pairs:
- dx = x1 - x2
- dy = y1 - y2
- dz = z1 - z2
- e -= (m1 * m2) / ((dx * dx + dy * dy + dz * dz) ** 0.5)
- for (r, [vx, vy, vz], m) in bodies:
- e += m * (vx * vx + vy * vy + vz * vz) / 2.
- print("%.9f" % e)
- def offset_momentum(ref, bodies):
- px = py = pz = 0.0
- for (r, [vx, vy, vz], m) in bodies:
- px -= vx * m
- py -= vy * m
- pz -= vz * m
- (r, v, m) = ref
- v[0] = px / m
- v[1] = py / m
- v[2] = pz / m
- def main(n, bodies=BODIES, ref='sun'):
- system = list(bodies.values())
- offset_momentum(bodies[ref], system)
- report_energy(system)
- system = advance(0.01, n, system)
- report_energy(system)
- [chris@Vonbrigdum 02-cython]$ cd ../03-pure-c/
- [chris@Vonbrigdum 03-pure-c]$ ls
- Makefile nbody.c
- [chris@Vonbrigdum 03-pure-c]$ cat nbody.c
- /* The Computer Language Benchmarks Game
- * http://benchmarksgame.alioth.debian.org/
- *
- * contributed by Christoph Bauer
- * slightly sped up by Petr Prokhorenkov
- */
- #include <math.h>
- #include <stdio.h>
- #include <stdlib.h>
- #define pi 3.141592653589793
- #define solar_mass (4 * pi * pi)
- #define days_per_year 365.24
- struct planet {
- double x, y, z;
- double vx, vy, vz;
- double mass;
- };
- /*
- * Here's one weird thing: inlining of this function
- * decreases performance by 25%. (I.e. do not compile this with -O3)
- * Advances with dt == 1.0
- */
- void advance(int nbodies, struct planet * bodies)
- {
- int i, j;
- for (i = 0; i < nbodies; i++) {
- struct planet * b = &(bodies[i]);
- for (j = i + 1; j < nbodies; j++) {
- struct planet * b2 = &(bodies[j]);
- double dx = b->x - b2->x;
- double dy = b->y - b2->y;
- double dz = b->z - b2->z;
- double inv_distance = 1.0/sqrt(dx * dx + dy * dy + dz * dz);
- double mag = inv_distance * inv_distance * inv_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;
- }
- }
- for (i = 0; i < nbodies; i++) {
- struct planet * b = &(bodies[i]);
- b->x += b->vx;
- b->y += b->vy;
- b->z += b->vz;
- }
- }
- double energy(int nbodies, struct planet * bodies)
- {
- double e;
- int i, j;
- e = 0.0;
- for (i = 0; i < nbodies; i++) {
- struct planet * b = &(bodies[i]);
- e += 0.5 * b->mass * (b->vx * b->vx + b->vy * b->vy + b->vz * b->vz);
- for (j = i + 1; j < nbodies; j++) {
- struct planet * b2 = &(bodies[j]);
- 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;
- }
- }
- return e;
- }
- void offset_momentum(int nbodies, struct planet * bodies)
- {
- double px = 0.0, py = 0.0, pz = 0.0;
- int i;
- for (i = 0; i < nbodies; i++) {
- px += bodies[i].vx * bodies[i].mass;
- py += bodies[i].vy * bodies[i].mass;
- pz += bodies[i].vz * bodies[i].mass;
- }
- bodies[0].vx = - px / solar_mass;
- bodies[0].vy = - py / solar_mass;
- bodies[0].vz = - pz / solar_mass;
- }
- #define NBODIES 5
- struct planet bodies[NBODIES] = {
- { /* sun */
- 0, 0, 0, 0, 0, 0, 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
- }
- };
- #define DT 1e-2
- #define RECIP_DT (1.0/DT)
- /*
- * Rescale certain properties of bodies. That allows doing
- * consequential advance()'s as if dt were equal to 1.0.
- *
- * When all advances done, rescale bodies back to obtain correct energy.
- */
- void scale_bodies(int nbodies, struct planet * bodies, double scale) {
- int i;
- for (i = 0; i < nbodies; i++) {
- bodies[i].mass *= scale*scale;
- bodies[i].vx *= scale;
- bodies[i].vy *= scale;
- bodies[i].vz *= scale;
- }
- }
- int main(int argc, char ** argv)
- {
- int n = atoi(argv[1]);
- int i;
- offset_momentum(NBODIES, bodies);
- printf ("%.9f\n", energy(NBODIES, bodies));
- scale_bodies(NBODIES, bodies, DT);
- for (i = 1; i <= n; i++) {
- advance(NBODIES, bodies);
- }
- scale_bodies(NBODIES, bodies, RECIP_DT);
- printf ("%.9f\n", energy(NBODIES, bodies));
- return 0;
- }
Add Comment
Please, Sign In to add comment