Goal is to figure out why nbodies in lua is 20x faster than Rune

benches/benches/nbodies.lua

@@ -0,0 +1,123 @@
+-- The Computer Language Benchmarks Game
+-- https://salsa.debian.org/benchmarksgame-team/benchmarksgame/
+-- contributed by Mike Pall
+-- modified by Geoff Leyland
+
+local sqrt = math.sqrt
+
+local PI = 3.141592653589793
+local SOLAR_MASS = 4 * PI * PI
+local DAYS_PER_YEAR = 365.24
+local bodies = {
+  { -- Sun
+    x = 0,
+    y = 0,
+    z = 0,
+    vx = 0,
+    vy = 0,
+    vz = 0,
+    mass = SOLAR_MASS
+  },
+  { -- Jupiter
+    x = 4.84143144246472090e+00,
+    y = -1.16032004402742839e+00,
+    z = -1.03622044471123109e-01,
+    vx = 1.66007664274403694e-03 * DAYS_PER_YEAR,
+    vy = 7.69901118419740425e-03 * DAYS_PER_YEAR,
+    vz = -6.90460016972063023e-05 * DAYS_PER_YEAR,
+    mass = 9.54791938424326609e-04 * SOLAR_MASS
+  },
+  { -- Saturn
+    x = 8.34336671824457987e+00,
+    y = 4.12479856412430479e+00,
+    z = -4.03523417114321381e-01,
+    vx = -2.76742510726862411e-03 * DAYS_PER_YEAR,
+    vy = 4.99852801234917238e-03 * DAYS_PER_YEAR,
+    vz = 2.30417297573763929e-05 * DAYS_PER_YEAR,
+    mass = 2.85885980666130812e-04 * SOLAR_MASS
+  },
+  { -- Uranus
+    x = 1.28943695621391310e+01,
+    y = -1.51111514016986312e+01,
+    z = -2.23307578892655734e-01,
+    vx = 2.96460137564761618e-03 * DAYS_PER_YEAR,
+    vy = 2.37847173959480950e-03 * DAYS_PER_YEAR,
+    vz = -2.96589568540237556e-05 * DAYS_PER_YEAR,
+    mass = 4.36624404335156298e-05 * SOLAR_MASS
+  },
+  { -- Neptune
+    x = 1.53796971148509165e+01,
+    y = -2.59193146099879641e+01,
+    z = 1.79258772950371181e-01,
+    vx = 2.68067772490389322e-03 * DAYS_PER_YEAR,
+    vy = 1.62824170038242295e-03 * DAYS_PER_YEAR,
+    vz = -9.51592254519715870e-05 * DAYS_PER_YEAR,
+    mass = 5.15138902046611451e-05 * SOLAR_MASS
+  }
+}
+
+local function advance(bodies, nbody, dt)
+  for i=1,nbody do
+    local bi = bodies[i]
+    local bix, biy, biz, bimass = bi.x, bi.y, bi.z, bi.mass
+    local bivx, bivy, bivz = bi.vx, bi.vy, bi.vz
+    for j=i+1,nbody do
+      local bj = bodies[j]
+      local dx, dy, dz = bix-bj.x, biy-bj.y, biz-bj.z
+      local mag = sqrt(dx*dx + dy*dy + dz*dz)
+      mag = dt / (mag * mag * mag)
+      local bm = bj.mass*mag
+      bivx = bivx - (dx * bm)
+      bivy = bivy - (dy * bm)
+      bivz = bivz - (dz * bm)
+      bm = bimass*mag
+      bj.vx = bj.vx + (dx * bm)
+      bj.vy = bj.vy + (dy * bm)
+      bj.vz = bj.vz + (dz * bm)
+    end
+    bi.vx = bivx
+    bi.vy = bivy
+    bi.vz = bivz
+    bi.x = bix + dt * bivx
+    bi.y = biy + dt * bivy
+    bi.z = biz + dt * bivz
+  end
+end
+
+local function energy(bodies, nbody)
+  local e = 0
+  for i=1,nbody do
+    local bi = bodies[i]
+    local vx, vy, vz, bim = bi.vx, bi.vy, bi.vz, bi.mass
+    e = e + (0.5 * bim * (vx*vx + vy*vy + vz*vz))
+    for j=i+1,nbody do
+      local bj = bodies[j]
+      local dx, dy, dz = bi.x-bj.x, bi.y-bj.y, bi.z-bj.z
+      local distance = sqrt(dx*dx + dy*dy + dz*dz)
+      e = e - ((bim * bj.mass) / distance)
+    end
+  end
+  return e
+end
+
+local function offsetMomentum(b, nbody)
+  local px, py, pz = 0, 0, 0
+  for i=1,nbody do
+    local bi = b[i]
+    local bim = bi.mass
+    px = px + (bi.vx * bim)
+    py = py + (bi.vy * bim)
+    pz = pz + (bi.vz * bim)
+  end
+  b[1].vx = -px / SOLAR_MASS
+  b[1].vy = -py / SOLAR_MASS
+  b[1].vz = -pz / SOLAR_MASS
+end
+
+local N = 20000
+local nbody = #bodies
+
+offsetMomentum(bodies, nbody)
+io.write( string.format("%0.9f",energy(bodies, nbody)), "\n")
+for i=1,N do advance(bodies, nbody, 0.01) end
+io.write( string.format("%0.9f",energy(bodies, nbody)), "\n")

benches/benches/nbodies.rn

