Algorithms

Eight compute kernels measured against the same eight in C, and six judge problems timed the way a judge times them.

algorithms.ppy holds eight kernels: sieve, Collatz, knapsack, edit distance, Floyd–Warshall, matmul, union-find, and Fermat. They run against algorithms.c, compiled with both gcc and clang. 15a–15f are six competitive-programming problems, each in its own folder. Each reads from standard input and is timed as a whole process, startup included.

Run it

python  algorithms.ppy
ppy run algorithms.ppy
gcc   -O3 algorithms.c -o algorithms_c     -lm && ./algorithms_c
clang -O3 algorithms.c -o algorithms_clang -lm && ./algorithms_clang

What it prints

python algorithms.ppy

sieve 2e6              195.0 ms   -> 148933
collatz 3e5           1209.1 ms   -> 442
knapsack 400x2e4       480.0 ms   -> 199600
edit 2000x2000         515.0 ms   -> 1846
floyd 220              473.4 ms   -> 558837
matmul 220             505.6 ms   -> 18883
union-find 5e5         155.6 ms   -> 250000
fermat 6e4              26.3 ms   -> 6114

ppy run algorithms.ppy

sieve 2e6                9.9 ms   -> 148933
collatz 3e5             41.8 ms   -> 442
knapsack 400x2e4         5.1 ms   -> 199600
edit 2000x2000           3.7 ms   -> 1846
floyd 220                3.2 ms   -> 558837
matmul 220               4.0 ms   -> 18883
union-find 5e5           3.3 ms   -> 250000
fermat 6e4               2.7 ms   -> 6114

gcc -O3 algorithms.c -o algorithms_c -lm && ./algorithms_c

sieve 2e6               15.9 ms   -> 148933
collatz 3e5             43.6 ms   -> 442
knapsack 400x2e4         2.6 ms   -> 199600
edit 2000x2000           3.8 ms   -> 1846
floyd 220                5.8 ms   -> 558837
matmul 220               2.5 ms   -> 18883
union-find 5e5           3.7 ms   -> 250000
fermat 6e4               1.8 ms   -> 6114

clang -O3 algorithms.c -o algorithms_clang -lm && ./algorithms_clang

sieve 2e6               18.7 ms   -> 148933
collatz 3e5             40.2 ms   -> 442
knapsack 400x2e4         2.0 ms   -> 199600
edit 2000x2000           3.6 ms   -> 1846
floyd 220                3.0 ms   -> 558837
matmul 220               3.9 ms   -> 18883
union-find 5e5           3.9 ms   -> 250000
fermat 6e4               1.5 ms   -> 6114

Each subfolder's README has its own commands, with a small input.txt so they run as written. bench.py generates the judge-sized inputs.

The eight kernels against C

Same machine, same session, kernel wall time: mean ± standard deviation over 7 runs, each a fresh process. Every row prints the same answer in every column.

kernel	plain	ppy run	ppy build --unsafe	C (gcc -O3)	C (clang -O3)
sieve 2e6	173.1 ± 2.9	9.2 ± 0.5	8.8 ± 0.4	13.9 ± 0.6	13.8 ± 0.3
collatz 3e5	1111.0 ± 13.1	39.9 ± 0.8	30.3 ± 0.4	41.4 ± 0.5	32.2 ± 0.3
knapsack 400×2e4	423.1 ± 2.7	5.0 ± 0.1	3.6 ± 0.2	2.4 ± 0.1	2.1 ± 0.3
edit 2000×2000	471.4 ± 6.5	2.9 ± 0.1	2.9 ± 0.3	3.6 ± 0.1	3.5 ± 0.2
floyd 220	436.3 ± 5.3	3.4 ± 0.5	3.1 ± 0.2	5.0 ± 0.1	2.8 ± 0.0
matmul 220	473.9 ± 5.7	3.6 ± 0.4	3.5 ± 0.1	2.1 ± 0.1	3.4 ± 0.1
union-find 5e5	149.6 ± 3.8	3.2 ± 0.2	3.2 ± 0.4	3.6 ± 0.2	3.8 ± 0.3
fermat 6e4	23.9 ± 0.3	2.3 ± 0.1	2.6 ± 0.2	1.8 ± 0.1	1.5 ± 0.1

Bold is the fastest cell in the row. The native path is 9× to 160× faster than plain CPython on these eight, and all three paths print identical answers.

The two PPy columns differ in integer semantics:

• ppy run keeps Python-integer semantics. Overflow is guarded and falls back to arbitrary precision.
• ppy build --unsafe is the wrap-semantics artifact. That is where collatz picks up its remaining 10 ms and knapsack its 1.4 ms. The other kernels are already guard-free in the loop and do not move.

