| OLD | NEW |
| (Empty) |
| 1 /* The Computer Language Benchmarks Game | |
| 2 * http://benchmarksgame.alioth.debian.org/ | |
| 3 * | |
| 4 * contributed by Christoph Bauer | |
| 5 * slightly sped up by Petr Prokhorenkov | |
| 6 */ | |
| 7 | |
| 8 #include <math.h> | |
| 9 #include <stdio.h> | |
| 10 #include <stdlib.h> | |
| 11 | |
| 12 #define pi 3.141592653589793 | |
| 13 #define solar_mass (4 * pi * pi) | |
| 14 #define days_per_year 365.24 | |
| 15 | |
| 16 struct planet { | |
| 17 double x, y, z; | |
| 18 double vx, vy, vz; | |
| 19 double mass; | |
| 20 }; | |
| 21 | |
| 22 /* | |
| 23 * Here's one weird thing: inlining of this function | |
| 24 * decreases performance by 25%. (I.e. do not compile this with -O3) | |
| 25 * Advances with dt == 1.0 | |
| 26 */ | |
| 27 void advance(int nbodies, struct planet * bodies) | |
| 28 { | |
| 29 int i, j; | |
| 30 | |
| 31 for (i = 0; i < nbodies; i++) { | |
| 32 struct planet * b = &(bodies[i]); | |
| 33 for (j = i + 1; j < nbodies; j++) { | |
| 34 struct planet * b2 = &(bodies[j]); | |
| 35 double dx = b->x - b2->x; | |
| 36 double dy = b->y - b2->y; | |
| 37 double dz = b->z - b2->z; | |
| 38 double inv_distance = 1.0/sqrt(dx * dx + dy * dy + dz * dz); | |
| 39 double mag = inv_distance * inv_distance * inv_distance; | |
| 40 b->vx -= dx * b2->mass * mag; | |
| 41 b->vy -= dy * b2->mass * mag; | |
| 42 b->vz -= dz * b2->mass * mag; | |
| 43 b2->vx += dx * b->mass * mag; | |
| 44 b2->vy += dy * b->mass * mag; | |
| 45 b2->vz += dz * b->mass * mag; | |
| 46 } | |
| 47 } | |
| 48 for (i = 0; i < nbodies; i++) { | |
| 49 struct planet * b = &(bodies[i]); | |
| 50 b->x += b->vx; | |
| 51 b->y += b->vy; | |
| 52 b->z += b->vz; | |
| 53 } | |
| 54 } | |
| 55 | |
| 56 double energy(int nbodies, struct planet * bodies) | |
| 57 { | |
| 58 double e; | |
| 59 int i, j; | |
| 60 | |
| 61 e = 0.0; | |
| 62 for (i = 0; i < nbodies; i++) { | |
| 63 struct planet * b = &(bodies[i]); | |
| 64 e += 0.5 * b->mass * (b->vx * b->vx + b->vy * b->vy + b->vz * b->vz); | |
| 65 for (j = i + 1; j < nbodies; j++) { | |
| 66 struct planet * b2 = &(bodies[j]); | |
| 67 double dx = b->x - b2->x; | |
| 68 double dy = b->y - b2->y; | |
| 69 double dz = b->z - b2->z; | |
| 70 double distance = sqrt(dx * dx + dy * dy + dz * dz); | |
| 71 e -= (b->mass * b2->mass) / distance; | |
| 72 } | |
| 73 } | |
| 74 return e; | |
| 75 } | |
| 76 | |
| 77 void offset_momentum(int nbodies, struct planet * bodies) | |
| 78 { | |
| 79 double px = 0.0, py = 0.0, pz = 0.0; | |
| 80 int i; | |
| 81 for (i = 0; i < nbodies; i++) { | |
| 82 px += bodies[i].vx * bodies[i].mass; | |
| 83 py += bodies[i].vy * bodies[i].mass; | |
| 84 pz += bodies[i].vz * bodies[i].mass; | |
| 85 } | |
| 86 bodies[0].vx = - px / solar_mass; | |
| 87 bodies[0].vy = - py / solar_mass; | |
| 88 bodies[0].vz = - pz / solar_mass; | |
| 89 } | |
| 90 | |
| 91 #define NBODIES 5 | |
| 92 struct planet bodies[NBODIES] = { | |
| 93 { /* sun */ | |
| 94 0, 0, 0, 0, 0, 0, solar_mass | |
| 95 }, | |
| 96 { /* jupiter */ | |
| 97 4.84143144246472090e+00, | |
| 98 -1.16032004402742839e+00, | |
| 99 -1.03622044471123109e-01, | |
| 100 1.66007664274403694e-03 * days_per_year, | |
| 101 7.69901118419740425e-03 * days_per_year, | |
| 102 -6.90460016972063023e-05 * days_per_year, | |
| 103 9.54791938424326609e-04 * solar_mass | |
| 104 }, | |
| 105 { /* saturn */ | |
| 106 8.34336671824457987e+00, | |
| 107 4.12479856412430479e+00, | |
| 108 -4.03523417114321381e-01, | |
| 109 -2.76742510726862411e-03 * days_per_year, | |
| 110 4.99852801234917238e-03 * days_per_year, | |
| 111 2.30417297573763929e-05 * days_per_year, | |
| 112 2.85885980666130812e-04 * solar_mass | |
| 113 }, | |
| 114 { /* uranus */ | |
| 115 1.28943695621391310e+01, | |
| 116 -1.51111514016986312e+01, | |
| 117 -2.23307578892655734e-01, | |
| 118 2.96460137564761618e-03 * days_per_year, | |
| 119 2.37847173959480950e-03 * days_per_year, | |
| 120 -2.96589568540237556e-05 * days_per_year, | |
| 121 4.36624404335156298e-05 * solar_mass | |
| 122 }, | |
| 123 { /* neptune */ | |
| 124 1.53796971148509165e+01, | |
| 125 -2.59193146099879641e+01, | |
| 126 1.79258772950371181e-01, | |
| 127 2.68067772490389322e-03 * days_per_year, | |
| 128 1.62824170038242295e-03 * days_per_year, | |
| 129 -9.51592254519715870e-05 * days_per_year, | |
| 130 5.15138902046611451e-05 * solar_mass | |
| 131 } | |
| 132 }; | |
| 133 | |
| 134 #define DT 1e-2 | |
| 135 #define RECIP_DT (1.0/DT) | |
| 136 | |
| 137 /* | |
| 138 * Rescale certain properties of bodies. That allows doing | |
| 139 * consequential advance()'s as if dt were equal to 1.0. | |
| 140 * | |
| 141 * When all advances done, rescale bodies back to obtain correct energy. | |
| 142 */ | |
| 143 void scale_bodies(int nbodies, struct planet * bodies, double scale) { | |
| 144 int i; | |
| 145 | |
| 146 for (i = 0; i < nbodies; i++) { | |
| 147 bodies[i].mass *= scale*scale; | |
| 148 bodies[i].vx *= scale; | |
| 149 bodies[i].vy *= scale; | |
| 150 bodies[i].vz *= scale; | |
| 151 } | |
| 152 } | |
| 153 | |
| 154 int main(int argc, char ** argv) | |
| 155 { | |
| 156 int n = atoi(argv[1]); | |
| 157 int i; | |
| 158 | |
| 159 offset_momentum(NBODIES, bodies); | |
| 160 printf ("%.9f\n", energy(NBODIES, bodies)); | |
| 161 scale_bodies(NBODIES, bodies, DT); | |
| 162 for (i = 1; i <= n; i++) { | |
| 163 advance(NBODIES, bodies); | |
| 164 } | |
| 165 scale_bodies(NBODIES, bodies, RECIP_DT); | |
| 166 printf ("%.9f\n", energy(NBODIES, bodies)); | |
| 167 return 0; | |
| 168 } | |
| 169 | |
| OLD | NEW |