TypedPython benchmarks: Cython (pure, typed) and numba columns

.gitignore

@@ -58,3 +58,10 @@ google-services.json
 # Local Claude Code settings: per machine, not part of the repository.
 .claude/
 .cache/
+
+# TypedPython benchmark harness: uv venv, Cython/numba build output and numba caches
+benchmarks/typedpython/.venv/
+benchmarks/typedpython/build/
+benchmarks/typedpython/target/
+benchmarks/typedpython/java-out/
+benchmarks/typedpython/tp-build/

benchmarks/typedpython/cython/binary_trees.pyx

@@ -0,0 +1,46 @@
+# cython: language_level=3, boundscheck=False, wraparound=False, cdivision=True, initializedcheck=False
+# Same algorithm as py/binary_trees.py: one heap object per node, reference counted, freed by
+# refcounting (not an arena). Semantics that differ from CPython: Node is a cdef class (no __dict__,
+# attributes not visible from Python), C long counters, recursion uses the C stack (no
+# RecursionError), no None checks on typed attributes.
+import sys
+
+
+cdef class Node:
+    cdef Node left
+    cdef Node right
+
+    def __init__(self, Node left, Node right):
+        self.left = left
+        self.right = right
+
+
+cdef Node make(int depth):
+    if depth > 0:
+        return Node(make(depth - 1), make(depth - 1))
+    return Node(None, None)
+
+
+cdef long check(Node node) noexcept:
+    cdef Node left = node.left
+    cdef Node right = node.right
+    if left is None or right is None:
+        return 1
+    return 1 + check(left) + check(right)
+
+
+def main():
+    cdef int n = int(sys.argv[1])
+    cdef int min_depth = 4
+    cdef int max_depth = max(min_depth + 2, n)
+    cdef int depth, i
+    cdef long iterations, chk, total = 0
+    cdef long stretch = check(make(max_depth + 1))
+    cdef Node long_lived = make(max_depth)
+    for depth in range(min_depth, max_depth + 1, 2):
+        iterations = 1 << (max_depth - depth + min_depth)
+        chk = 0
+        for i in range(iterations):
+            chk += check(make(depth))
+        total += chk
+    print(f"{stretch} {total} {check(long_lived)}")

benchmarks/typedpython/cython/build.py

@@ -0,0 +1,46 @@
+"""Build the Cython columns. Run with the harness venv's python (see ../README in run.py's docstring).
+
+    .venv/bin/python cython/build.py <out_dir>
+
+Writes <out_dir>/pure/<name>.so (the unmodified py/*.py compiled in pure Python mode, annotations
+as written) and <out_dir>/typed/<name>_typed.so (cython/*.pyx). Not part of any timed run.
+Flags: -O3 -ffp-contract=off (no fused multiply-add, so float results match CPython bit for bit).
+"""
+import shutil
+import sys
+from pathlib import Path
+
+from Cython.Build import cythonize
+from setuptools import Extension, setup
+
+HERE = Path(__file__).resolve().parent
+PY = HERE.parent / "py"
+NAMES = {"nbody": "nbody", "spectral_norm": "spectral_norm", "fannkuch": "fannkuch",
+         "binary_trees": "binary_trees", "wordfreq": "wordfreq"}
+CFLAGS = ["-O3", "-ffp-contract=off"]
+
+
+def build(out: Path, kind: str) -> None:
+    work = out / f"src-{kind}"
+    shutil.rmtree(work, ignore_errors=True)
+    work.mkdir(parents=True)
+    exts = []
+    for name in NAMES:
+        if kind == "pure":
+            src = work / f"{name}.py"
+            shutil.copy(PY / f"{name}.py", src)       # unmodified copy
+            mod = name
+        else:
+            mod = f"{name}_typed"
+            src = work / f"{mod}.pyx"
+            shutil.copy(HERE / f"{name}.pyx", src)
+        exts.append(Extension(mod, [str(src)], extra_compile_args=CFLAGS))
+    ext_modules = cythonize(exts, language_level=3, quiet=True, build_dir=str(work / "c"))
+    setup(name=f"tp-bench-{kind}", ext_modules=ext_modules, script_name="setup.py",
+          script_args=["build_ext", "--build-lib", str(out / kind), "--build-temp", str(work / "obj"), "-q"])
+
+
+if __name__ == "__main__":
+    target = Path(sys.argv[1]).resolve()
+    for k in ("pure", "typed"):
+        build(target, k)