--host-cpu is within the noise on all eight.

Guard hoisting ([tool.ppy.llvm] safeguards) closed most of the old gaps. A multiplied index like i * n + k proves its extreme cases once in a guard block ahead of the loop. The body then runs plain mul nsw with no side exits, which lets LLVM strength-reduce and vectorize.

Against both compilers:

• PPy wins sieve, edit distance, and union-find, and collatz once the artifact wraps.
• floyd beats gcc and loses to clang by a third of a millisecond.
• gcc and clang keep knapsack, matmul, and fermat, whose remaining guards live on data values no range can prove.

gcc 13.3 and clang 22.1, both -O3, no -march. Writes go through borrowed buffers, so the caller sees them.

The eight kernels against Numba, Mojo, and Codon

The same eight kernels are ported to three tools that compile Python-shaped code, in compare/: algorithms_numba.py, algorithms.mojo, algorithms_codon.py. The table shows kernel wall time inside the program, in milliseconds, mean and spread over seven fresh processes. Numba is timed after one warm call, so its compile is not in the number.

Here is Collatz as each of them spells it.

PPy is annotated Python, and ppy run keeps Python's integers: the multiply is overflow-checked and falls back to arbitrary precision.

def collatz_longest(limit: int) -> int:
    best: int = 0
    for start in range(1, limit):
        n: int = start
        steps: int = 0
        while n != 1:
            n = n // 2 if n % 2 == 0 else 3 * n + 1
            steps += 1
        best = max(best, steps)
    return best

Numba is the same function with @njit in place of the annotations and a NumPy array in place of array.array. Codon is the same source with List[int] annotations, run with codon run -release. Both wrap at 64 bits without a word.

@njit
def collatz_longest(limit):
    best = 0
    for start in range(1, limit):
        n = start
        ...

Mojo is typed by hand: Int for every index, Int64 for every element, var on every local. This is the plain List version without UnsafePointer:

def collatz_longest(limit: Int) -> Int64:
    var best: Int64 = 0
    for start in range(1, limit):
        var n = Int64(start)
        var steps: Int64 = 0
        while n != 1:
            n = n // 2 if n % 2 == 0 else 3 * n + 1
            steps += 1
        best = max(best, steps)
    return best

	ppy run	ppy build --unsafe	Numba @njit	Mojo	Codon
sieve 2e6	9.14 ± 0.67	9.04 ± 0.67	8.04 ± 0.25	11.42 ± 0.41	8.28 ± 0.38
collatz 3e5	41.56 ± 1.13	31.40 ± 0.55	33.14 ± 0.73	70.72 ± 0.97	33.38 ± 0.56
knapsack 400x2e4	5.08 ± 0.04	3.66 ± 0.09	4.48 ± 0.08	2.26 ± 0.11	5.90 ± 0.10
edit 2000x2000	3.02 ± 0.16	2.86 ± 0.05	1.76 ± 0.05	8.32 ± 0.20	2.46 ± 0.05
floyd 220	3.22 ± 0.04	3.10 ± 0.00	2.48 ± 0.08	3.15 ± 0.18	1.30 ± 0.34
matmul 220	3.64 ± 0.05	3.72 ± 0.08	3.48 ± 0.11	8.25 ± 0.16	3.90 ± 0.10
union-find 5e5	3.10 ± 0.07	3.18 ± 0.41	2.46 ± 0.05	6.68 ± 0.51	3.98 ± 0.32
fermat 6e4	2.36 ± 0.05	2.62 ± 0.04	1.74 ± 0.05	2.20 ± 0.03	1.52 ± 0.04

Read ppy build --unsafe against the other three: it is the wrap-semantics column. ppy run carries the overflow guards, which is where collatz pays its 10 ms.

• Numba leads on six of the eight by ten to thirty percent. It uses the same LLVM, a for over a NumPy array, and no guard of any kind.
• Mojo's plain List loop is the slowest on five kernels and the fastest on knapsack, where its bounds-free table indexing shows.

Numba 0.67.0 on CPython 3.12.13, Mojo 1.0.0 (-O3), Codon 0.19.6, PPy on CPython 3.13.13; Intel Core Ultra 9 386H.

The six problems

Each subfolder is one competitive-programming problem, the shapes a judge sets. It reads from standard input and answers on stdout, with its own input and output format, so nothing depends on an outside site being up. Nothing inside the programs is instrumented. The times below are wall time of the whole process, measured from outside, input and interpreter startup included. The C reference reads the same input with scanf.

	problem	plain	ppy build --unsafe	--standalone --unsafe	C (gcc)	C (clang)
