|
| 1 | +// based on: |
| 2 | +// http://shootout.alioth.debian.org/u32/benchmark.php?test=nbody&lang=java |
| 3 | + |
| 4 | +fn main() { |
| 5 | + |
| 6 | + let vec[int] inputs = vec( |
| 7 | + 50000 |
| 8 | + //these segfault :( |
| 9 | + //500000, |
| 10 | + //5000000, |
| 11 | + //50000000 |
| 12 | + ); |
| 13 | + |
| 14 | + let vec[Body.props] bodies = NBodySystem.MakeNBodySystem(); |
| 15 | + |
| 16 | + for (int n in inputs) { |
| 17 | + // TODO: make #fmt handle floats? |
| 18 | + log NBodySystem.energy(bodies); |
| 19 | + |
| 20 | + let int i = 0; |
| 21 | + while (i < n) { |
| 22 | + bodies = NBodySystem.advance(bodies, 0.01); |
| 23 | + i = i+1; |
| 24 | + } |
| 25 | + log NBodySystem.energy(bodies); |
| 26 | + } |
| 27 | +} |
| 28 | + |
| 29 | +// making a native call to sqrt |
| 30 | +native "rust" mod rustrt { |
| 31 | + fn squareroot(&float input, &mutable float output); |
| 32 | +} |
| 33 | + |
| 34 | +// Body.props is a record of floats, so |
| 35 | +// vec[Body.props] is a vector of records of floats |
| 36 | + |
| 37 | +mod NBodySystem { |
| 38 | + |
| 39 | + fn MakeNBodySystem() -> vec[Body.props] { |
| 40 | + // can't iterate over a record? how about a vector, then? |
| 41 | + let vec[Body.props] bodies = vec( |
| 42 | + // these each return a Body.props |
| 43 | + Body.sun(), |
| 44 | + Body.jupiter(), |
| 45 | + Body.saturn(), |
| 46 | + Body.uranus(), |
| 47 | + Body.neptune()); |
| 48 | + |
| 49 | + let float px = 0.0; |
| 50 | + let float py = 0.0; |
| 51 | + let float pz = 0.0; |
| 52 | + |
| 53 | + for (Body.props body in bodies) { |
| 54 | + px += body.vx * body.mass; |
| 55 | + py += body.vy * body.mass; |
| 56 | + pz += body.vz * body.mass; |
| 57 | + } |
| 58 | + bodies.(0) = Body.offsetMomentum(bodies.(0), px, py, pz); |
| 59 | + |
| 60 | + ret bodies; |
| 61 | + } |
| 62 | + |
| 63 | + fn advance(vec[Body.props] bodies, float dt) -> vec[Body.props] { |
| 64 | + for (Body.props ibody in bodies) { |
| 65 | + |
| 66 | + let Body.props iBody = ibody; |
| 67 | + |
| 68 | + for (Body.props jbody in bodies) { |
| 69 | + let float dx = iBody.x - jbody.x; |
| 70 | + let float dy = iBody.y - jbody.y; |
| 71 | + let float dz = iBody.z - jbody.z; |
| 72 | + |
| 73 | + let float dSquared = dx * dx + dy * dy + dz * dz; |
| 74 | + |
| 75 | + let float distance; |
| 76 | + rustrt.squareroot(dSquared, distance); |
| 77 | + let float mag = dt / (dSquared * distance); |
| 78 | + } |
| 79 | + } |
| 80 | + |
| 81 | + for (Body.props body in bodies) { |
| 82 | + body.x += dt * body.vx; |
| 83 | + body.y += dt * body.vy; |
| 84 | + body.z += dt * body.vz; |
| 85 | + } |
| 86 | + |
| 87 | + ret bodies; |
| 88 | + } |
| 89 | + |
| 90 | + fn energy(vec[Body.props] bodies) -> float { |
| 91 | + let float dx; |
| 92 | + let float dy; |
| 93 | + let float dz; |
| 94 | + let float distance; |
| 95 | + let float e = 0.0; |
| 96 | + |
| 97 | + for (Body.props ibody in bodies) { |
| 98 | + |
| 99 | + // do we need this? |
| 100 | + let Body.props iBody = ibody; |
| 101 | + |
| 102 | + e += 0.5 * iBody.mass * |
| 103 | + ( iBody.vx * iBody.vx |
| 104 | + + iBody.vy * iBody.vy |
| 105 | + + iBody.vz * iBody.vz ); |
| 106 | + |
| 107 | + for (Body.props jbody in bodies) { |
| 108 | + |
| 109 | + // do we need this? |
| 110 | + let Body.props jBody = jbody; |
| 111 | + |
| 112 | + dx = iBody.x - jBody.x; |
| 113 | + dy = iBody.y - jBody.y; |
| 114 | + dz = iBody.z - jBody.z; |
| 115 | + |
| 116 | + rustrt.squareroot(dx*dx + dy*dy + dz*dz, distance); |
| 117 | + e -= (iBody.mass * jBody.mass) / distance; |
| 118 | + } |
| 119 | + } |
| 120 | + ret e; |
| 121 | + } |
| 122 | + |
| 123 | +} |
| 124 | + |
| 125 | +mod Body { |
| 126 | + |
| 127 | + const float PI = 3.141592; |
| 128 | + const float SOLAR_MASS = 39.478417; // was 4 * PI * PI originally |
| 129 | + const float DAYS_PER_YEAR = 365.24; |
| 130 | + |
| 131 | + type props = rec(float x, |
| 132 | + float y, |
| 133 | + float z, |
| 134 | + float vx, |
| 135 | + float vy, |
| 136 | + float vz, |
| 137 | + float mass); |
| 138 | + |
| 139 | + fn jupiter() -> Body.props { |
| 140 | + // current limitation of the float lexer: decimal part has to |
| 141 | + // fit into a 32-bit int. |
| 142 | + |
| 143 | + let Body.props p; |
| 144 | + p.x = 4.841431e+00; |
| 145 | + p.y = -1.160320e+00; |
| 146 | + p.z = -1.036220e-01; |
| 147 | + p.vx = 1.660076e-03 * DAYS_PER_YEAR; |
| 148 | + p.vy = 7.699011e-03 * DAYS_PER_YEAR; |
| 149 | + p.vz = -6.904600e-05 * DAYS_PER_YEAR; |
| 150 | + p.mass = 9.547919e-04 * SOLAR_MASS; |
| 151 | + ret p; |
| 152 | + } |
| 153 | + |
| 154 | + fn saturn() -> Body.props { |
| 155 | + let Body.props p; |
| 156 | + p.x = 8.343366e+00; |
| 157 | + p.y = 4.124798e+00; |
| 158 | + p.z = -4.035234e-01; |
| 159 | + p.vx = -2.767425e-03 * DAYS_PER_YEAR; |
| 160 | + p.vy = 4.998528e-03 * DAYS_PER_YEAR; |
| 161 | + p.vz = 2.304172e-05 * DAYS_PER_YEAR; |
| 162 | + p.mass = 2.858859e-04 * SOLAR_MASS; |
| 163 | + ret p; |
| 164 | + } |
| 165 | + |
| 166 | + fn uranus() -> Body.props { |
| 167 | + let Body.props p; |
| 168 | + p.x = 1.289436e+01; |
| 169 | + p.y = -1.511115e+01; |
| 170 | + p.z = -2.233075e-01; |
| 171 | + p.vx = 2.964601e-03 * DAYS_PER_YEAR; |
| 172 | + p.vy = 2.378471e-03 * DAYS_PER_YEAR; |
| 173 | + p.vz = -2.965895e-05 * DAYS_PER_YEAR; |
| 174 | + p.mass = 4.366244e-05 * SOLAR_MASS; |
| 175 | + ret p; |
| 176 | + } |
| 177 | + |
| 178 | + fn neptune() -> Body.props { |
| 179 | + let Body.props p; |
| 180 | + p.x = 1.537969e+01; |
| 181 | + p.y = -2.591931e+01; |
| 182 | + p.z = 1.792587e-01; |
| 183 | + p.vx = 2.680677e-03 * DAYS_PER_YEAR; |
| 184 | + p.vy = 1.628241e-03 * DAYS_PER_YEAR; |
| 185 | + p.vz = -9.515922e-05 * DAYS_PER_YEAR; |
| 186 | + p.mass = 5.151389e-05 * SOLAR_MASS; |
| 187 | + ret p; |
| 188 | + } |
| 189 | + |
| 190 | + fn sun() -> Body.props { |
| 191 | + let Body.props p; |
| 192 | + p.mass = SOLAR_MASS; |
| 193 | + ret p; |
| 194 | + } |
| 195 | + |
| 196 | + fn offsetMomentum(Body.props props, |
| 197 | + float px, |
| 198 | + float py, |
| 199 | + float pz) -> Body.props { |
| 200 | + |
| 201 | + // TODO: should we create a new one or mutate the original? |
| 202 | + let Body.props p = props; |
| 203 | + p.vx = -px / SOLAR_MASS; |
| 204 | + p.vy = -py / SOLAR_MASS; |
| 205 | + p.vz = -pz / SOLAR_MASS; |
| 206 | + ret p; |
| 207 | + } |
| 208 | + |
| 209 | +} |
0 commit comments