benchmarks/typedpython/cython/fannkuch.pyx

@@ -0,0 +1,61 @@
+# cython: language_level=3, boundscheck=False, wraparound=False, cdivision=True, initializedcheck=False
+# Same algorithm as py/fannkuch.py. Semantics that differ from CPython: C int arrays and C long
+# counters (wrap instead of growing), no bounds checks, fixed capacity of 16 elements (n > 16 raises).
+import sys
+
+
+cdef (long, long) fannkuch(int n) noexcept:
+    cdef int perm1[16]
+    cdef int count[16]
+    cdef int perm[16]
+    cdef long max_flips = 0, checksum = 0, perm_count = 0, flips
+    cdef int r = n, i, k, lo, hi, t, p0
+    for i in range(n):
+        perm1[i] = i
+        count[i] = 0
+    while True:
+        while r != 1:
+            count[r - 1] = r
+            r -= 1
+        for i in range(n):
+            perm[i] = perm1[i]
+        flips = 0
+        k = perm[0]
+        while k != 0:
+            lo = 0
+            hi = k
+            while lo < hi:
+                t = perm[lo]
+                perm[lo] = perm[hi]
+                perm[hi] = t
+                lo += 1
+                hi -= 1
+            flips += 1
+            k = perm[0]
+        if flips > max_flips:
+            max_flips = flips
+        if perm_count % 2 == 0:
+            checksum += flips
+        else:
+            checksum -= flips
+        while True:
+            if r == n:
+                return checksum, max_flips
+            p0 = perm1[0]
+            for i in range(r):
+                perm1[i] = perm1[i + 1]
+            perm1[r] = p0
+            count[r] -= 1
+            if count[r] > 0:
+                break
+            r += 1
+        perm_count += 1
+
+
+def main():
+    cdef int n = int(sys.argv[1])
+    if n > 16:
+        raise ValueError("fannkuch.pyx supports n <= 16")
+    cdef long checksum, max_flips
+    checksum, max_flips = fannkuch(n)
+    print(f"{checksum} {max_flips}")

benchmarks/typedpython/cython/nbody.pyx