15a	N-Queens	142.0 ms	48.7 ms	5.7 ms	4.8 ms	5.4 ms
15b	shortest path	4858.5 ms	4426.9 ms	111.7 ms	142.2 ms	134.3 ms
15c	substring search	292.4 ms	53.3 ms	—	9.3 ms	9.0 ms
15d	range sums	1571.7 ms	1524.4 ms	30.1 ms	50.8 ms	50.9 ms
15e	longest increasing subsequence	541.6 ms	111.5 ms	41.5 ms	57.7 ms	53.9 ms
15f	counting inversions	603.7 ms	98.4 ms	38.6 ms	43.5 ms	44.2 ms

Every cell is the mean of five runs, recorded in measurements.json with the machine it was measured on. bench.py reproduces it, and scripts/refresh.py says when a number here has drifted. Both fail if the paths stop agreeing on the answer. Bold is faster than both C references.

ppy run is left out because it compiles before it runs. That adds a flat two seconds or so to every row, and it is the development path rather than the one to submit.

ppy build --unsafe and --standalone --unsafe

• ppy build --unsafe is the hybrid. The kernels are native, but the glue around them (main, the buffers, print of a Python int) is the optimized Python the build wrote. So the binary embeds an interpreter and imports the runtime before the program begins. That is ~35 ms, and on the smaller problems it is most of what separates the column from C. Where main reads its input line by line (1.2 million edge lines in 15b, 400 thousand command lines in 15d), each ppy.input[...]() is a call into the runtime from the interpreter, about a microsecond apiece. Those reads are most of the plain and ppy build --unsafe columns there.
• --standalone is a binary with no CPython in it at all. ldd shows libc and nothing else, and ppy.scan[int]() lowers to the same buffered scan of standard input that scanf does. The whole N-Queens executable is 17.0 KB against the C one's 16.1 KB.

Four of the five standalone rows beat both C references for one reason: ppy.scan reads into memory faster than scanf parses. N-Queens reads a single integer and then computes, so there is nothing for a faster reader to win back, and it stays behind.

Reading input

n = ppy.input[int]()                 # one line, as int(input()) reads it
a, b = ppy.input[tuple[int, int]]()  # one line, exactly two fields
row = ppy.input[Buffer[int]]()       # one line of integers, into a buffer
values = ppy.scan[Buffer[int]](n)    # n integer tokens, whatever lines they are on

ppy.input[T]() reads one line the way the builtin input() does and types the result from T. ppy.scan[T]() reads tokens across lines. Both buffer reads go straight into memory rather than building a Python object per field: 12.6 ms for 500k integers against 49.8 ms for sys.stdin.read().split() and 20.6 ms for C's scanf.

The conversion writes the line reads for you, and each reads the line the original read:

• int(input()) becomes ppy.input[int]().
• a, b = map(int, input().split()) becomes the tuple read.
• array.array("q", map(int, input().split())) becomes the buffer line read.

A module that also touches sys.stdin keeps the input it has, since the typed reader owns the file descriptor.

Limitations

The standalone column covers five of the six problems. A standalone build needs everything main reaches to be native, and substring search reads its text as a token, for which ppy.read_token has no standalone lowering yet. The standalone variants trade try/except and array.array for the subset and live in standalone/. bench.py builds them from there and holds them to the same answer as every other column.

Where the code comes from

algorithms.ppy is hand-written; there is no .py source and no conversion step. algorithms.c is the same eight kernels hand-written in C, with the same workloads and the same printed answers.

Five of the six problem solutions are written as ordinary Python and converted: ppy convert <name>.py --promote-buffers writes the .ppy beside it, and examples/verify_conversions.py checks that the committed file is exactly that. 15c is hand-written, because the character buffer it wants has no plain-Python spelling that converts to it. The six problem subfolders say in their own READMEs which of them are generated.

Read on: Reading input · Native lowering · Performance · Buffers and JIT

Problems

• 15a_nqueens
• 15b_dijkstra
• 15c_kmp
• 15d_segment_tree
• 15e_lis
• 15f_input

15_algorithms/algorithms.ppy

import array
import time

import ppy
from ppy import Buffer


@ppy.opt(3)
def sieve_count(flags: Buffer[int], limit: int) -> int:
    for i in range(limit):
        flags[i] = 1
    flags[0] = 0
    if limit > 1:
        flags[1] = 0
    p: int = 2
    while p * p < limit:
        if flags[p] == 1:
            multiple: int = p * p
            while multiple < limit:
                flags[multiple] = 0
                multiple += p
        p += 1
    total: int = 0
    for i in range(limit):
        total += flags[i]
    return total


