NnBody
N-Body
See http://openeuphoria.org/forum/119406.wc
N-Body Euphoria
-- The Computer Language Shootout Benchmarks -- http://shootout.alioth.debian.org/ -- -- Converted to Euphoria by Jason Gade -- Optimized by Matt Lewis -- run: eui nbody.ex N without warning without type_check include get.e constant PI = 3.141592653589793, SOLAR_MASS = 4 * PI * PI, DAYS_PER_YEAR = 365.24 -- struct planet constant name = 1, x = 2, y = 3, z = 4, vx = 5, vy = 6, vz = 7, mass = 8 -- end struct sequence bodies bodies = { { "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 } } constant NBODIES = length(bodies) sequence NAME = {}, X = {}, Y = {}, Z = {}, VX = {}, VY = {}, VZ = {}, MASS = {} for i = 1 to NBODIES do for f = name to mass do switch f do case name then NAME = append( NAME, bodies[i][f] ) case x then X &= bodies[i][f] case y then Y &= bodies[i][f] case z then Z &= bodies[i][f] case vx then VX &= bodies[i][f] case vy then VY &= bodies[i][f] case vz then VZ &= bodies[i][f] case mass then MASS &= bodies[i][f] end switch end for end for procedure advance(atom dt) atom dx, dy, dz, distance, mag, mass, distance_2, mass_mag for i = 1 to NBODIES do for j = i + 1 to NBODIES do dx = X[i] - X[j] dy = Y[i] - Y[j] dz = Z[i] - Z[j] distance_2 = dx*dx distance_2 += dy*dy distance_2 += dz*dz distance = sqrt(distance_2) distance *= distance_2 mag = dt / distance mass = MASS[j] mass_mag = mass * mag VX[i] -= dx * mass_mag VY[i] -= dy * mass_mag VZ[i] -= dz * mass_mag mass = MASS[i] mass_mag = mass * mag VX[j] += dx * mass_mag VY[j] += dy * mass_mag VZ[j] += dz * mass_mag end for X[i] += dt * VX[i] Y[i] += dt * VY[i] Z[i] += dt * VZ[i] end for end procedure -- advance function energy() atom e, dx, dy, dz, distance, mass_i e = 0.0 for i = 1 to NBODIES do mass_i = MASS[i] e += 0.5 * mass_i * (VX[i]*VX[i] + VY[i]*VY[i] + VZ[i]*VZ[i]) for j = i + 1 to NBODIES do dx = X[i] - X[j] dy = Y[i] - Y[j] dz = Z[i] - Z[j] distance = sqrt(dx*dx + dy*dy + dz*dz) e -= (mass_i*MASS[j])/distance end for end for return e end function -- energy procedure offset_momentum() atom px, py, pz px = 0.0 py = 0.0 pz = 0.0 for i = 1 to NBODIES do px += VX[i] * MASS[i] py += VY[i] * MASS[i] pz += VZ[i] * MASS[i] end for VX[1] = - px / SOLAR_MASS VY[1] = - py / SOLAR_MASS VZ[1] = - pz / SOLAR_MASS end procedure -- offset_momentum procedure main(sequence argv) object n if length(argv) > 2 then n = value(argv[3]) n = n[2] else n = 1000 end if offset_momentum() printf(1, "%.9f\n", energy()) for i = 1 to n do advance(0.01) end for printf(1, "%.9f\n", energy()) end procedure -- main main(command_line())
N-Body Euphoria (Alternate)
-- The Computer Language Shootout Benchmarks -- http://benchmarksgame.alioth.debian.org -- -- Converted to Euphoria by Jason Gade -- Version 6: vectorized position and velocity operations -- run: exu nbody.ex N -- with profile without warning without type_check include get.e constant PI = 3.141592653589793, SOLAR_MASS = 4 * PI * PI, DAYS_PER_YEAR = 365.24 -- point fields constant X = 1, Y = 2, Z = 3 -- body fields constant POS = 1, VELOCITY = 2, MASS = 3 sequence sun = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, SOLAR_MASS} sequence 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} sequence bodies = {sun, jupiter, saturn, uranus, neptune} procedure offset_momentum() for i = 1 to length(bodies) do bodies[1][VELOCITY] -= bodies[i][VELOCITY] * bodies[i][MASS] / SOLAR_MASS end for end procedure function energy() atom e = 0 sequence dx, vx for i = 1 to length(bodies) do vx = bodies[i][VELOCITY] * bodies[i][VELOCITY] e += 0.5 * bodies[i][MASS] * (vx[X] + vx[Y] + vx[Z]) for j = i+1 to length(bodies) do dx = bodies[i][POS] - bodies[j][POS] dx *= dx e -= bodies[i][MASS] * bodies[j][MASS] / sqrt(dx[X] + dx[Y] + dx[Z]) end for end for return e end function procedure advance(atom dt) sequence dx, dx2 atom dx2sum, distance, mag for i = 1 to length(bodies) do for j = i+1 to length(bodies) do dx = bodies[i][POS] - bodies[j][POS] dx2 = dx * dx dx2sum = dx2[X] + dx2[Y] + dx2[Z] distance = sqrt(dx2sum) mag = dt / (distance * dx2sum) bodies[i][VELOCITY] -= dx * bodies[j][MASS] * mag bodies[j][VELOCITY] += dx * bodies[i][MASS] * mag end for bodies[i][POS] += dt * bodies[i][VELOCITY] end for end procedure procedure main(sequence argv) object n if length(argv) > 2 then n = value(argv[3]) else n = {GET_FAIL, 0} end if if n[1] = GET_SUCCESS then n = n[2] else n = 1000 end if offset_momentum() printf(1, "%.9f\n", energy()) for i = 1 to n do advance(0.01) end for printf(1, "%.9f\n", energy()) end procedure main(command_line())
N-Body Python
# The Computer Language Benchmarks Game # http://shootout.alioth.debian.org/ # # originally by Kevin Carson # modified by Tupteq, Fredrik Johansson, and Daniel Nanz # modified by Maciej Fijalkowski # 2to3 import sys 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 PI = 3.14159265358979323 SOLAR_MASS = 4 * PI * PI DAYS_PER_YEAR = 365.24 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) } SYSTEM = list(BODIES.values()) PAIRS = combinations(SYSTEM) def advance(dt, n, bodies=SYSTEM, pairs=PAIRS): for i in range(n): for (([x1, y1, z1], v1, m1), ([x2, y2, z2], v2, m2)) in pairs: dx = x1 - x2 dy = y1 - y2 dz = z1 - z2 mag = dt * ((dx * dx + dy * dy + dz * dz) ** (-1.5)) b1m = m1 * mag b2m = m2 * mag v1[0] -= dx * b2m v1[1] -= dy * b2m v1[2] -= dz * b2m v2[0] += dx * b1m v2[1] += dy * b1m v2[2] += dz * b1m for (r, [vx, vy, vz], m) in bodies: r[0] += dt * vx r[1] += dt * vy r[2] += dt * vz def report_energy(bodies=SYSTEM, pairs=PAIRS, e=0.0): 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=SYSTEM, px=0.0, py=0.0, 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, ref='sun'): offset_momentum(BODIES[ref]) report_energy() advance(0.01, n) report_energy() if __name__ == '__main__': main(int(sys.argv[1]))
N-Body Perl
# The Computer Language Shootout # http://shootout.alioth.debian.org/ # # contributed by Christoph Bauer # converted into Perl by M?rton Papp # fixed and cleaned up by Danny Sauer # optimized by Jesse Millikan use constant PI => 3.141592653589793; use constant SOLAR_MASS => (4 * PI * PI); use constant DAYS_PER_YEAR => 365.24; # Globals for arrays... Oh well. # Almost every iteration is a range, so I keep the last index rather than a count. my (@xs, @ys, @zs, @vxs, @vys, @vzs, @mass, $last); sub advance($) { my ($dt) = @_; my ($mm, $mm2, $j, $dx, $dy, $dz, $distance, $mag); # This is faster in the outer loop... for (0..$last) { # But not in the inner loop. Strange. for ($j = $_ + 1; $j < $last + 1; $j++) { $dx = $xs[$_] - $xs[$j]; $dy = $ys[$_] - $ys[$j]; $dz = $zs[$_] - $zs[$j]; $distance = sqrt($dx * $dx + $dy * $dy + $dz * $dz); $mag = $dt / ($distance * $distance * $distance); $mm = $mass[$_] * $mag; $mm2 = $mass[$j] * $mag; $vxs[$_] -= $dx * $mm2; $vxs[$j] += $dx * $mm; $vys[$_] -= $dy * $mm2; $vys[$j] += $dy * $mm; $vzs[$_] -= $dz * $mm2; $vzs[$j] += $dz * $mm; } # We're done with planet $_ at this point # This could be done in a seperate loop, but it's slower $xs[$_] += $dt * $vxs[$_]; $ys[$_] += $dt * $vys[$_]; $zs[$_] += $dt * $vzs[$_]; } } sub energy { my ($e, $i, $dx, $dy, $dz, $distance); $e = 0.0; for $i (0..$last) { $e += 0.5 * $mass[$i] * ($vxs[$i] * $vxs[$i] + $vys[$i] * $vys[$i] + $vzs[$i] * $vzs[$i]); for ($i + 1..$last) { $dx = $xs[$i] - $xs[$_]; $dy = $ys[$i] - $ys[$_]; $dz = $zs[$i] - $zs[$_]; $distance = sqrt($dx * $dx + $dy * $dy + $dz * $dz); $e -= ($mass[$i] * $mass[$_]) / $distance; } } return $e; } sub offset_momentum { my ($px, $py, $pz) = (0.0, 0.0, 0.0); for (0..$last) { $px += $vxs[$_] * $mass[$_]; $py += $vys[$_] * $mass[$_]; $pz += $vzs[$_] * $mass[$_]; } $vxs[0] = - $px / SOLAR_MASS; $vys[0] = - $py / SOLAR_MASS; $vzs[0] = - $pz / SOLAR_MASS; } # @ns = ( sun, jupiter, saturn, uranus, neptune ) @xs = (0, 4.84143144246472090e+00, 8.34336671824457987e+00, 1.28943695621391310e+01, 1.53796971148509165e+01); @ys = (0, -1.16032004402742839e+00, 4.12479856412430479e+00, -1.51111514016986312e+01, -2.59193146099879641e+01); @zs = (0, -1.03622044471123109e-01, -4.03523417114321381e-01, -2.23307578892655734e-01, 1.79258772950371181e-01); @vxs = map {$_ * DAYS_PER_YEAR} (0, 1.66007664274403694e-03, -2.76742510726862411e-03, 2.96460137564761618e-03, 2.68067772490389322e-03); @vys = map {$_ * DAYS_PER_YEAR} (0, 7.69901118419740425e-03, 4.99852801234917238e-03, 2.37847173959480950e-03, 1.62824170038242295e-03); @vzs = map {$_ * DAYS_PER_YEAR} (0, -6.90460016972063023e-05, 2.30417297573763929e-05, -2.96589568540237556e-05, -9.51592254519715870e-05); @mass = map {$_ * SOLAR_MASS} (1, 9.54791938424326609e-04, 2.85885980666130812e-04, 4.36624404335156298e-05, 5.15138902046611451e-05); $last = @xs - 1; offset_momentum(); printf ("%.9f\n", energy()); my $n = $ARGV[0]; # This does not, in fact, consume N*4 bytes of memory for (1..$n){ advance(0.01); } printf ("%.9f\n", energy());