2.0 3.0 > doc > benchmark > nbody
Previous  Next  Edit  Rename  Undo  Search  Administration
 Documentation
N-Body
This benchmark performs an N-body simulation of the Jovian planets. It is repeated ten times.

## Results

The execution time is the user time added to the system time, as returned by the bash "time" command.

Python Perl Gambas Gambas + JIT
Execution time (s) 22.4 33.3 20.8 3.07
vs. Python 1 1.44 0.93 0.14
vs. Perl 0.69 1 0.65 0.09
vs. Gambas 1.07 1.55 1 0.62
vs. Gambas + JIT 7.30 10.85 6.78 1

## Python source code

```#!/usr/bin/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)

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():
offset_momentum(BODIES['sun'])
for t in range(10):
report_energy()
for n in range(100000):

report_energy()

if __name__ == '__main__':
main()
```

## Perl source code

```#!/usr/bin/perl -w

# 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);

{
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();

for (1..10)
{
printf ("%.9f\n", energy());

# This does not, in fact, consume N*4 bytes of memory
for (1..100000){
}

}

printf ("%.9f\n", energy());
```

## Gambas source code

```#!/usr/bin/env gbs3

Class Body

Static Private SOLAR_MASS As Float = 4 * Pi * Pi
Private Const DAYS_PER_YEAR As Float = 365.24

Public X As Float
Public Y As Float
Public Z As Float
Public VX As Float
Public VY As Float
Public VZ As Float
Public Mass As Float

Static Public Sub Jupiter() As Body

Dim P As New Body
p.x = 4.84143144246472090e+00
p.y = -1.16032004402742839e+00
p.z = -1.03622044471123109e-01
p.vx = 1.66007664274403694e-03 * DAYS_PER_YEAR
p.vy = 7.69901118419740425e-03 * DAYS_PER_YEAR
p.vz = -6.90460016972063023e-05 * DAYS_PER_YEAR
p.mass = 9.54791938424326609e-04 * SOLAR_MASS
return p

End

Static Public Sub Saturn() As Body

Dim P As New Body
p.x = 8.34336671824457987e+00
p.y = 4.12479856412430479e+00
p.z = -4.03523417114321381e-01
p.vx = -2.76742510726862411e-03 * DAYS_PER_YEAR
p.vy = 4.99852801234917238e-03 * DAYS_PER_YEAR
p.vz = 2.30417297573763929e-05 * DAYS_PER_YEAR
p.mass = 2.85885980666130812e-04 * SOLAR_MASS
return p

End

Static Public Sub Uranus() As Body

Dim P As New Body

p.x = 1.28943695621391310e+01
p.y = -1.51111514016986312e+01
p.z = -2.23307578892655734e-01
p.vx = 2.96460137564761618e-03 * DAYS_PER_YEAR
p.vy = 2.37847173959480950e-03 * DAYS_PER_YEAR
p.vz = -2.96589568540237556e-05 * DAYS_PER_YEAR
p.mass = 4.36624404335156298e-05 * SOLAR_MASS
return p

End

Static Public Sub Neptune() As Body

Dim P As New Body

p.x = 1.53796971148509165e+01
p.y = -2.59193146099879641e+01
p.z = 1.79258772950371181e-01
p.vx = 2.68067772490389322e-03 * DAYS_PER_YEAR
p.vy = 1.62824170038242295e-03 * DAYS_PER_YEAR
p.vz = -9.51592254519715870e-05 * DAYS_PER_YEAR
p.mass = 5.15138902046611451e-05 * SOLAR_MASS
return p

End

Static Public Sub Sun() As Body

Dim P As New Body

p.mass = SOLAR_MASS
return p

End

Public Sub OffsetMomentum(px As Float, py As Float, pz As Float) As Body

vx = -px / SOLAR_MASS
vy = -py / SOLAR_MASS
vz = -pz / SOLAR_MASS

Return Me

End

End Class

Class NBodySystem

Private Bodies As Body[]

Public Sub _new()

Dim PX, PY, PZ As Float
Dim I As Integer

Bodies = [ Body.Sun(), Body.Jupiter(), Body.Saturn(), Body.Uranus(), Body.Neptune() ]

For I = 0 To Bodies.Max
PX += Bodies[I].vx * Bodies[I].mass
PY += Bodies[I].vy * Bodies[I].mass
PZ += Bodies[I].vz * Bodies[I].mass
Next

Bodies[0].offsetMomentum(PX, PY, PZ)

End

Dim I, J As Integer
Dim iBody, jBody As Body
Dim dx, dy, dz As Float
Dim dSquared, fMag As Float
Dim iMass, jMass, iMag, jMag As Float

For I = 0 To Bodies.Max
iBody = Bodies[I]
iMass = iBody.mass

For J = I + 1 To Bodies.Max
jBody = Bodies[J]
jMass = jBody.mass

dx = iBody.x - jBody.x
dy = iBody.y - jBody.y
dz = iBody.z - jBody.z

dSquared = dx * dx + dy * dy + dz * dz
fMag = dt / (dSquared * Sqr(dSquared))
iMag = iMass * fMag
jMag = jMass * fMag

iBody.vx -= dx * jMag
iBody.vy -= dy * jMag
iBody.vz -= dz * jMag

jBody.vx += dx * iMag
jBody.vy += dy * iMag
jBody.vz += dz * iMag
Next
Next

For Each iBody in Bodies
iBody.x += dt * iBody.vx
iBody.y += dt * iBody.vy
iBody.z += dt * iBody.vz
Next

End

Public Sub Energy() As Float

Dim dx, dy, dz, distance, E As Float
Dim iBody, jBody As Body
Dim I, J As Integer

For I = 0 To Bodies.Max

iBody = bodies[i]
E += 0.5 * iBody.mass * (iBody.vx * iBody.vx + iBody.vy * iBody.vy + iBody.vz * iBody.vz)

For J = I + 1 To Bodies.Max

jBody = Bodies[J]
dx = iBody.x - jBody.x
dy = iBody.y - jBody.y
dz = iBody.z - jBody.z

distance = Sqr(dx*dx + dy*dy + dz*dz)
E -= (iBody.mass * jBody.mass) / distance

Next
Next

return E

End

End Class

Dim S As New NBodySystem
Dim I, N As Integer

For N = 1 To 10

Print S.Energy()

For I = 1 To 100000