@ppy.pure
@ppy.opt(3)
def collatz_longest(limit: int) -> int:
    best: int = 0
    for start in range(1, limit):
        n: int = start
        steps: int = 0
        while n != 1:
            if n % 2 == 0:
                n = n // 2
            else:
                n = 3 * n + 1
            steps += 1
        best = max(best, steps)
    return best


@ppy.opt(3)
def knapsack(weights: Buffer[int], values: Buffer[int], table: Buffer[int], capacity: int) -> int:
    for i in range(capacity + 1):
        table[i] = 0
    for item in range(len(weights)):
        weight: int = weights[item]
        value: int = values[item]
        room: int = capacity
        while room >= weight:
            candidate: int = table[room - weight] + value
            if candidate > table[room]:
                table[room] = candidate
            room -= 1
    return table[capacity]


@ppy.opt(3)
def edit_distance(a: Buffer[int], b: Buffer[int], row: Buffer[int]) -> int:
    width: int = len(b)
    for j in range(width + 1):
        row[j] = j
    for i in range(len(a)):
        previous: int = row[0]
        row[0] = i + 1
        for j in range(width):
            current: int = row[j + 1]
            cost: int = 0
            if a[i] != b[j]:
                cost = 1
            best: int = previous + cost
            best = min(best, row[j] + 1)
            best = min(best, current + 1)
            row[j + 1] = best
            previous = current
    return row[width]


@ppy.opt(3)
def floyd_warshall(dist: Buffer[int], n: int) -> int:
    for k in range(n):
        for i in range(n):
            through: int = dist[i * n + k]
            if through < 1000000000:
                for j in range(n):
                    candidate: int = through + dist[k * n + j]
                    if candidate < dist[i * n + j]:
                        dist[i * n + j] = candidate
    total: int = 0
    for i in range(n * n):
        if dist[i] < 1000000000:
            total += dist[i]
    return total


@ppy.opt(3)
def matmul(a: Buffer[float], b: Buffer[float], out: Buffer[float], n: int) -> float:
    for i in range(n):
        for j in range(n):
            total: float = 0.0
            for k in range(n):
                total += a[i * n + k] * b[k * n + j]
            out[i * n + j] = total
    return out[n * n - 1]


@ppy.opt(3)
def union_find(parent: Buffer[int], edges: Buffer[int], n: int) -> int:
    for i in range(n):
        parent[i] = i
    pairs: int = len(edges) // 2
    for e in range(pairs):
        a: int = edges[e * 2]
        b: int = edges[e * 2 + 1]
        while parent[a] != a:
            parent[a] = parent[parent[a]]
            a = parent[a]
        while parent[b] != b:
            parent[b] = parent[parent[b]]
            b = parent[b]
        if a != b:
            parent[a] = b
    components: int = 0
    for i in range(n):
        if parent[i] == i:
            components += 1
    return components


@ppy.pure
@ppy.opt(3)
def modpow(base: int, exponent: int, modulus: int) -> int:
    result: int = 1
    b: int = base % modulus
    e: int = exponent
    while e > 0:
        if e % 2 == 1:
            result = result * b % modulus
        b = b * b % modulus
        e = e // 2
    return result


@ppy.pure
@ppy.opt(3)
def count_primes_fermat(limit: int) -> int:
    found: int = 0
    for n in range(3, limit, 2):
        if modpow(2, n - 1, n) == 1:
            found += 1
    return found


def run(label: str, seconds: float, answer: int) -> None:
    print(f"{label:<18s} {seconds * 1000.0:9.1f} ms   -> {answer}")


def main() -> None:
    flags = array.array("q", [0] * 2000000)
    start: float = time.perf_counter()
    primes: int = sieve_count(flags, 2000000)
    run("sieve 2e6", time.perf_counter() - start, primes)

    start = time.perf_counter()
    longest: int = collatz_longest(300000)
    run("collatz 3e5", time.perf_counter() - start, longest)

    weights = array.array("q", [(i * 7919) % 97 + 1 for i in range(400)])
    values = array.array("q", [(i * 104729) % 1000 + 1 for i in range(400)])
    table = array.array("q", [0] * 20001)
    start = time.perf_counter()
    best: int = knapsack(weights, values, table, 20000)
    run("knapsack 400x2e4", time.perf_counter() - start, best)

    a = array.array("q", [(i * 31) % 26 for i in range(2000)])
    b = array.array("q", [(i * 17) % 26 for i in range(2000)])
    row = array.array("q", [0] * 2001)
    start = time.perf_counter()
    distance: int = edit_distance(a, b, row)
    run("edit 2000x2000", time.perf_counter() - start, distance)

    size: int = 220
    dist = array.array("q", [0] * (size * size))
    for i in range(size):
        for j in range(size):
            if i == j:
                dist[i * size + j] = 0
            else:
                dist[i * size + j] = (i * 7 + j * 13) % 100 + 1
    start = time.perf_counter()
    total: int = floyd_warshall(dist, size)
    run("floyd 220", time.perf_counter() - start, total)

    n: int = 220
    left = array.array("d", [float((i * 31) % 17) for i in range(n * n)])
    right = array.array("d", [float((i * 13) % 23) for i in range(n * n)])
    out = array.array("d", [0.0] * (n * n))
    start = time.perf_counter()
    corner: float = matmul(left, right, out, n)
    run("matmul 220", time.perf_counter() - start, int(corner))

    nodes: int = 500000
    parent = array.array("q", [0] * nodes)
    edges = array.array("q", [(i * 7919) % nodes for i in range(2000000)])
    start = time.perf_counter()
    components: int = union_find(parent, edges, nodes)
    run("union-find 5e5", time.perf_counter() - start, components)

    start = time.perf_counter()
    fermat: int = count_primes_fermat(60000)
    run("fermat 6e4", time.perf_counter() - start, fermat)