@@ -0,0 +1,80 @@
+# cython: language_level=3, boundscheck=False, wraparound=False, cdivision=True, initializedcheck=False
+# Same algorithm as py/nbody.py. Semantics that differ from CPython: float division does not raise
+# ZeroDivisionError (cdivision), no bounds checks (arrays are C arrays), loop counters are C long.
+from libc.math cimport sqrt
+import sys
+
+cdef double PI = 3.141592653589793
+cdef double SOLAR_MASS = 4 * PI * PI
+cdef double DAYS_PER_YEAR = 365.24
+
+
+cdef void offset_momentum(double* vx, double* vy, double* vz, double* mass) noexcept nogil:
+    cdef double px = 0.0, py = 0.0, pz = 0.0
+    cdef int i
+    for i in range(5):
+        px += vx[i] * mass[i]
+        py += vy[i] * mass[i]
+        pz += vz[i] * mass[i]
+    vx[0] = -px / SOLAR_MASS
+    vy[0] = -py / SOLAR_MASS
+    vz[0] = -pz / SOLAR_MASS
+
+
+cdef void advance(double dt, long steps, double* x, double* y, double* z,
+                  double* vx, double* vy, double* vz, double* mass) noexcept nogil:
+    cdef long s
+    cdef int i, j
+    cdef double dx, dy, dz, d2, mag
+    for s in range(steps):
+        for i in range(5):
+            for j in range(i + 1, 5):
+                dx = x[i] - x[j]
+                dy = y[i] - y[j]
+                dz = z[i] - z[j]
+                d2 = dx * dx + dy * dy + dz * dz
+                mag = dt / (d2 * sqrt(d2))
+                vx[i] -= dx * mass[j] * mag
+                vy[i] -= dy * mass[j] * mag
+                vz[i] -= dz * mass[j] * mag
+                vx[j] += dx * mass[i] * mag
+                vy[j] += dy * mass[i] * mag
+                vz[j] += dz * mass[i] * mag
+        for i in range(5):
+            x[i] += dt * vx[i]
+            y[i] += dt * vy[i]
+            z[i] += dt * vz[i]
+
+
+cdef double energy(double* x, double* y, double* z, double* vx, double* vy, double* vz,
+                   double* mass) noexcept nogil:
+    cdef double e = 0.0, dx, dy, dz
+    cdef int i, j
+    for i in range(5):
+        e += 0.5 * mass[i] * (vx[i] * vx[i] + vy[i] * vy[i] + vz[i] * vz[i])
+        for j in range(i + 1, 5):
+            dx = x[i] - x[j]
+            dy = y[i] - y[j]
+            dz = z[i] - z[j]
+            e -= mass[i] * mass[j] / sqrt(dx * dx + dy * dy + dz * dz)
+    return e
+
+
+def main():
+    cdef long n = int(sys.argv[1])
+    cdef double[5] x = [0.0, 4.84143144246472090e+00, 8.34336671824457987e+00, 1.28943695621391310e+01, 1.53796971148509165e+01]
+    cdef double[5] y = [0.0, -1.16032004402742839e+00, 4.12479856412430479e+00, -1.51111514016986312e+01, -2.59193146099879641e+01]
+    cdef double[5] z = [0.0, -1.03622044471123109e-01, -4.03523417114321381e-01, -2.23307578892655734e-01, 1.79258772950371181e-01]
+    cdef double[5] vx = [0.0, 1.66007664274403694e-03 * DAYS_PER_YEAR, -2.76742510726862411e-03 * DAYS_PER_YEAR,
+                         2.96460137564761618e-03 * DAYS_PER_YEAR, 2.68067772490389322e-03 * DAYS_PER_YEAR]
+    cdef double[5] vy = [0.0, 7.69901118419740425e-03 * DAYS_PER_YEAR, 4.99852801234917238e-03 * DAYS_PER_YEAR,
+                         2.37847173959480950e-03 * DAYS_PER_YEAR, 1.62824170038242295e-03 * DAYS_PER_YEAR]
+    cdef double[5] vz = [0.0, -6.90460016972063023e-05 * DAYS_PER_YEAR, 2.30417297573763929e-05 * DAYS_PER_YEAR,
+                         -2.96589568540237556e-05 * DAYS_PER_YEAR, -9.51592254519715870e-05 * DAYS_PER_YEAR]
+    cdef double[5] mass = [SOLAR_MASS, 9.54791938424326609e-04 * SOLAR_MASS, 2.85885980666130812e-04 * SOLAR_MASS,
+                           4.36624404335156298e-05 * SOLAR_MASS, 5.15138902046611451e-05 * SOLAR_MASS]
+    offset_momentum(vx, vy, vz, mass)
+    cdef double before = energy(x, y, z, vx, vy, vz, mass)
+    advance(0.01, n, x, y, z, vx, vy, vz, mass)
+    cdef double after = energy(x, y, z, vx, vy, vz, mass)
+    print(f"{before:.9f} {after:.9f}")

benchmarks/typedpython/cython/spectral_norm.pyx