@@ -0,0 +1,150 @@
+const PI = 3.141592653589793;
+const SOLAR_MASS = 4.0 * PI * PI;
+const DAYS_PER_YEAR = 365.24;
+
+fn advance(bodies, dt) {
+    for i in 0..bodies.len() {
+        let bi = bodies[i];
+        let bix = bi.x;
+        let biy = bi.y;
+        let biz = bi.z;
+        let bimass = bi.mass;
+        let bivx = bi.vx;
+        let bivy = bi.vy;
+        let bivz = bi.vz;
+
+        for j in i + 1..bodies.len() {
+            let bj = bodies[j];
+            let dx = bix - bj.x;
+            let dy = biy - bj.y;
+            let dz = biz - bj.z;
+            let mag = (dx * dx + dy * dy + dz * dz).sqrt();
+            mag = dt / (mag * mag * mag);
+            let bm = bj.mass * mag;
+            bivx = bivx - (dx * bm);
+            bivy = bivy - (dy * bm);
+            bivz = bivz - (dz * bm);
+            bm = bimass * mag;
+            bj.vx = bj.vx + (dx * bm);
+            bj.vy = bj.vy + (dy * bm);
+            bj.vz = bj.vz + (dz * bm);
+        }
+
+        bi.vx = bivx;
+        bi.vy = bivy;
+        bi.vz = bivz;
+        bi.x = bix + dt * bivx;
+        bi.y = biy + dt * bivy;
+        bi.z = biz + dt * bivz;
+    }
+}
+
+fn energy(bodies) {
+    let e = 0.0;
+
+    for i in 0..bodies.len() {
+        let bi = bodies[i];
+        let vx = bi.vx;
+        let vy = bi.vy;
+        let vz = bi.vz;
+        let bim = bi.mass;
+
+        e = e + (0.5 * bim * (vx * vx + vy * vy + vz * vz));
+
+        for j in i + 1..bodies.len() {
+            let bj = bodies[j];
+            let dx = bi.x - bj.x;
+            let dy = bi.y - bj.y;
+            let dz = bi.z - bj.z;
+            let distance = (dx * dx + dy * dy + dz * dz).sqrt();
+            e = e - ((bim * bj.mass) / distance);
+        }
+    }
+
+    e
+}
+
+fn offset_momentum(bodies) {
+    let px = 0.0;
+    let py = 0.0;
+    let pz = 0.0;
+
+    for i in 0..bodies.len() {
+        let bi = bodies[i];
+        let bim = bi.mass;
+        px = px + (bi.vx * bim);
+        py = py + (bi.vy * bim);
+        pz = pz + (bi.vz * bim);
+    }
+
+    bodies[0].vx = -px / SOLAR_MASS;
+    bodies[0].vy = -py / SOLAR_MASS;
+    bodies[0].vz = -pz / SOLAR_MASS;
+}
+
+pub fn main() {
+    let bodies = [
+        // Sun
+        #{
+            x: 0.0,
+            y: 0.0,
+            z: 0.0,
+            vx: 0.0,
+            vy: 0.0,
+            vz: 0.0,
+            mass: SOLAR_MASS,
+        },
+        // Jupiter
+        #{
+            x: 4.84143144246472090e+00,
+            y: -1.16032004402742839e+00,
+            z: -1.03622044471123109e-01,
+            vx: 1.66007664274403694e-03 * DAYS_PER_YEAR,
+            vy: 7.69901118419740425e-03 * DAYS_PER_YEAR,
+            vz: -6.90460016972063023e-05 * DAYS_PER_YEAR,
+            mass: 9.54791938424326609e-04 * SOLAR_MASS,
+        },
+        // Saturn
+        #{
+            x: 8.34336671824457987e+00,
+            y: 4.12479856412430479e+00,
+            z: -4.03523417114321381e-01,
+            vx: -2.76742510726862411e-03 * DAYS_PER_YEAR,
+            vy: 4.99852801234917238e-03 * DAYS_PER_YEAR,
+            vz: 2.30417297573763929e-05 * DAYS_PER_YEAR,
+            mass: 2.85885980666130812e-04 * SOLAR_MASS,
+        },
+        // Uranus
+        #{
+            x: 1.28943695621391310e+01,
+            y: -1.51111514016986312e+01,
+            z: -2.23307578892655734e-01,
+            vx: 2.96460137564761618e-03 * DAYS_PER_YEAR,
+            vy: 2.37847173959480950e-03 * DAYS_PER_YEAR,
+            vz: -2.96589568540237556e-05 * DAYS_PER_YEAR,
+            mass: 4.36624404335156298e-05 * SOLAR_MASS,
+        },
+        // Neptune
+        #{
+            x: 1.53796971148509165e+01,
+            y: -2.59193146099879641e+01,
+            z: 1.79258772950371181e-01,
+            vx: 2.68067772490389322e-03 * DAYS_PER_YEAR,
+            vy: 1.62824170038242295e-03 * DAYS_PER_YEAR,
+            vz: -9.51592254519715870e-05 * DAYS_PER_YEAR,
+            mass: 5.15138902046611451e-05 * SOLAR_MASS,
+        }
+    ];
+
+    let N = 20000;
+    let nbody = bodies.len();
+
+    offset_momentum(bodies);
+    println!("{}", energy(bodies));
+
+    for i in 1..N {
+        advance(bodies, 0.01);
+    }
+
+    println!("{}", energy(bodies));
+}