main()

Counterpart programs

The programs the comparison above measured, each written the way its tool expects. The PPy one is first.

algorithms.mojo (Mojo)

"""The eight kernels of `algorithms.ppy` in Mojo: the same loops over Lists, typed by hand."""

from std.time import perf_counter_ns


def sieve_count(mut flags: List[Int64], limit: Int) -> Int64:
    for i in range(limit):
        flags[i] = 1
    flags[0] = 0
    if limit > 1:
        flags[1] = 0
    var p = 2
    while p * p < limit:
        if flags[p] == 1:
            var multiple = p * p
            while multiple < limit:
                flags[multiple] = 0
                multiple += p
        p += 1
    var total: Int64 = 0
    for i in range(limit):
        total += flags[i]
    return total


def collatz_longest(limit: Int) -> Int64:
    var best: Int64 = 0
    for start in range(1, limit):
        var n = Int64(start)
        var steps: Int64 = 0
        while n != 1:
            if n % 2 == 0:
                n = n // 2
            else:
                n = 3 * n + 1
            steps += 1
        best = max(best, steps)
    return best


def knapsack(weights: List[Int64], values: List[Int64], mut table: List[Int64], capacity: Int) -> Int64:
    for i in range(capacity + 1):
        table[i] = 0
    for item in range(len(weights)):
        var weight = Int(weights[item])
        var value = values[item]
        var room = capacity
        while room >= weight:
            var candidate = table[room - weight] + value
            if candidate > table[room]:
                table[room] = candidate
            room -= 1
    return table[capacity]


def edit_distance(a: List[Int64], b: List[Int64], mut row: List[Int64]) -> Int64:
    var width = len(b)
    for j in range(width + 1):
        row[j] = Int64(j)
    for i in range(len(a)):
        var previous = row[0]
        row[0] = Int64(i + 1)
        for j in range(width):
            var current = row[j + 1]
            var cost: Int64 = 0
            if a[i] != b[j]:
                cost = 1
            var best = previous + cost
            best = min(best, row[j] + 1)
            best = min(best, current + 1)
            row[j + 1] = best
            previous = current
    return row[width]


def floyd_warshall(mut dist: List[Int64], n: Int) -> Int64:
    for k in range(n):
        for i in range(n):
            var through = dist[i * n + k]
            if through < 1000000000:
                for j in range(n):
                    var candidate = through + dist[k * n + j]
                    if candidate < dist[i * n + j]:
                        dist[i * n + j] = candidate
    var total: Int64 = 0
    for i in range(n * n):
        if dist[i] < 1000000000:
            total += dist[i]
    return total


def matmul(a: List[Float64], b: List[Float64], mut out: List[Float64], n: Int) -> Float64:
    for i in range(n):
        for j in range(n):
            var total: Float64 = 0.0
            for k in range(n):
                total += a[i * n + k] * b[k * n + j]
            out[i * n + j] = total
    return out[n * n - 1]


def union_find(mut parent: List[Int64], edges: List[Int64], n: Int) -> Int64:
    for i in range(n):
        parent[i] = Int64(i)
    var pairs = len(edges) // 2
    for e in range(pairs):
        var a = Int(edges[e * 2])
        var b = Int(edges[e * 2 + 1])
        while Int(parent[a]) != a:
            parent[a] = parent[Int(parent[a])]
            a = Int(parent[a])
        while Int(parent[b]) != b:
            parent[b] = parent[Int(parent[b])]
            b = Int(parent[b])
        if a != b:
            parent[a] = Int64(b)
    var components: Int64 = 0
    for i in range(n):
        if Int(parent[i]) == i:
            components += 1
    return components


