Wiki Diff NnBody, revision #2 to tip
== N-Body
See http://openeuphoria.org/forum/119406.wc
=== N-Body Euphoria
<eucode>
-- 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())
</eucode>
=== N-Body Euphoria (Alternate)
<eucode>
-- 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())
</eucode>
=== 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());
}}}
See http://openeuphoria.org/forum/119406.wc
=== N-Body Euphoria
<eucode>
-- 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())
</eucode>
=== N-Body Euphoria (Alternate)
<eucode>
-- 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())
</eucode>
=== 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());
}}}