@@ -0,0 +1,61 @@
+# cython: language_level=3, boundscheck=False, wraparound=False, cdivision=True, initializedcheck=False
+# Same algorithm as py/spectral_norm.py. Semantics that differ from CPython: loop and index
+# arithmetic is C long (wraps instead of growing; not reachable at these sizes), cdivision (no
+# ZeroDivisionError), no bounds checks (malloc'd C arrays).
+from libc.math cimport sqrt
+from libc.stdlib cimport malloc, free
+import sys
+
+
+cdef inline double eval_a(long i, long j) noexcept nogil:
+    return 1.0 / ((i + j) * (i + j + 1) // 2 + i + 1)
+
+
+cdef void mul_av(long n, double* v, double* out) noexcept nogil:
+    cdef long i, j
+    cdef double s
+    for i in range(n):
+        s = 0.0
+        for j in range(n):
+            s += eval_a(i, j) * v[j]
+        out[i] = s
+
+
+cdef void mul_atv(long n, double* v, double* out) noexcept nogil:
+    cdef long i, j
+    cdef double s
+    for i in range(n):
+        s = 0.0
+        for j in range(n):
+            s += eval_a(j, i) * v[j]
+        out[i] = s
+
+
+cdef void mul_atav(long n, double* v, double* out, double* tmp) noexcept nogil:
+    mul_av(n, v, tmp)
+    mul_atv(n, tmp, out)
+
+
+def main():
+    cdef long n = int(sys.argv[1])
+    cdef long i
+    cdef double* u = <double*> malloc(n * sizeof(double))
+    cdef double* v = <double*> malloc(n * sizeof(double))
+    cdef double* tmp = <double*> malloc(n * sizeof(double))
+    cdef double vbv = 0.0, vv = 0.0
+    try:
+        for i in range(n):
+            u[i] = 1.0
+            v[i] = 0.0
+            tmp[i] = 0.0
+        for i in range(10):
+            mul_atav(n, u, v, tmp)
+            mul_atav(n, v, u, tmp)
+        for i in range(n):
+            vbv += u[i] * v[i]
+            vv += v[i] * v[i]
+    finally:
+        free(u)
+        free(v)
+        free(tmp)
+    print(f"{sqrt(vbv / vv):.9f}")

benchmarks/typedpython/cython/wordfreq.pyx

@@ -0,0 +1,54 @@
+# cython: language_level=3, boundscheck=False, wraparound=False, cdivision=True, initializedcheck=False
+# Same algorithm as py/wordfreq.py: the counts stay a Python dict[str, int] (the algorithm's data
+# structure), only the arithmetic is typed and words are built in a C buffer. Semantics that differ
+# from CPython: C long arithmetic (wraps; the LCG product stays below 2**47 and the weighted sum
+# below 2**63 at benchmark size), cdivision.
+import sys
+
+cdef long VOCAB = 50000
+cdef long MOD = 2147483647
+
+
+cdef inline long next_state(long s) noexcept nogil:
+    return (s * 48271) % MOD
+
+
+cdef str word_of(long word_id):
+    cdef char buf[32]
+    cdef int pos = 32
+    cdef long x = word_id + 17576
+    while x > 0:
+        pos -= 1
+        buf[pos] = <char>(97 + x % 26)
+        x //= 26
+    return (<bytes>buf[pos:32]).decode("ascii")
+
+
+def main():
+    cdef long n = int(sys.argv[1])
+    cdef dict counts = {}
+    cdef long s = 12345, a, b, i, c, weighted = 0, best_count
+    cdef str word, w, best_word
+    cdef list taken = []
+    cdef list parts = []
+    for i in range(n):
+        s = next_state(s)
+        a = s % VOCAB
+        s = next_state(s)
+        b = s % VOCAB
+        word = word_of(a * b // VOCAB)
+        counts[word] = <long>counts.get(word, 0) + 1
+    for w, c in counts.items():
+        weighted += c * c
+    for i in range(5):
+        best_word = ""
+        best_count = -1
+        for w, c in counts.items():
+            if w in taken:
+                continue
+            if c > best_count or (c == best_count and w < best_word):
+                best_word = w
+                best_count = c
+        taken.append(best_word)
+        parts.append(f"{best_word}:{best_count}")
+    print(f"{len(counts)} {weighted} {','.join(parts)}")