def modpow(base: Int64, exponent: Int64, modulus: Int64) -> Int64:
    var result: Int64 = 1
    var b = base % modulus
    var e = exponent
    while e > 0:
        if e % 2 == 1:
            result = result * b % modulus
        b = b * b % modulus
        e = e // 2
    return result


def count_primes_fermat(limit: Int) -> Int64:
    var found: Int64 = 0
    for n in range(3, limit, 2):
        if modpow(2, Int64(n - 1), Int64(n)) == 1:
            found += 1
    return found


def report(label: String, started: Int, answer: Int64):
    var ms = Float64(perf_counter_ns() - started) / 1e6
    print(label, ms, "ms   ->", answer)


def main():
    var flags = List[Int64](length=2000000, fill=0)
    var started = perf_counter_ns()
    var primes = sieve_count(flags, 2000000)
    report("sieve 2e6", started, primes)

    started = perf_counter_ns()
    var longest = collatz_longest(300000)
    report("collatz 3e5", started, longest)

    var weights = List[Int64](capacity=400)
    var values = List[Int64](capacity=400)
    for i in range(400):
        weights.append(Int64((i * 7919) % 97 + 1))
        values.append(Int64((i * 104729) % 1000 + 1))
    var table = List[Int64](length=20001, fill=0)
    started = perf_counter_ns()
    var best = knapsack(weights, values, table, 20000)
    report("knapsack 400x2e4", started, best)

    var a = List[Int64](capacity=2000)
    var b = List[Int64](capacity=2000)
    for i in range(2000):
        a.append(Int64((i * 31) % 26))
        b.append(Int64((i * 17) % 26))
    var row = List[Int64](length=2001, fill=0)
    started = perf_counter_ns()
    var distance = edit_distance(a, b, row)
    report("edit 2000x2000", started, distance)

    var size = 220
    var dist = List[Int64](length=size * size, fill=0)
    for i in range(size):
        for j in range(size):
            if i == j:
                dist[i * size + j] = 0
            else:
                dist[i * size + j] = Int64((i * 7 + j * 13) % 100 + 1)
    started = perf_counter_ns()
    var total = floyd_warshall(dist, size)
    report("floyd 220", started, total)

    var n = 220
    var left = List[Float64](capacity=n * n)
    var right = List[Float64](capacity=n * n)
    for i in range(n * n):
        left.append(Float64((i * 31) % 17))
        right.append(Float64((i * 13) % 23))
    var out = List[Float64](length=n * n, fill=0.0)
    started = perf_counter_ns()
    var corner = matmul(left, right, out, n)
    report("matmul 220", started, Int64(corner))

    var nodes = 500000
    var parent = List[Int64](length=nodes, fill=0)
    var edges = List[Int64](capacity=2000000)
    for i in range(2000000):
        edges.append(Int64((i * 7919) % nodes))
    started = perf_counter_ns()
    var components = union_find(parent, edges, nodes)
    report("union-find 5e5", started, components)

    started = perf_counter_ns()
    var fermat = count_primes_fermat(60000)
    report("fermat 6e4", started, fermat)

algorithms_codon.py (Python)

"""The eight kernels of `algorithms.ppy` under Codon: the same loops over typed lists."""

from time import time


def sieve_count(flags: List[int], limit: int) -> int:
    for i in range(limit):
        flags[i] = 1
    flags[0] = 0
    if limit > 1:
        flags[1] = 0
    p = 2
    while p * p < limit:
        if flags[p] == 1:
            multiple = p * p
            while multiple < limit:
                flags[multiple] = 0
                multiple += p
        p += 1
    total = 0
    for i in range(limit):
        total += flags[i]
    return total


def collatz_longest(limit: int) -> int:
    best = 0
    for start in range(1, limit):
        n = start
        steps = 0
        while n != 1:
            if n % 2 == 0:
                n = n // 2
            else:
                n = 3 * n + 1
            steps += 1
        best = max(best, steps)
    return best


def knapsack(weights: List[int], values: List[int], table: List[int], capacity: int) -> int:
    for i in range(capacity + 1):
        table[i] = 0
    for item in range(len(weights)):
        weight = weights[item]
        value = values[item]
        room = capacity
        while room >= weight:
            candidate = table[room - weight] + value
            if candidate > table[room]:
                table[room] = candidate
            room -= 1
    return table[capacity]


def edit_distance(a: List[int], b: List[int], row: List[int]) -> int:
    width = len(b)
    for j in range(width + 1):
        row[j] = j
    for i in range(len(a)):
        previous = row[0]
        row[0] = i + 1
        for j in range(width):
            current = row[j + 1]
            cost = 0
            if a[i] != b[j]:
                cost = 1
            best = previous + cost
            best = min(best, row[j] + 1)
            best = min(best, current + 1)
            row[j + 1] = best
            previous = current
    return row[width]


