examples/shootout/nbody.pir

   1 #!./parrot
   2 # Copyright (C) 2006-2009, Parrot Foundation.
   3 # $Id$
   4 # The Computer Language Shootout
   5 # http://shootout.alioth.debian.org/
   6 #
   7 # ./parrot -R jit
   8 # Contributed by Joshua Isom
   9 # speed up  from 1m25 to 6s by Leopold Toetsch
  10 # changed default value to 1000 to match shootout default (karl)
  11
  12 .const int x = 0
  13 .const int y = 1
  14 .const int z = 2
  15 .const int vx = 3
  16 .const int vy = 4
  17 .const int vz = 5
  18 .const int m = 6
  19
  20 # save on repetition of code
  21 .macro InitBodies (bodies, i, num1, num2, num3, num4, num5, num6, num7)
  22         $N0 = .num4 * 365.24
  23         $N1 = .num5 * 365.24
  24         $N2 = .num6 * 365.24
  25         $N3 = .num7 * 39.478417604357428
  26         $P0 = new 'FixedFloatArray'
  27         $P0 = 7
  28         .bodies[.i] = $P0
  29         .bodies[.i; x] = .num1
  30         .bodies[.i; y] = .num2
  31         .bodies[.i; z] = .num3
  32         .bodies[.i; vx] = $N0
  33         .bodies[.i; vy] = $N1
  34         .bodies[.i; vz] = $N2
  35         .bodies[.i; m] = $N3
  36 .endm
  37
  38 .sub main :main
  39         .param pmc argv
  40         .local int argc, n, nbodies
  41         argc = argv
  42         n = 1000
  43         unless argc > 1 goto argsok
  44         $S0 = argv[1]
  45         n = $S0
  46 argsok:
  47         .local pmc bodies
  48         bodies = new 'FixedPMCArray'
  49         bodies = 5
  50         # Sun
  51         .InitBodies(bodies, 0, 0, 0, 0, 0, 0, 0, 1)
  52         # Jupiter
  53         .InitBodies(bodies, 1,
  54          4.84143144246472090e+00,
  55         -1.16032004402742839e+00,
  56         -1.03622044471123109e-01,
  57          1.66007664274403694e-03,
  58          7.69901118419740425e-03,
  59         -6.90460016972063023e-05,
  60          9.54791938424326609e-04)
  61         # Saturn
  62         .InitBodies(bodies, 2,
  63          8.34336671824457987e+00,
  64          4.12479856412430479e+00,
  65         -4.03523417114321381e-01,
  66         -2.76742510726862411e-03,
  67          4.99852801234917238e-03,
  68          2.30417297573763929e-05,
  69          2.85885980666130812e-04)
  70         # Uranus
  71         .InitBodies(bodies, 3,
  72          1.28943695621391310e+01,
  73         -1.51111514016986312e+01,
  74         -2.23307578892655734e-01,
  75          2.96460137564761618e-03,
  76          2.37847173959480950e-03,
  77         -2.96589568540237556e-05,
  78          4.36624404335156298e-05)
  79         # Neptune
  80         .InitBodies(bodies, 4,
  81          1.53796971148509165e+01,
  82         -2.59193146099879641e+01,
  83          1.79258772950371181e-01,
  84          2.68067772490389322e-03,
  85          1.62824170038242295e-03,
  86         -9.51592254519715870e-05,
  87          5.15138902046611451e-05)
  88
  89         nbodies = bodies
  90         offset_momentum(nbodies, bodies)
  91         $N0 = energy(nbodies, bodies)
  92         .local pmc spf
  93         spf = new 'FixedFloatArray'
  94         spf = 1
  95         spf[0] = $N0
  96         $S0 = sprintf "%.9f\n", spf
  97         print $S0
  98         .local int i
  99         i = 0
 100 beginfor:
 101         unless i < n goto endfor
 102                 advance(nbodies, bodies, 0.01)
 103         inc i
 104         goto beginfor
 105 endfor:
 106
 107         $N0 = energy(nbodies, bodies)
 108         spf[0] = $N0
 109         $S0 = sprintf "%.9f\n", spf
 110         print $S0
 111         exit 0
 112 .end
 113
 114 .sub offset_momentum
 115         .param int nbodies
 116         .param pmc bodies
 117         .local num px, py, pz
 118         px = 0.0
 119         py = 0.0
 120         pz = 0.0
 121         .local int i
 122         i = 0
 123 beginfor:
 124         unless i < nbodies goto endfor
 125         $N0 = bodies[i; vx]
 126         $N1 = bodies[i; m]
 127         $N0 *= $N1
 128         px += $N0
 129         $N0 = bodies[i; vy]
 130         $N1 = bodies[i; m]
 131         $N0 *= $N1
 132         py += $N0
 133         $N0 = bodies[i; vz]
 134         $N1 = bodies[i; m]
 135         $N0 *= $N1
 136         pz += $N0
 137         inc i
 138         goto beginfor
 139 endfor:
 140         # bodies[0].vx = - px / solar_mass;
 141         px /= -39.478417604357428
 142         bodies[0; vx] = px
 143         py /= -39.478417604357428
 144         bodies[0; vy] = py
 145         pz /= -39.478417604357428
 146         bodies[0; vz] = pz
 147 .end
 148
 149 .sub energy
 150         .param int nbodies
 151         .param pmc bodies
 152         .local num e, tmp
 153         .local int i, j
 154         e = 0.0
 155         i = 0
 156 beginfor_0:
 157         unless i < nbodies goto endfor_0
 158
 159         # e += 0.5 * b->mass * (b->vx * b->vx + b->vy * b->vy + b->vz * b->vz);
 160         $N0 = bodies[i; m] # mass
 161         $N0 *= 0.5
 162
 163         $N1 = bodies[i; vx] # vx
 164         $N3 = pow $N1, 2.0
 165
 166         $N1 = bodies[i; vy] # vy
 167         $N2 = pow $N1, 2.0
 168         $N3 += $N2
 169
 170         $N1 = bodies[i; vz] # vz
 171         $N2 = pow $N1, 2.0
 172         $N3 += $N2
 173
 174         $N0 *= $N3
 175
 176         e += $N0
 177
 178         j = i + 1
 179 beginfor_1:
 180         unless j < nbodies goto endfor_1
 181         .local num dx, dy, dz, distance
 182
 183         # dx = b->x - b2->x;
 184         $N0 = bodies[i; x]
 185         $N1 = bodies[j; x]
 186         dx = $N0 - $N1
 187
 188         # dy = b->y - b2->y;
 189         $N0 = bodies[i; y]
 190         $N1 = bodies[j; y]
 191         dy = $N0 - $N1
 192
 193         # dz = b->z - b2->z;
 194         $N0 = bodies[i; z]
 195         $N1 = bodies[j; z]
 196         dz = $N0 - $N1
 197
 198         # distance = sqrt(dx * dx + dy * dy + dz * dz);
 199         $N0 = dx * dx
 200         $N1 = dy * dy
 201         $N2 = dz * dz
 202         $N0 += $N1
 203         $N0 += $N2
 204         distance = sqrt $N0
 205
 206         # e -= (b->mass * b2->mass) / distance;
 207         $N0 = bodies[i; m]
 208         $N1 = bodies[j; m]
 209         $N0 *= $N1
 210         $N0 /= distance
 211         e -= $N0
 212
 213         inc j
 214         goto beginfor_1
 215 endfor_1:
 216
 217         inc i
 218         goto beginfor_0
 219 endfor_0:
 220         .return(e)
 221 .end
 222
 223 .sub advance
 224         .param int nbodies
 225         .param pmc bodies
 226         .param num dt
 227         .local int i, j
 228         .local num dx, dy, dz, distance, mag
 229         .local num bx, by, bz, bm, bvx, bvy, bvz
 230         .local num b2x, b2y, b2z, b2m
 231         .local pmc bi, bj
 232         i = 0
 233 beginfor_0:
 234         unless i < nbodies goto endfor_0
 235         bi = bodies[i]
 236         bx = bi[x]
 237         by = bi[y]
 238         bz = bi[z]
 239         bm = bi[m]
 240         bvx = bi[vx]
 241         bvy = bi[vy]
 242         bvz = bi[vz]
 243
 244         j = i + 1
 245         beginfor_1:
 246                 unless j < nbodies goto endfor_1
 247                 bj = bodies[j]
 248                 b2x = bj[x]
 249                 b2y = bj[y]
 250                 b2z = bj[z]
 251                 b2m = bj[m]
 252
 253                 # dx = b->x - b2->x;
 254                 dx = bx - b2x
 255                 # dy = b->y - b2->y;
 256                 dy = by - b2y
 257                 # dz = b->z - b2->z;
 258                 dz = bz - b2z
 259
 260                 # distance = sqrt(dx * dx + dy * dy + dz * dz);
 261                 $N0 = dx * dx
 262                 $N1 = dy * dy
 263                 $N0 += $N1
 264                 $N1 = dz * dz
 265                 $N0 += $N1
 266                 distance = sqrt $N0
 267
 268                 # mag = dt / (distance * distance * distance);
 269                 $N0 = distance * distance
 270                 $N0 *= distance
 271                 mag = dt / $N0
 272
 273                 # b->vx -= dx * b2->mass * mag;
 274                 $N0 = dx * b2m
 275                 $N0 *= mag
 276                 bvx -= $N0
 277
 278                 # b->vy -= dy * b2->mass * mag;
 279                 $N0 = dy * b2m
 280                 $N0 *= mag
 281                 bvy -= $N0
 282
 283                 # b->vz -= dz * b2->mass * mag;
 284                 $N0 = dz * b2m
 285                 $N0 *= mag
 286                 bvz -= $N0
 287
 288                 # b2->vx += dx * b->mass * mag;
 289                 $N0 = dx * bm
 290                 $N0 *= mag
 291                 $N1 = bj[vx]
 292                 $N1 += $N0
 293                 bj[vx] = $N1
 294
 295                 # b2->vy += dy * b->mass * mag;
 296                 $N0 = dy * bm
 297                 $N0 *= mag
 298                 $N1 = bj[vy]
 299                 $N1 += $N0
 300                 bj[vy] = $N1
 301
 302                 # b2->vz += dz * b->mass * mag;
 303                 $N0 = dz * bm
 304                 $N0 *= mag
 305                 $N1 = bj[vz]
 306                 $N1 += $N0
 307                 bj[vz] = $N1
 308
 309                 inc j
 310                 goto beginfor_1
 311         endfor_1:
 312         bi[vx] = bvx
 313         bi[vy] = bvy
 314         bi[vz] = bvz
 315
 316         inc i
 317         goto beginfor_0
 318 endfor_0:
 319
 320         i = 0
 321 beginfor_2:
 322         unless i < nbodies goto endfor_2
 323         bi = bodies[i]
 324         # b->x += dt * b->vx;
 325         $N0 = bi[vx]
 326         $N1 = dt * $N0
 327         $N0 = bi[x]
 328         $N0 += $N1
 329         bi[x] = $N0
 330
 331         # b->y += dt * b->vy;
 332         $N0 = bi[vy]
 333         $N1 = dt * $N0
 334         $N0 = bi[y]
 335         $N0 += $N1
 336         bi[y] = $N0
 337
 338         # b->z += dt * b->vz;
 339         $N0 = bi[vz]
 340         $N1 = dt * $N0
 341         $N0 = bi[z]
 342         $N0 += $N1
 343         bi[z] = $N0
 344
 345         inc i
 346         goto beginfor_2
 347 endfor_2:
 348
 349 .end
 350
 351 # Local Variables:
 352 #   mode: pir
 353 #   fill-column: 100
 354 # End:
 355 # vim: expandtab shiftwidth=4 ft=pir: