Chris Collins

Lab · Playground

N-body: a small solar system, timed

The same five-body gravity simulation in fifteen languages, each timing itself — a hands-on look at how much the runtime matters when the program is all floating-point arithmetic in a tight loop.

Code blocks on this page are editable — change them and press Run.

A star and four planets, pulling on each other by Newton’s law of gravitation. Each step works out the force between every pair of bodies (ten pairs), nudges their velocities, then moves them along — the simplest integrator that keeps an orbit from spiralling away, "kick then drift" (semi-implicit Euler). The program prints the total energy of the system before and after 100 000 steps of 0.001 time units. Energy should be conserved, so the two numbers agreeing to six significant figures says the simulation is behaving; the last few digits drifting is the integrator’s error.

The physics is the excuse. The real subject is speed: this is all floating-point arithmetic in a tight loop, the thing language runtimes differ on most. Every block times its own simulation loop (not the start-up, not the runtime download) and prints elapsed:. Press Run on a few and compare — then raise STEPS and watch which ones keep up.

All fifteen print the same two energies to twelve decimal places, because they do the same IEEE 754 operations in the same order (the Mandelbrot page has more on why that matters). The units are chosen so that the gravitational constant is 1, and the star’s velocity is set at the start so the system’s total momentum is zero and it does not drift across the screen.

Indicative timings

One run of each block — the simulation loop only, as the block itself reports it — in headless Chromium 151 on an Apple M4 Max, with every runtime already downloaded. These are one run on one machine; your numbers will differ, and the ratios matter more than the milliseconds. Rust, Go and Kotlin are compiled and run on their official playgrounds' servers, not in your browser, so their times measure someone else’s computer.

Language Runs as Simulation loop

C++

Clang → WebAssembly

4 ms

JavaScript

V8’s JIT

6 ms

C

Clang → WebAssembly

7 ms

OCaml

js_of_ocaml → JavaScript

18 ms

Fortran

LFortran → WebAssembly

102 ms

Lua

Lua 5.4 interpreter (WebAssembly)

189 ms

Python

CPython 3.12 (Pyodide, WebAssembly)

462 ms

Clojure

Scittle interpreter (JavaScript)

733 ms

PHP

PHP 8.4 (WebAssembly)

1 401 ms

Perl

Perl 5 (WebAssembly)

1 458 ms

R

webR (WebAssembly)

3 530 ms

Ruby

CRuby 3.4 (ruby.wasm)

4 085 ms

Kotlin

JVM on api.kotlinlang.org

55 ms

Rust

debug build on play.rust-lang.org

73 ms

Go

go.dev playground

0 ms (its clock is fake — see below)

The C, C++ and JavaScript times are too short to rank against each other: at a few milliseconds, the timer and the first few hundred steps before the JIT warms up are most of it. Raise STEPS to 10 000 000 to separate them (and leave Ruby and R at 100 000).

Compiled to WebAssembly

C
#include <math.h>
#include <stdio.h>
#include <time.h>

#define STEPS 100000
#define DT 0.001
#define N 5

typedef struct { double x, y, z, vx, vy, vz, m; } Body;

/* A star and four planets, in units where G = 1. */
static Body bodies[N] = {
    {  0.0,  0.0,  0.0,   0.0,   0.0,  0.0,  1.0     },
    {  1.0,  0.0,  0.0,   0.0,   1.0,  0.01, 0.001   },
    {  0.0, -2.0,  0.05,  0.7,   0.0,  0.0,  0.0003  },
    { -3.5,  0.0,  0.0,   0.0,  -0.53, 0.02, 0.00005 },
    {  0.0,  5.0, -0.1,  -0.45,  0.0,  0.0,  0.00004 },
};

static void offset_momentum(void) {
    double px = 0, py = 0, pz = 0;
    for (int i = 0; i < N; i++) {
        px += bodies[i].m * bodies[i].vx;
        py += bodies[i].m * bodies[i].vy;
        pz += bodies[i].m * bodies[i].vz;
    }
    bodies[0].vx = -px / bodies[0].m;
    bodies[0].vy = -py / bodies[0].m;
    bodies[0].vz = -pz / bodies[0].m;
}

static double energy(void) {
    double e = 0;
    for (int i = 0; i < N; i++) {
        Body *a = &bodies[i];
        e += 0.5 * a->m * (a->vx * a->vx + a->vy * a->vy + a->vz * a->vz);
        for (int j = i + 1; j < N; j++) {
            Body *b = &bodies[j];
            double dx = a->x - b->x, dy = a->y - b->y, dz = a->z - b->z;
            e -= a->m * b->m / sqrt(dx * dx + dy * dy + dz * dz);
        }
    }
    return e;
}

static void advance(void) {
    for (int i = 0; i < N; i++) {
        Body *a = &bodies[i];
        for (int j = i + 1; j < N; j++) {
            Body *b = &bodies[j];
            double dx = a->x - b->x, dy = a->y - b->y, dz = a->z - b->z;
            double d2 = dx * dx + dy * dy + dz * dz;
            double mag = DT / (d2 * sqrt(d2));
            a->vx -= dx * b->m * mag; a->vy -= dy * b->m * mag; a->vz -= dz * b->m * mag;
            b->vx += dx * a->m * mag; b->vy += dy * a->m * mag; b->vz += dz * a->m * mag;
        }
    }
    for (int i = 0; i < N; i++) {
        bodies[i].x += DT * bodies[i].vx;
        bodies[i].y += DT * bodies[i].vy;
        bodies[i].z += DT * bodies[i].vz;
    }
}

static double now_ms(void) {
    struct timespec t;
    timespec_get(&t, TIME_UTC);  /* standard C11; clock_gettime is POSIX */
    return t.tv_sec * 1e3 + t.tv_nsec / 1e6;
}

int main(void) {
    offset_momentum();
    printf("bodies: %d, steps: %d\n", N, STEPS);
    printf("energy before: %.12f\n", energy());
    double t0 = now_ms();
    for (int s = 0; s < STEPS; s++) advance();
    double ms = now_ms() - t0;
    printf("energy after:  %.12f\n", energy());
    printf("elapsed: %.0f ms\n", ms);
    return 0;
}

The baseline: a fixed array of structs, pointers to the two bodies of each pair, and nothing allocated in the loop. Compiled by Clang at -O2 to WebAssembly, then run by the browser’s WebAssembly engine — which compiles it again, to machine code.

C++
#include <array>
#include <chrono>
#include <cmath>
#include <cstdio>

constexpr int STEPS = 100000;
constexpr double DT = 0.001;

struct Vec {
    double x, y, z;
    Vec operator-(const Vec &o) const { return {x - o.x, y - o.y, z - o.z}; }
    Vec operator*(double k) const { return {x * k, y * k, z * k}; }
    Vec &operator+=(const Vec &o) { x += o.x; y += o.y; z += o.z; return *this; }
    Vec &operator-=(const Vec &o) { x -= o.x; y -= o.y; z -= o.z; return *this; }
    double norm2() const { return x * x + y * y + z * z; }
};

struct Body { Vec pos, vel; double m; };

// A star and four planets, in units where G = 1.
std::array<Body, 5> bodies{{
    {{0, 0, 0}, {0, 0, 0}, 1.0},
    {{1, 0, 0}, {0, 1, 0.01}, 0.001},
    {{0, -2, 0.05}, {0.7, 0, 0}, 0.0003},
    {{-3.5, 0, 0}, {0, -0.53, 0.02}, 0.00005},
    {{0, 5, -0.1}, {-0.45, 0, 0}, 0.00004},
}};

double energy() {
    double e = 0;
    for (size_t i = 0; i < bodies.size(); ++i) {
        e += 0.5 * bodies[i].m * bodies[i].vel.norm2();
        for (size_t j = i + 1; j < bodies.size(); ++j)
            e -= bodies[i].m * bodies[j].m / std::sqrt((bodies[i].pos - bodies[j].pos).norm2());
    }
    return e;
}

void advance() {
    for (size_t i = 0; i < bodies.size(); ++i)
        for (size_t j = i + 1; j < bodies.size(); ++j) {
            Body &a = bodies[i], &b = bodies[j];
            Vec d = a.pos - b.pos;
            double d2 = d.norm2();
            double mag = DT / (d2 * std::sqrt(d2));
            a.vel -= d * b.m * mag;
            b.vel += d * a.m * mag;
        }
    for (Body &b : bodies) b.pos += b.vel * DT;
}

int main() {
    Vec p{0, 0, 0};
    for (const Body &b : bodies) p += b.vel * b.m;
    Body &star = bodies[0];
    star.vel = {-p.x / star.m, -p.y / star.m, -p.z / star.m};

    std::printf("bodies: %zu, steps: %d\n", bodies.size(), STEPS);
    std::printf("energy before: %.12f\n", energy());
    auto t0 = std::chrono::steady_clock::now();
    for (int s = 0; s < STEPS; ++s) advance();
    std::chrono::duration<double, std::milli> ms = std::chrono::steady_clock::now() - t0;
    std::printf("energy after:  %.12f\n", energy());
    std::printf("elapsed: %.0f ms\n", ms.count());
}

A small Vec type with operators makes the physics read like the formula, and costs nothing: the compiler inlines it all back into the arithmetic the C version spells out. (d * b.m * mag is (d * b.m) * mag, the same order as the C, which is why the energies match to the last digit.)

Fortran
program nbody
  implicit none
  integer, parameter :: n = 5, steps = 100000
  real(8), parameter :: dt = 0.001d0
  ! a star and four planets, in units where G = 1: position, velocity, mass
  real(8) :: pos(3, n), vel(3, n), m(n)
  real(8) :: p(3)
  integer :: i, s
  integer(8) :: t0, t1, rate

  pos = reshape([0d0, 0d0, 0d0,   1d0, 0d0, 0d0,   0d0, -2d0, 0.05d0, &
                 -3.5d0, 0d0, 0d0,   0d0, 5d0, -0.1d0], [3, n])
  vel = reshape([0d0, 0d0, 0d0,   0d0, 1d0, 0.01d0,   0.7d0, 0d0, 0d0, &
                 0d0, -0.53d0, 0.02d0,   -0.45d0, 0d0, 0d0], [3, n])
  m = [1d0, 0.001d0, 0.0003d0, 0.00005d0, 0.00004d0]

  p = 0
  do i = 1, n
    p = p + m(i) * vel(:, i)
  end do
  vel(:, 1) = -p / m(1)

  write (*, '(a, i0, a, i0)') 'bodies: ', n, ', steps: ', steps
  write (*, '(a, f15.12)') 'energy before: ', energy()
  call system_clock(t0, rate)
  do s = 1, steps
    call advance()
  end do
  call system_clock(t1)
  write (*, '(a, f15.12)') 'energy after:  ', energy()
  write (*, '(a, i0, a)') 'elapsed: ', (t1 - t0) * 1000 / rate, ' ms'

contains

  function energy() result(e)
    real(8) :: e, d(3)
    integer :: i, j
    e = 0
    do i = 1, n
      e = e + 0.5d0 * m(i) * (vel(1, i)**2 + vel(2, i)**2 + vel(3, i)**2)
      do j = i + 1, n
        d = pos(:, i) - pos(:, j)
        e = e - m(i) * m(j) / sqrt(d(1)**2 + d(2)**2 + d(3)**2)
      end do
    end do
  end function energy

  subroutine advance()
    real(8) :: d(3), d2, mag
    integer :: i, j
    do i = 1, n
      do j = i + 1, n
        d = pos(:, i) - pos(:, j)
        d2 = d(1)**2 + d(2)**2 + d(3)**2
        mag = dt / (d2 * sqrt(d2))
        vel(:, i) = vel(:, i) - d * m(j) * mag
        vel(:, j) = vel(:, j) + d * m(i) * mag
      end do
    end do
    pos = pos + dt * vel
  end subroutine advance

end program nbody

Fortran was made for exactly this. Positions and velocities are 3 × 5 arrays, one column per body (Fortran is column-major, so a body’s x, y and z sit next to each other in memory), and vel(:, i) - d * m(j) * mag is a whole-vector expression, as is pos = pos + dt * vel. LFortran compiles it through LLVM to WebAssembly in the page. It came out an order of magnitude behind Clang’s C here — LFortran is a much younger compiler, and the whole-array expressions on three-element vectors are a lot of small loops.

Just-in-time compiled

JavaScript
const STEPS = 100000;
const DT = 0.001;

// A star and four planets, in units where G = 1.
const bodies = [
  { x: 0, y: 0, z: 0, vx: 0, vy: 0, vz: 0, m: 1 },
  { x: 1, y: 0, z: 0, vx: 0, vy: 1, vz: 0.01, m: 0.001 },
  { x: 0, y: -2, z: 0.05, vx: 0.7, vy: 0, vz: 0, m: 0.0003 },
  { x: -3.5, y: 0, z: 0, vx: 0, vy: -0.53, vz: 0.02, m: 0.00005 },
  { x: 0, y: 5, z: -0.1, vx: -0.45, vy: 0, vz: 0, m: 0.00004 },
];

function offsetMomentum() {
  let px = 0, py = 0, pz = 0;
  for (const b of bodies) { px += b.m * b.vx; py += b.m * b.vy; pz += b.m * b.vz; }
  const star = bodies[0];
  star.vx = -px / star.m; star.vy = -py / star.m; star.vz = -pz / star.m;
}

function energy() {
  let e = 0;
  for (let i = 0; i < bodies.length; i++) {
    const a = bodies[i];
    e += 0.5 * a.m * (a.vx * a.vx + a.vy * a.vy + a.vz * a.vz);
    for (let j = i + 1; j < bodies.length; j++) {
      const b = bodies[j];
      const dx = a.x - b.x, dy = a.y - b.y, dz = a.z - b.z;
      e -= a.m * b.m / Math.sqrt(dx * dx + dy * dy + dz * dz);
    }
  }
  return e;
}

function advance() {
  for (let i = 0; i < bodies.length; i++) {
    const a = bodies[i];
    for (let j = i + 1; j < bodies.length; j++) {
      const b = bodies[j];
      const dx = a.x - b.x, dy = a.y - b.y, dz = a.z - b.z;
      const d2 = dx * dx + dy * dy + dz * dz;
      const mag = DT / (d2 * Math.sqrt(d2));
      a.vx -= dx * b.m * mag; a.vy -= dy * b.m * mag; a.vz -= dz * b.m * mag;
      b.vx += dx * a.m * mag; b.vy += dy * a.m * mag; b.vz += dz * a.m * mag;
    }
  }
  for (const b of bodies) { b.x += DT * b.vx; b.y += DT * b.vy; b.z += DT * b.vz; }
}

offsetMomentum();
console.log(`bodies: ${bodies.length}, steps: ${STEPS}`);
console.log(`energy before: ${energy().toFixed(12)}`);
const t0 = performance.now();
for (let s = 0; s < STEPS; s++) advance();
const ms = performance.now() - t0;
console.log(`energy after:  ${energy().toFixed(12)}`);
console.log(`elapsed: ${ms.toFixed(0)} ms`);

Plain objects with the same seven fields, which the JavaScript engine notices: after a few thousand steps V8 compiles advance to machine code specialised for objects of exactly this shape, with the doubles unboxed. Add a field to just one body, or store a number as a string, and watch the time go up.

OCaml
let steps = 100000
let dt = 0.001

type body = {
  mutable x : float; mutable y : float; mutable z : float;
  mutable vx : float; mutable vy : float; mutable vz : float;
  m : float;
}

let body x y z vx vy vz m = { x; y; z; vx; vy; vz; m }

(* A star and four planets, in units where G = 1. *)
let bodies = [|
  body 0. 0. 0. 0. 0. 0. 1.;
  body 1. 0. 0. 0. 1. 0.01 0.001;
  body 0. (-2.) 0.05 0.7 0. 0. 0.0003;
  body (-3.5) 0. 0. 0. (-0.53) 0.02 0.00005;
  body 0. 5. (-0.1) (-0.45) 0. 0. 0.00004;
|]
let n = Array.length bodies

let offset_momentum () =
  let px = ref 0. and py = ref 0. and pz = ref 0. in
  Array.iter (fun b ->
      px := !px +. b.m *. b.vx; py := !py +. b.m *. b.vy; pz := !pz +. b.m *. b.vz) bodies;
  let star = bodies.(0) in
  star.vx <- -. !px /. star.m; star.vy <- -. !py /. star.m; star.vz <- -. !pz /. star.m

let energy () =
  let e = ref 0. in
  for i = 0 to n - 1 do
    let a = bodies.(i) in
    e := !e +. 0.5 *. a.m *. (a.vx *. a.vx +. a.vy *. a.vy +. a.vz *. a.vz);
    for j = i + 1 to n - 1 do
      let b = bodies.(j) in
      let dx = a.x -. b.x and dy = a.y -. b.y and dz = a.z -. b.z in
      e := !e -. a.m *. b.m /. sqrt (dx *. dx +. dy *. dy +. dz *. dz)
    done
  done;
  !e

let advance () =
  for i = 0 to n - 1 do
    let a = bodies.(i) in
    for j = i + 1 to n - 1 do
      let b = bodies.(j) in
      let dx = a.x -. b.x and dy = a.y -. b.y and dz = a.z -. b.z in
      let d2 = dx *. dx +. dy *. dy +. dz *. dz in
      let mag = dt /. (d2 *. sqrt d2) in
      a.vx <- a.vx -. dx *. b.m *. mag; a.vy <- a.vy -. dy *. b.m *. mag; a.vz <- a.vz -. dz *. b.m *. mag;
      b.vx <- b.vx +. dx *. a.m *. mag; b.vy <- b.vy +. dy *. a.m *. mag; b.vz <- b.vz +. dz *. a.m *. mag
    done
  done;
  Array.iter (fun b -> b.x <- b.x +. dt *. b.vx; b.y <- b.y +. dt *. b.vy; b.z <- b.z +. dt *. b.vz) bodies

let () =
  offset_momentum ();
  Printf.printf "bodies: %d, steps: %d\n" n steps;
  Printf.printf "energy before: %.12f\n" (energy ());
  let t0 = Sys.time () in
  for _ = 1 to steps do advance () done;
  let ms = (Sys.time () -. t0) *. 1000. in
  Printf.printf "energy after:  %.12f\n" (energy ());
  Printf.printf "elapsed: %.0f ms\n" ms

A record with mutable fields — OCaml makes you say which fields can change, and ← is the only way to change them. Natively, OCaml compiles to fast machine code; here the toplevel translates each phrase to JavaScript (js_of_ocaml) and V8 compiles that, so it lands near the JavaScript version.