def floyd_warshall(dist: List[int], n: int) -> int:
    for k in range(n):
        for i in range(n):
            through = dist[i * n + k]
            if through < 1000000000:
                for j in range(n):
                    candidate = through + dist[k * n + j]
                    if candidate < dist[i * n + j]:
                        dist[i * n + j] = candidate
    total = 0
    for i in range(n * n):
        if dist[i] < 1000000000:
            total += dist[i]
    return total


def matmul(a: List[float], b: List[float], out: List[float], n: int) -> float:
    for i in range(n):
        for j in range(n):
            total = 0.0
            for k in range(n):
                total += a[i * n + k] * b[k * n + j]
            out[i * n + j] = total
    return out[n * n - 1]


def union_find(parent: List[int], edges: List[int], n: int) -> int:
    for i in range(n):
        parent[i] = i
    pairs = len(edges) // 2
    for e in range(pairs):
        a = edges[e * 2]
        b = edges[e * 2 + 1]
        while parent[a] != a:
            parent[a] = parent[parent[a]]
            a = parent[a]
        while parent[b] != b:
            parent[b] = parent[parent[b]]
            b = parent[b]
        if a != b:
            parent[a] = b
    components = 0
    for i in range(n):
        if parent[i] == i:
            components += 1
    return components


def modpow(base: int, exponent: int, modulus: int) -> int:
    result = 1
    b = base % modulus
    e = exponent
    while e > 0:
        if e % 2 == 1:
            result = result * b % modulus
        b = b * b % modulus
        e = e // 2
    return result


def count_primes_fermat(limit: int) -> int:
    found = 0
    for n in range(3, limit, 2):
        if modpow(2, n - 1, n) == 1:
            found += 1
    return found


def run(label: str, seconds: float, answer: int):
    print(label, round(seconds * 1000.0, 1), "ms   ->", answer)


def main():
    flags = [0 for _ in range(2000000)]
    start = time()
    primes = sieve_count(flags, 2000000)
    run("sieve 2e6", time() - start, primes)

    start = time()
    longest = collatz_longest(300000)
    run("collatz 3e5", time() - start, longest)

    weights = [(i * 7919) % 97 + 1 for i in range(400)]
    values = [(i * 104729) % 1000 + 1 for i in range(400)]
    table = [0 for _ in range(20001)]
    start = time()
    best = knapsack(weights, values, table, 20000)
    run("knapsack 400x2e4", time() - start, best)

    a = [(i * 31) % 26 for i in range(2000)]
    b = [(i * 17) % 26 for i in range(2000)]
    row = [0 for _ in range(2001)]
    start = time()
    distance = edit_distance(a, b, row)
    run("edit 2000x2000", time() - start, distance)

    size = 220
    dist = [0 for _ in range(size * size)]
    for i in range(size):
        for j in range(size):
            dist[i * size + j] = 0 if i == j else (i * 7 + j * 13) % 100 + 1
    start = time()
    total = floyd_warshall(dist, size)
    run("floyd 220", time() - start, total)

    n = 220
    left = [float((i * 31) % 17) for i in range(n * n)]
    right = [float((i * 13) % 23) for i in range(n * n)]
    out = [0.0 for _ in range(n * n)]
    start = time()
    corner = matmul(left, right, out, n)
    run("matmul 220", time() - start, int(corner))

    nodes = 500000
    parent = [0 for _ in range(nodes)]
    edges = [(i * 7919) % nodes for i in range(2000000)]
    start = time()
    components = union_find(parent, edges, nodes)
    run("union-find 5e5", time() - start, components)

    start = time()
    fermat = count_primes_fermat(60000)
    run("fermat 6e4", time() - start, fermat)


main()

algorithms_numba.py (Python)

"""The eight kernels of `algorithms.ppy` under Numba's `@njit`: the same loops over NumPy arrays."""

import time

import numpy as np
from numba import njit


@njit
def sieve_count(flags, limit):
    for i in range(limit):
        flags[i] = 1
    flags[0] = 0
    if limit > 1:
        flags[1] = 0
    p = 2
    while p * p < limit:
        if flags[p] == 1:
            multiple = p * p
            while multiple < limit:
                flags[multiple] = 0
                multiple += p
        p += 1
    total = 0
    for i in range(limit):
        total += flags[i]
    return total


@njit
def collatz_longest(limit):
    best = 0
    for start in range(1, limit):
        n = start
        steps = 0
        while n != 1:
            if n % 2 == 0:
                n = n // 2
            else:
                n = 3 * n + 1
            steps += 1
        best = max(best, steps)
    return best


