diff --git a/benchmark/nbody.gbs b/benchmark/nbody.gbs new file mode 100755 index 000000000..0660342a8 --- /dev/null +++ b/benchmark/nbody.gbs @@ -0,0 +1,205 @@ +#!/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 + + Public Sub Advance(dt As Float) + + 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 + S.Advance(0.01) + Next + +Next + +Print S.Energy() + diff --git a/benchmark/nbody.pl b/benchmark/nbody.pl new file mode 100755 index 000000000..c4c3f3ffb --- /dev/null +++ b/benchmark/nbody.pl @@ -0,0 +1,112 @@ +#!/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); + +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(); + +for (1..10) +{ + printf ("%.9f\n", energy()); + + # This does not, in fact, consume N*4 bytes of memory + for (1..100000){ + advance(0.01); + } + +} + +printf ("%.9f\n", energy()); diff --git a/benchmark/nbody.py b/benchmark/nbody.py new file mode 100755 index 000000000..f5924cf8a --- /dev/null +++ b/benchmark/nbody.py @@ -0,0 +1,120 @@ +#!/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) + + +def advance(dt, bodies=SYSTEM, pairs=PAIRS): + + 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): + advance(0.01) + + report_energy() + +if __name__ == '__main__': + main() \ No newline at end of file