Clojure
(def steps 100000)
(def dt 0.001)

;; A star and four planets, in units where G = 1: [x y z vx vy vz m]
(def bodies
  (to-array-2d [[0 0 0 0 0 0 1]
                [1 0 0 0 1 0.01 0.001]
                [0 -2 0.05 0.7 0 0 0.0003]
                [-3.5 0 0 0 -0.53 0.02 0.00005]
                [0 5 -0.1 -0.45 0 0 0.00004]]))
(def n (alength bodies))

(defn offset-momentum! []
  (let [p (fn [k] (reduce + 0 (map #(* (aget % 6) (aget % k)) bodies)))
        star (aget bodies 0)]
    (doseq [k [3 4 5]] (aset star k (/ (- (p k)) (aget star 6))))))

(defn energy []
  (loop [i 0 e 0]
    (if (= i n)
      e
      (let [a (aget bodies i)
            e (+ e (* 0.5 (aget a 6) (+ (* (aget a 3) (aget a 3)) (* (aget a 4) (aget a 4)) (* (aget a 5) (aget a 5)))))]
        (recur (inc i)
               (loop [j (inc i) e e]
                 (if (= j n)
                   e
                   (let [b (aget bodies j)
                         dx (- (aget a 0) (aget b 0)) dy (- (aget a 1) (aget b 1)) dz (- (aget a 2) (aget b 2))]
                     (recur (inc j) (- e (/ (* (aget a 6) (aget b 6)) (Math/sqrt (+ (* dx dx) (* dy dy) (* dz dz))))))))))))))

(defn advance! []
  (dotimes [i n]
    (let [a (aget bodies i)]
      (doseq [j (range (inc i) n)]
        (let [b (aget bodies j)
              dx (- (aget a 0) (aget b 0)) dy (- (aget a 1) (aget b 1)) dz (- (aget a 2) (aget b 2))
              d2 (+ (* dx dx) (* dy dy) (* dz dz))
              mag (/ dt (* d2 (Math/sqrt d2)))]
          (doseq [[k d] [[3 dx] [4 dy] [5 dz]]]
            (aset a k (- (aget a k) (* d (aget b 6) mag)))
            (aset b k (+ (aget b k) (* d (aget a 6) mag))))))))
  (doseq [b bodies k [0 1 2]]
    (aset b k (+ (aget b k) (* dt (aget b (+ k 3)))))))

(offset-momentum!)
(println (str "bodies: " n ", steps: " steps))
(println (str "energy before: " (.toFixed (energy) 12)))
(let [t0 (js/performance.now)]
  (dotimes [_ steps] (advance!))
  (let [ms (- (js/performance.now) t0)]
    (println (str "energy after:  " (.toFixed (energy) 12)))
    (println (str "elapsed: " (.toFixed ms 0) " ms"))))

Clojure’s own data is immutable, so a simulation that changes five bodies 100 000 times reaches for mutable Java(Script) arrays — to-array-2d, aget, aset — which is what idiomatic performance-minded Clojure does too. The catch here is the runtime: Scittle is an interpreter for Clojure written in ClojureScript, so every aget is a function call walked by the interpreter. Even so it beat several of the interpreters below, because Scittle itself runs on V8’s JIT.

Interpreters compiled to WebAssembly

Python
import time
from math import sqrt

STEPS = 100_000
DT = 0.001

# A star and four planets, in units where G = 1: [x, y, z, vx, vy, vz, m]
bodies = [
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0],
    [1.0, 0.0, 0.0, 0.0, 1.0, 0.01, 0.001],
    [0.0, -2.0, 0.05, 0.7, 0.0, 0.0, 0.0003],
    [-3.5, 0.0, 0.0, 0.0, -0.53, 0.02, 0.00005],
    [0.0, 5.0, -0.1, -0.45, 0.0, 0.0, 0.00004],
]
pairs = [(bodies[i], bodies[j]) for i in range(len(bodies)) for j in range(i + 1, len(bodies))]


def offset_momentum():
    px = py = pz = 0.0
    for _, _, _, vx, vy, vz, m in bodies:
        px += m * vx
        py += m * vy
        pz += m * vz
    star = bodies[0]
    star[3], star[4], star[5] = -px / star[6], -py / star[6], -pz / star[6]


def energy():
    e = 0.0
    for i, (x, y, z, vx, vy, vz, m) in enumerate(bodies):
        e += 0.5 * m * (vx * vx + vy * vy + vz * vz)
        for x2, y2, z2, _, _, _, m2 in bodies[i + 1:]:
            dx, dy, dz = x - x2, y - y2, z - z2
            e -= m * m2 / sqrt(dx * dx + dy * dy + dz * dz)
    return e


def advance(steps):
    for _ in range(steps):
        for a, b in pairs:
            dx, dy, dz = a[0] - b[0], a[1] - b[1], a[2] - b[2]
            d2 = dx * dx + dy * dy + dz * dz
            mag = DT / (d2 * sqrt(d2))
            a[3] -= dx * b[6] * mag
            a[4] -= dy * b[6] * mag
            a[5] -= dz * b[6] * mag
            b[3] += dx * a[6] * mag
            b[4] += dy * a[6] * mag
            b[5] += dz * a[6] * mag
        for b in bodies:
            b[0] += DT * b[3]
            b[1] += DT * b[4]
            b[2] += DT * b[5]


offset_momentum()
print(f"bodies: {len(bodies)}, steps: {STEPS}")
print(f"energy before: {energy():.12f}")
t0 = time.perf_counter()
advance(STEPS)
ms = (time.perf_counter() - t0) * 1000
print(f"energy after:  {energy():.12f}")
print(f"elapsed: {ms:.0f} ms")

Each body is a plain list and the ten pairs are worked out once, up front, so the inner loop does no index arithmetic. Every + is still a dynamically typed operation on boxed float objects — which is why this is one or two orders of magnitude behind the compiled versions, and why numerical Python hands loops like this to NumPy, or to C.

Ruby
STEPS = 100_000
DT = 0.001

Body = Struct.new(:x, :y, :z, :vx, :vy, :vz, :m)

# A star and four planets, in units where G = 1.
BODIES = [
  Body.new(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0),
  Body.new(1.0, 0.0, 0.0, 0.0, 1.0, 0.01, 0.001),
  Body.new(0.0, -2.0, 0.05, 0.7, 0.0, 0.0, 0.0003),
  Body.new(-3.5, 0.0, 0.0, 0.0, -0.53, 0.02, 0.00005),
  Body.new(0.0, 5.0, -0.1, -0.45, 0.0, 0.0, 0.00004),
]
PAIRS = BODIES.combination(2).to_a

def offset_momentum
  px = py = pz = 0.0
  BODIES.each { |b| px += b.m * b.vx; py += b.m * b.vy; pz += b.m * b.vz }
  star = BODIES[0]
  star.vx, star.vy, star.vz = -px / star.m, -py / star.m, -pz / star.m
end

def energy
  e = 0.0
  BODIES.each_with_index do |a, i|
    e += 0.5 * a.m * (a.vx * a.vx + a.vy * a.vy + a.vz * a.vz)
    BODIES[(i + 1)..].each do |b|
      dx, dy, dz = a.x - b.x, a.y - b.y, a.z - b.z
      e -= a.m * b.m / Math.sqrt(dx * dx + dy * dy + dz * dz)
    end
  end
  e
end

def advance
  PAIRS.each do |a, b|
    dx, dy, dz = a.x - b.x, a.y - b.y, a.z - b.z
    d2 = dx * dx + dy * dy + dz * dz
    mag = DT / (d2 * Math.sqrt(d2))
    a.vx -= dx * b.m * mag; a.vy -= dy * b.m * mag; a.vz -= dz * b.m * mag
    b.vx += dx * a.m * mag; b.vy += dy * a.m * mag; b.vz += dz * a.m * mag
  end
  BODIES.each { |b| b.x += DT * b.vx; b.y += DT * b.vy; b.z += DT * b.vz }
end

offset_momentum
puts "bodies: #{BODIES.size}, steps: #{STEPS}"
puts format("energy before: %.12f", energy)
t0 = Process.clock_gettime(Process::CLOCK_MONOTONIC)
STEPS.times { advance }
ms = (Process.clock_gettime(Process::CLOCK_MONOTONIC) - t0) * 1000
puts format("energy after:  %.12f", energy)
puts format("elapsed: %.0f ms", ms)

A Struct gives named fields without writing a class, and combination(2) produces the ten pairs in the same order as the nested loops. Ruby has a JIT (YJIT) natively, but not in WebAssembly, so this is the bytecode interpreter.

Lua
local STEPS, DT = 100000, 0.001
local sqrt = math.sqrt

-- A star and four planets, in units where G = 1.
local bodies = {
  { x = 0, y = 0, z = 0, vx = 0, vy = 0, vz = 0, m = 1 },
  { x = 1, y = 0, z = 0, vx = 0, vy = 1, vz = 0.01, m = 0.001 },
  { x = 0, y = -2, z = 0.05, vx = 0.7, vy = 0, vz = 0, m = 0.0003 },
  { x = -3.5, y = 0, z = 0, vx = 0, vy = -0.53, vz = 0.02, m = 0.00005 },
  { x = 0, y = 5, z = -0.1, vx = -0.45, vy = 0, vz = 0, m = 0.00004 },
}
local n = #bodies

local function offset_momentum()
  local px, py, pz = 0.0, 0.0, 0.0
  for _, b in ipairs(bodies) do
    px = px + b.m * b.vx; py = py + b.m * b.vy; pz = pz + b.m * b.vz
  end
  local star = bodies[1]
  star.vx, star.vy, star.vz = -px / star.m, -py / star.m, -pz / star.m
end

local function energy()
  local e = 0.0
  for i = 1, n do
    local a = bodies[i]
    e = e + 0.5 * a.m * (a.vx * a.vx + a.vy * a.vy + a.vz * a.vz)
    for j = i + 1, n do
      local b = bodies[j]
      local dx, dy, dz = a.x - b.x, a.y - b.y, a.z - b.z
      e = e - a.m * b.m / sqrt(dx * dx + dy * dy + dz * dz)
    end
  end
  return e
end

local function advance()
  for i = 1, n do
    local a = bodies[i]
    for j = i + 1, n do
      local b = bodies[j]
      local dx, dy, dz = a.x - b.x, a.y - b.y, a.z - b.z
      local d2 = dx * dx + dy * dy + dz * dz
      local mag = DT / (d2 * sqrt(d2))
      a.vx = a.vx - dx * b.m * mag; a.vy = a.vy - dy * b.m * mag; a.vz = a.vz - dz * b.m * mag
      b.vx = b.vx + dx * a.m * mag; b.vy = b.vy + dy * a.m * mag; b.vz = b.vz + dz * a.m * mag
    end
  end
  for i = 1, n do
    local b = bodies[i]
    b.x = b.x + DT * b.vx; b.y = b.y + DT * b.vy; b.z = b.z + DT * b.vz
  end
end

offset_momentum()
print(string.format("bodies: %d, steps: %d", n, STEPS))
print(string.format("energy before: %.12f", energy()))
local t0 = os.clock()
for _ = 1, STEPS do advance() end
local ms = (os.clock() - t0) * 1000
print(string.format("energy after:  %.12f", energy()))
print(string.format("elapsed: %.0f ms", ms))

Lua’s one data structure, the table, does duty as the struct, and local sqrt = math.sqrt saves a table lookup on every call — the classic Lua speed tip. Lua’s interpreter is small and quick; LuaJIT would be another order of magnitude faster, but it cannot run in WebAssembly.

PHP
<?php
const STEPS = 100000;
const DT = 0.001;

// A star and four planets, in units where G = 1: [x, y, z, vx, vy, vz, m]
$bodies = [
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0],
    [1.0, 0.0, 0.0, 0.0, 1.0, 0.01, 0.001],
    [0.0, -2.0, 0.05, 0.7, 0.0, 0.0, 0.0003],
    [-3.5, 0.0, 0.0, 0.0, -0.53, 0.02, 0.00005],
    [0.0, 5.0, -0.1, -0.45, 0.0, 0.0, 0.00004],
];
$n = count($bodies);

function offset_momentum(array &$bodies): void {
    $px = $py = $pz = 0.0;
    foreach ($bodies as $b) {
        $px += $b[6] * $b[3]; $py += $b[6] * $b[4]; $pz += $b[6] * $b[5];
    }
    $bodies[0][3] = -$px / $bodies[0][6];
    $bodies[0][4] = -$py / $bodies[0][6];
    $bodies[0][5] = -$pz / $bodies[0][6];
}

function energy(array $bodies): float {
    $e = 0.0;
    $n = count($bodies);
    for ($i = 0; $i < $n; $i++) {
        [$x, $y, $z, $vx, $vy, $vz, $m] = $bodies[$i];
        $e += 0.5 * $m * ($vx * $vx + $vy * $vy + $vz * $vz);
        for ($j = $i + 1; $j < $n; $j++) {
            $dx = $x - $bodies[$j][0]; $dy = $y - $bodies[$j][1]; $dz = $z - $bodies[$j][2];
            $e -= $m * $bodies[$j][6] / sqrt($dx * $dx + $dy * $dy + $dz * $dz);
        }
    }
    return $e;
}

function advance(array &$bodies, int $n): void {
    for ($i = 0; $i < $n; $i++) {
        $a = &$bodies[$i];
        for ($j = $i + 1; $j < $n; $j++) {
            $b = &$bodies[$j];
            $dx = $a[0] - $b[0]; $dy = $a[1] - $b[1]; $dz = $a[2] - $b[2];
            $d2 = $dx * $dx + $dy * $dy + $dz * $dz;
            $mag = DT / ($d2 * sqrt($d2));
            $a[3] -= $dx * $b[6] * $mag; $a[4] -= $dy * $b[6] * $mag; $a[5] -= $dz * $b[6] * $mag;
            $b[3] += $dx * $a[6] * $mag; $b[4] += $dy * $a[6] * $mag; $b[5] += $dz * $a[6] * $mag;
            unset($b);
        }
        unset($a);
    }
    foreach ($bodies as &$b) {
        $b[0] += DT * $b[3]; $b[1] += DT * $b[4]; $b[2] += DT * $b[5];
    }
}

offset_momentum($bodies);
printf("bodies: %d, steps: %d\n", $n, STEPS);
printf("energy before: %.12f\n", energy($bodies));
$t0 = hrtime(true);
for ($s = 0; $s < STEPS; $s++) advance($bodies, $n);
$ms = (hrtime(true) - $t0) / 1e6;
printf("energy after:  %.12f\n", energy($bodies));
printf("elapsed: %.0f ms\n", $ms);

PHP arrays are values: passing one to a function copies it (lazily), so the functions that change the bodies take &$bodies by reference, and $a = &$bodies[$i] aliases a body rather than copying it. The unset calls break those aliases before the variable is reused — forgetting one is a classic PHP bug.

Perl
use strict;
use warnings;
use Time::HiRes qw(time);

my $STEPS = 100_000;
my $DT    = 0.001;

# A star and four planets, in units where G = 1: [x, y, z, vx, vy, vz, m]
my @bodies = (
    [0, 0, 0, 0, 0, 0, 1],
    [1, 0, 0, 0, 1, 0.01, 0.001],
    [0, -2, 0.05, 0.7, 0, 0, 0.0003],
    [-3.5, 0, 0, 0, -0.53, 0.02, 0.00005],
    [0, 5, -0.1, -0.45, 0, 0, 0.00004],
);
my @pairs = map { my $i = $_; map { [$bodies[$i], $bodies[$_]] } $i + 1 .. $#bodies } 0 .. $#bodies;

sub offset_momentum {
    my ($px, $py, $pz) = (0, 0, 0);
    for my $b (@bodies) {
        $px += $b->[6] * $b->[3]; $py += $b->[6] * $b->[4]; $pz += $b->[6] * $b->[5];
    }
    my $star = $bodies[0];
    @$star[3, 4, 5] = (-$px / $star->[6], -$py / $star->[6], -$pz / $star->[6]);
}

sub energy {
    my $e = 0;
    for my $i (0 .. $#bodies) {
        my ($x, $y, $z, $vx, $vy, $vz, $m) = @{ $bodies[$i] };
        $e += 0.5 * $m * ($vx * $vx + $vy * $vy + $vz * $vz);
        for my $b (@bodies[$i + 1 .. $#bodies]) {
            my ($dx, $dy, $dz) = ($x - $b->[0], $y - $b->[1], $z - $b->[2]);
            $e -= $m * $b->[6] / sqrt($dx * $dx + $dy * $dy + $dz * $dz);
        }
    }
    return $e;
}

sub advance {
    for my $p (@pairs) {
        my ($a, $b) = @$p;
        my ($dx, $dy, $dz) = ($a->[0] - $b->[0], $a->[1] - $b->[1], $a->[2] - $b->[2]);
        my $d2 = $dx * $dx + $dy * $dy + $dz * $dz;
        my $mag = $DT / ($d2 * sqrt($d2));
        $a->[3] -= $dx * $b->[6] * $mag; $a->[4] -= $dy * $b->[6] * $mag; $a->[5] -= $dz * $b->[6] * $mag;
        $b->[3] += $dx * $a->[6] * $mag; $b->[4] += $dy * $a->[6] * $mag; $b->[5] += $dz * $a->[6] * $mag;
    }
    for my $b (@bodies) {
        $b->[0] += $DT * $b->[3]; $b->[1] += $DT * $b->[4]; $b->[2] += $DT * $b->[5];
    }
}

offset_momentum();
printf "bodies: %d, steps: %d\n", scalar @bodies, $STEPS;
printf "energy before: %.12f\n", energy();
my $t0 = time;
advance() for 1 .. $STEPS;
my $ms = (time - $t0) * 1000;
printf "energy after:  %.12f\n", energy();
printf "elapsed: %.0f ms\n", $ms;

Array references ($b→[3]) for each body, a list of pairs built with nested map, and an array slice assignment (@$star[3, 4, 5] = …) for the momentum offset. Time::HiRes replaces the built-in time, which only counts whole seconds.

R
steps <- 100000
dt <- 0.001

# A star and four planets, in units where G = 1: one row per body
pos <- rbind(c(0, 0, 0), c(1, 0, 0), c(0, -2, 0.05), c(-3.5, 0, 0), c(0, 5, -0.1))
vel <- rbind(c(0, 0, 0), c(0, 1, 0.01), c(0.7, 0, 0), c(0, -0.53, 0.02), c(-0.45, 0, 0))
m <- c(1, 0.001, 0.0003, 0.00005, 0.00004)
n <- length(m)

p <- c(0, 0, 0)
for (i in 1:n) p <- p + m[i] * vel[i, ]
vel[1, ] <- -p / m[1]

energy <- function() {
  e <- 0
  for (i in 1:n) {
    e <- e + 0.5 * m[i] * (vel[i, 1] * vel[i, 1] + vel[i, 2] * vel[i, 2] + vel[i, 3] * vel[i, 3])
    for (j in seq_len(n - i) + i) {
      d <- pos[i, ] - pos[j, ]
      e <- e - m[i] * m[j] / sqrt(d[1] * d[1] + d[2] * d[2] + d[3] * d[3])
    }
  }
  e
}

advance <- function() {
  for (i in 1:(n - 1)) {
    for (j in (i + 1):n) {
      d <- pos[i, ] - pos[j, ]
      d2 <- d[1] * d[1] + d[2] * d[2] + d[3] * d[3]
      mag <- dt / (d2 * sqrt(d2))
      vel[i, ] <<- vel[i, ] - d * m[j] * mag
      vel[j, ] <<- vel[j, ] + d * m[i] * mag
    }
  }
  pos <<- pos + dt * vel
}

cat(sprintf("bodies: %d, steps: %d\n", n, steps))
cat(sprintf("energy before: %.12f\n", energy()))
t0 <- proc.time()[["elapsed"]]
for (s in seq_len(steps)) advance()
ms <- (proc.time()[["elapsed"]] - t0) * 1000
cat(sprintf("energy after:  %.12f\n", energy()))
cat(sprintf("elapsed: %.0f ms\n", ms))

R is vectorised, so each body’s position and velocity is a row of a matrix and the update works on three components at once. But the pairs loop is a loop, every vel[i, ] <← copies a row back into a global matrix, and R’s interpreter pays heavily for both. It is the language’s weak spot, and the reason serious R simulations are written with whole-matrix operations or in C++ via Rcpp.

Compiled on a server

Rust
use std::time::Instant;

const STEPS: usize = 100_000;
const DT: f64 = 0.001;

#[derive(Clone, Copy)]
struct Body { pos: [f64; 3], vel: [f64; 3], m: f64 }

fn energy(bodies: &[Body]) -> f64 {
    let mut e = 0.0;
    for (i, a) in bodies.iter().enumerate() {
        e += 0.5 * a.m * (a.vel[0] * a.vel[0] + a.vel[1] * a.vel[1] + a.vel[2] * a.vel[2]);
        for b in &bodies[i + 1..] {
            let d = [a.pos[0] - b.pos[0], a.pos[1] - b.pos[1], a.pos[2] - b.pos[2]];
            e -= a.m * b.m / (d[0] * d[0] + d[1] * d[1] + d[2] * d[2]).sqrt();
        }
    }
    e
}

fn advance(bodies: &mut [Body]) {
    for i in 0..bodies.len() {
        let (head, tail) = bodies.split_at_mut(i + 1);
        let a = &mut head[i];
        for b in tail.iter_mut() {
            let d = [a.pos[0] - b.pos[0], a.pos[1] - b.pos[1], a.pos[2] - b.pos[2]];
            let d2 = d[0] * d[0] + d[1] * d[1] + d[2] * d[2];
            let mag = DT / (d2 * d2.sqrt());
            for k in 0..3 {
                a.vel[k] -= d[k] * b.m * mag;
                b.vel[k] += d[k] * a.m * mag;
            }
        }
    }
    for b in bodies.iter_mut() {
        for k in 0..3 {
            b.pos[k] += DT * b.vel[k];
        }
    }
}

fn main() {
    // A star and four planets, in units where G = 1.
    let mut bodies = [
        Body { pos: [0.0, 0.0, 0.0], vel: [0.0, 0.0, 0.0], m: 1.0 },
        Body { pos: [1.0, 0.0, 0.0], vel: [0.0, 1.0, 0.01], m: 0.001 },
        Body { pos: [0.0, -2.0, 0.05], vel: [0.7, 0.0, 0.0], m: 0.0003 },
        Body { pos: [-3.5, 0.0, 0.0], vel: [0.0, -0.53, 0.02], m: 0.00005 },
        Body { pos: [0.0, 5.0, -0.1], vel: [-0.45, 0.0, 0.0], m: 0.00004 },
    ];
    let mut p = [0.0; 3];
    for b in &bodies {
        for k in 0..3 {
            p[k] += b.m * b.vel[k];
        }
    }
    let star = &mut bodies[0];
    for k in 0..3 {
        star.vel[k] = -p[k] / star.m;
    }

    println!("bodies: {}, steps: {}", bodies.len(), STEPS);
    println!("energy before: {:.12}", energy(&bodies));
    let t0 = Instant::now();
    for _ in 0..STEPS {
        advance(&mut bodies);
    }
    let ms = t0.elapsed().as_secs_f64() * 1000.0;
    println!("energy after:  {:.12}", energy(&bodies));
    println!("elapsed: {ms:.0} ms");
}

The borrow checker will not let you hold two &mut references into one array, so split_at_mut hands out the body i and the bodies after it as two separate mutable slices — Rust’s way of proving the pair are different bodies. The playground builds in debug mode, so this time is several times what cargo run --release would give.

Go
package main

import (
	"fmt"
	"math"
	"time"
)

const steps = 100000
const dt = 0.001

type body struct{ x, y, z, vx, vy, vz, m float64 }

func energy(bodies []body) float64 {
	e := 0.0
	for i := range bodies {
		a := &bodies[i]
		e += 0.5 * a.m * (a.vx*a.vx + a.vy*a.vy + a.vz*a.vz)
		for j := i + 1; j < len(bodies); j++ {
			b := &bodies[j]
			dx, dy, dz := a.x-b.x, a.y-b.y, a.z-b.z
			e -= a.m * b.m / math.Sqrt(dx*dx+dy*dy+dz*dz)
		}
	}
	return e
}

func advance(bodies []body) {
	for i := range bodies {
		a := &bodies[i]
		for j := i + 1; j < len(bodies); j++ {
			b := &bodies[j]
			dx, dy, dz := a.x-b.x, a.y-b.y, a.z-b.z
			d2 := dx*dx + dy*dy + dz*dz
			mag := dt / (d2 * math.Sqrt(d2))
			a.vx -= dx * b.m * mag
			a.vy -= dy * b.m * mag
			a.vz -= dz * b.m * mag
			b.vx += dx * a.m * mag
			b.vy += dy * a.m * mag
			b.vz += dz * a.m * mag
		}
	}
	for i := range bodies {
		b := &bodies[i]
		b.x += dt * b.vx
		b.y += dt * b.vy
		b.z += dt * b.vz
	}
}

func main() {
	// A star and four planets, in units where G = 1.
	bodies := []body{
		{0, 0, 0, 0, 0, 0, 1},
		{1, 0, 0, 0, 1, 0.01, 0.001},
		{0, -2, 0.05, 0.7, 0, 0, 0.0003},
		{-3.5, 0, 0, 0, -0.53, 0.02, 0.00005},
		{0, 5, -0.1, -0.45, 0, 0, 0.00004},
	}
	var px, py, pz float64
	for _, b := range bodies {
		px += b.m * b.vx
		py += b.m * b.vy
		pz += b.m * b.vz
	}
	bodies[0].vx, bodies[0].vy, bodies[0].vz = -px/bodies[0].m, -py/bodies[0].m, -pz/bodies[0].m

	fmt.Printf("bodies: %d, steps: %d\n", len(bodies), steps)
	fmt.Printf("energy before: %.12f\n", energy(bodies))
	t0 := time.Now()
	for s := 0; s < steps; s++ {
		advance(bodies)
	}
	ms := float64(time.Since(t0).Microseconds()) / 1000
	fmt.Printf("energy after:  %.12f\n", energy(bodies))
	fmt.Printf("elapsed: %.0f ms\n", ms)
}

Pointers into a slice (a := &bodies[i]) update the bodies in place; ranging over bodies by value would update copies. (On an ARM machine, Go may fuse a.vx - dx*b.m*mag into one multiply-add instruction, and the last digits of the energy can then differ — the Mandelbrot page shows the float64(…) conversions that prevent it. The playground’s servers are x86, where it does not happen.) One surprise: the Go playground runs programs with a fake clock that only moves when the program sleeps, so elapsed always says 0 ms there. Run it on your own machine for a real number.

Kotlin
import java.util.Locale
import kotlin.math.sqrt

const val STEPS = 100_000
const val DT = 0.001

class Body(var x: Double, var y: Double, var z: Double,
           var vx: Double, var vy: Double, var vz: Double, val m: Double)

// A star and four planets, in units where G = 1.
val bodies = arrayOf(
    Body(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0),
    Body(1.0, 0.0, 0.0, 0.0, 1.0, 0.01, 0.001),
    Body(0.0, -2.0, 0.05, 0.7, 0.0, 0.0, 0.0003),
    Body(-3.5, 0.0, 0.0, 0.0, -0.53, 0.02, 0.00005),
    Body(0.0, 5.0, -0.1, -0.45, 0.0, 0.0, 0.00004),
)

fun energy(): Double {
    var e = 0.0
    for (i in bodies.indices) {
        val a = bodies[i]
        e += 0.5 * a.m * (a.vx * a.vx + a.vy * a.vy + a.vz * a.vz)
        for (j in i + 1 until bodies.size) {
            val b = bodies[j]
            val dx = a.x - b.x; val dy = a.y - b.y; val dz = a.z - b.z
            e -= a.m * b.m / sqrt(dx * dx + dy * dy + dz * dz)
        }
    }
    return e
}

fun advance() {
    for (i in bodies.indices) {
        val a = bodies[i]
        for (j in i + 1 until bodies.size) {
            val b = bodies[j]
            val dx = a.x - b.x; val dy = a.y - b.y; val dz = a.z - b.z
            val d2 = dx * dx + dy * dy + dz * dz
            val mag = DT / (d2 * sqrt(d2))
            a.vx -= dx * b.m * mag; a.vy -= dy * b.m * mag; a.vz -= dz * b.m * mag
            b.vx += dx * a.m * mag; b.vy += dy * a.m * mag; b.vz += dz * a.m * mag
        }
    }
    for (b in bodies) { b.x += DT * b.vx; b.y += DT * b.vy; b.z += DT * b.vz }
}

fun f12(x: Double) = String.format(Locale.ROOT, "%.12f", x)

fun main() {
    val star = bodies[0]
    star.vx = -bodies.sumOf { it.m * it.vx } / star.m
    star.vy = -bodies.sumOf { it.m * it.vy } / star.m
    star.vz = -bodies.sumOf { it.m * it.vz } / star.m

    println("bodies: ${bodies.size}, steps: $STEPS")
    println("energy before: ${f12(energy())}")
    val t0 = System.nanoTime()
    repeat(STEPS) { advance() }
    val ms = (System.nanoTime() - t0) / 1e6
    println("energy after:  ${f12(energy())}")
    println("elapsed: ${ms.toLong()} ms")
}

A class with var properties, sumOf for the momentum, and repeat for the loop. The JVM starts slowly and then compiles the hot loop to machine code (HotSpot’s JIT), so a run this short measures partly the warm-up. Locale.ROOT keeps the decimal point a point on a server set up for, say, Germany.

What it shows

The ranking holds no surprises in its shape — compiled code first, JIT-compiled close behind, interpreters a long way back — but the size of the gaps is the lesson. For a loop of pure arithmetic over a handful of objects:

  • WebAssembly from C, C++ and Fortran is close to native speed. The browser compiles the WebAssembly to machine code, and there is nothing dynamic left to check.

  • JavaScript is within a small factor of that, because V8 sees the same object shape every time and generates specialised machine code. OCaml, translated to JavaScript, rides on the same JIT.

  • Interpreters (Python, Ruby, Lua, PHP, Perl, R) pay for every operation: look up the type, unbox the float, do the sum, box the result. Lua’s small interpreter is the quick one, Python is next; R, built for whole vectors rather than scalar loops, and Ruby, whose JIT (YJIT) cannot run in WebAssembly, are the slow ones — 500 to 1000 times the C.

  • An interpreter written in JavaScript (Clojure in Scittle) does better than you might expect, because V8 compiles the interpreter itself to machine code. The same Clojure compiled for the JVM would be far faster again.

And the energies agree to twelve places in every language, which is the less obvious lesson: the answer does not depend on the language, only the time it takes to get it.