@njit
def knapsack(weights, values, table, capacity):
    for i in range(capacity + 1):
        table[i] = 0
    for item in range(len(weights)):
        weight = weights[item]
        value = values[item]
        room = capacity
        while room >= weight:
            candidate = table[room - weight] + value
            if candidate > table[room]:
                table[room] = candidate
            room -= 1
    return table[capacity]


@njit
def edit_distance(a, b, row):
    width = len(b)
    for j in range(width + 1):
        row[j] = j
    for i in range(len(a)):
        previous = row[0]
        row[0] = i + 1
        for j in range(width):
            current = row[j + 1]
            cost = 0
            if a[i] != b[j]:
                cost = 1
            best = previous + cost
            best = min(best, row[j] + 1)
            best = min(best, current + 1)
            row[j + 1] = best
            previous = current
    return row[width]


@njit
def floyd_warshall(dist, n):
    for k in range(n):
        for i in range(n):
            through = dist[i * n + k]
            if through < 1000000000:
                for j in range(n):
                    candidate = through + dist[k * n + j]
                    if candidate < dist[i * n + j]:
                        dist[i * n + j] = candidate
    total = 0
    for i in range(n * n):
        if dist[i] < 1000000000:
            total += dist[i]
    return total


@njit
def matmul(a, b, out, n):
    for i in range(n):
        for j in range(n):
            total = 0.0
            for k in range(n):
                total += a[i * n + k] * b[k * n + j]
            out[i * n + j] = total
    return out[n * n - 1]


@njit
def union_find(parent, edges, n):
    for i in range(n):
        parent[i] = i
    pairs = len(edges) // 2
    for e in range(pairs):
        a = edges[e * 2]
        b = edges[e * 2 + 1]
        while parent[a] != a:
            parent[a] = parent[parent[a]]
            a = parent[a]
        while parent[b] != b:
            parent[b] = parent[parent[b]]
            b = parent[b]
        if a != b:
            parent[a] = b
    components = 0
    for i in range(n):
        if parent[i] == i:
            components += 1
    return components


@njit
def modpow(base, exponent, modulus):
    result = 1
    b = base % modulus
    e = exponent
    while e > 0:
        if e % 2 == 1:
            result = result * b % modulus
        b = b * b % modulus
        e = e // 2
    return result


@njit
def count_primes_fermat(limit):
    found = 0
    for n in range(3, limit, 2):
        if modpow(2, n - 1, n) == 1:
            found += 1
    return found


def run(label, seconds, answer):
    print(f"{label:<18s} {seconds * 1000.0:9.1f} ms   -> {answer}")


def timed(label, call):
    call()  # compile
    start = time.perf_counter()
    answer = call()
    run(label, time.perf_counter() - start, answer)


def main():
    flags = np.zeros(2000000, dtype=np.int64)
    timed("sieve 2e6", lambda: sieve_count(flags, 2000000))
    timed("collatz 3e5", lambda: collatz_longest(300000))
    weights = np.array([(i * 7919) % 97 + 1 for i in range(400)], dtype=np.int64)
    values = np.array([(i * 104729) % 1000 + 1 for i in range(400)], dtype=np.int64)
    table = np.zeros(20001, dtype=np.int64)
    timed("knapsack 400x2e4", lambda: knapsack(weights, values, table, 20000))
    a = np.array([(i * 31) % 26 for i in range(2000)], dtype=np.int64)
    b = np.array([(i * 17) % 26 for i in range(2000)], dtype=np.int64)
    row = np.zeros(2001, dtype=np.int64)
    timed("edit 2000x2000", lambda: edit_distance(a, b, row))
    size = 220
    dist = np.zeros(size * size, dtype=np.int64)
    for i in range(size):
        for j in range(size):
            dist[i * size + j] = 0 if i == j else (i * 7 + j * 13) % 100 + 1
    fresh = dist.copy()
    floyd_warshall(fresh, size)  # compile on a copy: the kernel relaxes in place
    start = time.perf_counter()
    total = floyd_warshall(dist, size)
    run("floyd 220", time.perf_counter() - start, total)
    n = 220
    left = np.array([float((i * 31) % 17) for i in range(n * n)])
    right = np.array([float((i * 13) % 23) for i in range(n * n)])
    out = np.zeros(n * n)
    timed("matmul 220", lambda: int(matmul(left, right, out, n)))
    nodes = 500000
    parent = np.zeros(nodes, dtype=np.int64)
    edges = np.array([(i * 7919) % nodes for i in range(2000000)], dtype=np.int64)
    timed("union-find 5e5", lambda: union_find(parent, edges, nodes))
    timed("fermat 6e4", lambda: count_primes_fermat(60000))


main()

Source: examples/15_algorithms.