Buffers and JIT

Borrowed memory, reassociation, and specialization each set how fast a numeric function runs, and each is chosen once in a signature or a decorator.

Run it

python  buffers_and_jit.ppy
ppy     buffers_and_jit.ppy
ppy run buffers_and_jit.ppy

Buffer[T] is borrowed, a list is copied

@ppy.pure
@ppy.opt(3)
def total_view(values: Buffer[float]) -> float:
    result: float = 0.0
    for i in range(len(values)):
        result += values[i]
    return result

total_list takes a list[float] and is copied into a buffer on every call. total_view takes a Buffer[float] (a memoryview over an array.array) and reads the caller's memory in place. Same loop, same sum, about 10× apart on 8,192 elements, because one of them pays for 8,192 boxed floats per call.

A buffer is borrowed unless the signature says otherwise. The callee may not keep it (E1612) or write through it without Mut[...] (E1613).

Reassociation is opt-in

dot_strict and dot_relaxed are the same dot product. The strict one keeps CPython's accumulation order and matches it bit for bit. The @ppy.fastmath one lets LLVM vectorize the reduction and differs in the last bits; the program prints the gap. You make the choice per function, and the default is the exact one.

@ppy.jit specializes on the values it sees

@ppy.jit(threshold=4, max_specializations=4)
@ppy.pure
@ppy.opt(3)
def digest_jit(values: Buffer[int], modulus: int) -> int:

After four calls with the same modulus, a version specialized to that constant is compiled, and % 1000003 becomes multiply-and-shift instead of a division. The guard on the value is in the generated C wrapper. A call with a different modulus falls back to the generic version or compiles another specialization, up to four.

Compared with Numba, Cython, NumPy, and C

The three kernels over 8,192 elements: a sum, a dot product in order and reassociated, and a modular digest. Each call is timed as the mean of 2,000 calls, best of five rounds, milliseconds per call over five processes. The programs are in compare/: kernels_bench.ppy, kernels_numba.py, kernels_cython.pyx, kernels_numpy.py, kernels.c. At this size a call is a few microseconds, so the table is as much about the call boundary as about the loop.

PPy is the loop as written, with a Buffer[float] borrowed from an array.array and @ppy.fastmath where reassociation is allowed:

@ppy.pure
@ppy.opt(3)
def dot(a: Buffer[float], b: Buffer[float]) -> float:
    result: float = 0.0
    for i in range(len(a)):
        result += a[i] * b[i]
    return result

Numba is the same loop under @njit, with no annotations, NumPy arrays, and fastmath=True for the relaxed dot. Its dispatcher types the arguments on every call:

@njit
def dot(a, b):
    result = 0.0
    for i in range(len(a)):
        result += a[i] * b[i]
    return result

Cython is typed memoryviews and cdef locals in a .pyx, built into an extension by cythonize -i. boundscheck=False is what makes it a C loop, and there is no relaxed dot to write:

cpdef double dot(double[::1] a, double[::1] b):
    cdef double result = 0.0
    cdef Py_ssize_t i
    for i in range(a.shape[0]):
        result += a[i] * b[i]
    return result

NumPy has no loop: np.dot is BLAS, np.sum is pairwise, and (counts % m).sum() builds a temporary.

C is the loop with no Python around it, timed inside the process: the floor.

	PPy ppy run	Numba @njit	Cython	NumPy	C (the loop alone)
total	0.0035 ± 0.0000	0.0036 ± 0.0000	0.0036 ± 0.0001	0.0021 ± 0.0001	0.0034 ± 0.0001
dot	0.0035 ± 0.0000	0.0036 ± 0.0000	0.0036 ± 0.0000	0.0012 ± 0.0001	0.0035 ± 0.0001
dot_relaxed	0.0011 ± 0.0001	0.0011 ± 0.0000	0.0036 ± 0.0001	0.0012 ± 0.0001	0.0034 ± 0.0000
digest	0.0072 ± 0.0000	0.0112 ± 0.0001	0.0078 ± 0.0001	0.0199 ± 0.0002	0.0045 ± 0.0001

The ordered sums sit on the C loop in PPy, Numba, and Cython alike, a microsecond of call above it. The reassociated dot is where @ppy.fastmath and Numba's fastmath=True let the vectorizer in, and Cython, with no such switch, stays scalar.

NumPy's dot is BLAS and the fastest row. Its digest pays for an 8,192-element temporary and is the slowest. The C digest is the loop's floor; what stands between it and the others is the price of a call through the interpreter.

Intel Core Ultra 9 386H; Numba 0.67.0, Cython 3.3.0, NumPy 2.5.3 on CPython 3.12.13, gcc 13.3, PPy on CPython 3.13.13, from a checkout on a native filesystem.

What it prints

python buffers_and_jit.ppy

list[float] copied        104.387 us   
Buffer[float] borrowed    158.837 us   0.66x, same sum: True
dot, strict order         277.649 us   
dot, @ppy.fastmath        257.140 us   1.08x, differs by 0.00e+00
digest, generic           263.962 us   
digest, @ppy.jit          260.730 us   1.01x, same: True

ppy buffers_and_jit.ppy

list[float] copied        102.127 us   
Buffer[float] borrowed    162.701 us   0.63x, same sum: True
dot, strict order         261.217 us   
dot, @ppy.fastmath        259.543 us   1.01x, differs by 0.00e+00
digest, generic           270.856 us   
digest, @ppy.jit          268.401 us   1.01x, same: True

ppy run buffers_and_jit.ppy

list[float] copied         11.790 us   
Buffer[float] borrowed      3.726 us   3.16x, same sum: True
dot, strict order           4.196 us   
dot, @ppy.fastmath          0.849 us   4.94x, differs by 4.07e-10
digest, generic            11.656 us   
digest, @ppy.jit            8.874 us   1.31x, same: True

Read on: Directives and markers · Algorithms

buffers_and_jit.ppy is hand-written; there is no .py source and no conversion step.

12_buffers_and_jit/buffers_and_jit.ppy

import array
import time

import ppy
from ppy import Buffer


@ppy.pure
@ppy.opt(3)
def total_list(values: list[float]) -> float:
    result: float = 0.0
    for i in range(len(values)):
        result += values[i]
    return result


@ppy.pure
@ppy.opt(3)
def total_view(values: Buffer[float]) -> float:
    result: float = 0.0
    for i in range(len(values)):
        result += values[i]
    return result


@ppy.pure
@ppy.opt(3)
def dot_strict(a: Buffer[float], b: Buffer[float]) -> float:
    result: float = 0.0
    for i in range(len(a)):
        result += a[i] * b[i]
    return result


@ppy.pure
@ppy.fastmath
@ppy.opt(3)
def dot_relaxed(a: Buffer[float], b: Buffer[float]) -> float:
    result: float = 0.0
    for i in range(len(a)):
        result += a[i] * b[i]
    return result


@ppy.jit(threshold=4, max_specializations=4)
@ppy.pure
@ppy.opt(3)
def digest_jit(values: Buffer[int], modulus: int) -> int:
    result: int = 0
    for i in range(len(values)):
        result += values[i] % modulus
    return result


@ppy.pure
@ppy.opt(3)
def digest_generic(values: Buffer[int], modulus: int) -> int:
    result: int = 0
    for i in range(len(values)):
        result += values[i] % modulus
    return result


def report(label: str, us: float, note: str) -> None:
    print(f"{label:24s} {us:8.3f} us   {note}")


def main() -> None:
    size: int = 8192
    raw: list[float] = [float(i) * 0.001 for i in range(size)]
    x = array.array("d", raw)
    y = array.array("d", [float(i) * 0.002 for i in range(size)])
    counts = array.array("q", [i * 7919 for i in range(size)])

    rounds: int = 2000
    for _i in range(400):
        total_list(raw)
        total_view(x)
        dot_strict(x, y)
        dot_relaxed(x, y)
        digest_jit(counts, 1000003)
        digest_generic(counts, 1000003)

    start: float = time.perf_counter()
    for _i in range(rounds):
        copied = total_list(raw)
    copied_us: float = (time.perf_counter() - start) / rounds * 1e6

    start = time.perf_counter()
    for _i in range(rounds):
        borrowed = total_view(x)
    borrowed_us: float = (time.perf_counter() - start) / rounds * 1e6

    report("list[float] copied", copied_us, "")
    ratio = copied_us / borrowed_us
    report("Buffer[float] borrowed", borrowed_us, f"{ratio:.2f}x, same sum: {copied == borrowed}")

    start = time.perf_counter()
    for _i in range(rounds):
        strict = dot_strict(x, y)
    strict_us: float = (time.perf_counter() - start) / rounds * 1e6

    start = time.perf_counter()
    for _i in range(rounds):
        relaxed = dot_relaxed(x, y)
    relaxed_us: float = (time.perf_counter() - start) / rounds * 1e6

    report("dot, strict order", strict_us, "")
    ratio = strict_us / relaxed_us
    gap = abs(strict - relaxed)
    report("dot, @ppy.fastmath", relaxed_us, f"{ratio:.2f}x, differs by {gap:.2e}")

    start = time.perf_counter()
    for _i in range(rounds):
        generic = digest_generic(counts, 1000003)
    generic_us: float = (time.perf_counter() - start) / rounds * 1e6

    start = time.perf_counter()
    for _i in range(rounds):
        specialized = digest_jit(counts, 1000003)
    specialized_us: float = (time.perf_counter() - start) / rounds * 1e6

    report("digest, generic", generic_us, "")
    ratio = generic_us / specialized_us
    report("digest, @ppy.jit", specialized_us, f"{ratio:.2f}x, same: {generic == specialized}")


if __name__ == "__main__":
    main()

Counterpart programs

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

kernels_bench.ppy (PPy)

"""The kernels of `buffers_and_jit.ppy` on 8,192 elements, microseconds per call: the PPY side."""

import array
import time

import ppy
from ppy import Buffer


@ppy.pure
@ppy.opt(3)
def total(values: Buffer[float]) -> float:
    result: float = 0.0
    for i in range(len(values)):
        result += values[i]
    return result


@ppy.pure
@ppy.opt(3)
def dot(a: Buffer[float], b: Buffer[float]) -> float:
    result: float = 0.0
    for i in range(len(a)):
        result += a[i] * b[i]
    return result


@ppy.pure
@ppy.fastmath
@ppy.opt(3)
def dot_relaxed(a: Buffer[float], b: Buffer[float]) -> float:
    result: float = 0.0
    for i in range(len(a)):
        result += a[i] * b[i]
    return result


@ppy.jit(threshold=4, max_specializations=4)
@ppy.pure
@ppy.opt(3)
def digest(values: Buffer[int], modulus: int) -> int:
    result: int = 0
    for i in range(len(values)):
        result += values[i] % modulus
    return result


def main() -> None:
    size = 8192
    x = array.array("d", [float(i) * 0.001 for i in range(size)])
    y = array.array("d", [float(i) * 0.002 for i in range(size)])
    counts = array.array("q", [i * 7919 for i in range(size)])
    for _ in range(50):
        total(x)
        dot(x, y)
        dot_relaxed(x, y)
        digest(counts, 1000003)
    rounds = 2000
    best = [1e9, 1e9, 1e9, 1e9]
    for _ in range(5):
        started = time.perf_counter()
        for _i in range(rounds):
            total(x)
        best[0] = min(best[0], (time.perf_counter() - started) * 1000.0 / rounds)
        started = time.perf_counter()
        for _i in range(rounds):
            dot(x, y)
        best[1] = min(best[1], (time.perf_counter() - started) * 1000.0 / rounds)
        started = time.perf_counter()
        for _i in range(rounds):
            dot_relaxed(x, y)
        best[2] = min(best[2], (time.perf_counter() - started) * 1000.0 / rounds)
        started = time.perf_counter()
        for _i in range(rounds):
            digest(counts, 1000003)
        best[3] = min(best[3], (time.perf_counter() - started) * 1000.0 / rounds)
    labels = ["total", "dot", "dot_relaxed", "digest"]
    for label, took in zip(labels, best, strict=True):
        print(f"# {label}: {took:.4f} ms")
    relaxed = round(dot_relaxed(x, y), 3)
    print(round(total(x), 6), round(dot(x, y), 6), relaxed, digest(counts, 1000003))


main()

kernels.c (C)

/* The same kernels in C, the loops alone: what the machine does with no call boundary. */
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

static double total(const double *values, int64_t n) {
    double result = 0.0;
    for (int64_t i = 0; i < n; i++) result += values[i];
    return result;
}

static double dot(const double *a, const double *b, int64_t n) {
    double result = 0.0;
    for (int64_t i = 0; i < n; i++) result += a[i] * b[i];
    return result;
}

static int64_t digest(const int64_t *values, int64_t n, int64_t modulus) {
    int64_t result = 0;
    for (int64_t i = 0; i < n; i++) result += values[i] % modulus;
    return result;
}

/* `round(v, places)` printed the way Python prints it: no trailing zeros. */
static const char *shortest(char *out, double value, int places) {
    snprintf(out, 64, "%.*f", places, value);
    char *end = out + strlen(out) - 1;
    while (*end == '0') *end-- = 0;
    if (*end == '.') *end = 0;
    return out;
}

static double now(void) {
    struct timespec t;
    clock_gettime(CLOCK_MONOTONIC, &t);
    return (double)t.tv_sec + (double)t.tv_nsec * 1e-9;
}

#define TIMED(label, expr, rounds, sink)                                           \
    do {                                                                           \
        double best = 1e9;                                                         \
        for (int r = 0; r < 5; r++) {                                              \
            double started = now();                                                \
            for (int i = 0; i < (rounds); i++) { sink = (expr); asm volatile("" : : "g"(&sink) : "memory"); } \
            double took = (now() - started) / (rounds);                            \
            if (took < best) best = took;                                          \
        }                                                                          \
        printf("# %s: %.4f ms\n", label, best * 1000.0);                           \
    } while (0)

int main(void) {
    const int64_t n = 8192;
    double *x = malloc(n * sizeof *x), *y = malloc(n * sizeof *y);
    int64_t *counts = malloc(n * sizeof *counts);
    for (int64_t i = 0; i < n; i++) {
        x[i] = (double)i * 0.001;
        y[i] = (double)i * 0.002;
        counts[i] = i * 7919;
    }
    volatile double d;
    volatile int64_t k;
    TIMED("total", total(x, n), 2000, d);
    TIMED("dot", dot(x, y, n), 2000, d);
    TIMED("dot_relaxed", dot(x, y, n), 2000, d);
    TIMED("digest", digest(counts, n, 1000003), 2000, k);
    char a[64], b[64], c[64];
    printf("%s %s %s %lld\n", shortest(a, total(x, n), 6), shortest(b, dot(x, y, n), 6),
           shortest(c, dot(x, y, n), 3), (long long)digest(counts, n, 1000003));
    free(x); free(y); free(counts);
    return 0;
}

kernels_cython.pyx (Cython)

# cython: language_level=3, boundscheck=False, wraparound=False
"""The same kernels in Cython: typed memoryviews and C locals, built by `cythonize`."""

import time


cpdef double total(double[::1] values):
    cdef double result = 0.0
    cdef Py_ssize_t i
    for i in range(values.shape[0]):
        result += values[i]
    return result


cpdef double dot(double[::1] a, double[::1] b):
    cdef double result = 0.0
    cdef Py_ssize_t i
    for i in range(a.shape[0]):
        result += a[i] * b[i]
    return result


cpdef long long digest(long long[::1] values, long long modulus):
    cdef long long result = 0
    cdef Py_ssize_t i
    for i in range(values.shape[0]):
        result += values[i] % modulus
    return result


def timed(label, run, rounds):
    best = 1e9
    answer = None
    for _ in range(5):
        started = time.perf_counter()
        for _i in range(rounds):
            answer = run()
        best = min(best, (time.perf_counter() - started) / rounds)
    print(f"# {label}: {best * 1000:.4f} ms")
    return answer


def main():
    import array
    size = 8192
    x = array.array("d", [float(i) * 0.001 for i in range(size)])
    y = array.array("d", [float(i) * 0.002 for i in range(size)])
    counts = array.array("q", [i * 7919 for i in range(size)])
    timed("total", lambda: total(x), 2000)
    timed("dot", lambda: dot(x, y), 2000)
    timed("dot_relaxed", lambda: dot(x, y), 2000)
    timed("digest", lambda: digest(counts, 1000003), 2000)
    print(round(total(x), 6), round(dot(x, y), 6), round(dot(x, y), 3), digest(counts, 1000003))


if __name__ == "__main__":
    main()

kernels_numba.py (Python)

"""The same kernels under Numba: `@njit` over NumPy arrays."""

import time

import numpy as np
from numba import njit


@njit(cache=False)
def total(values):
    result = 0.0
    for i in range(len(values)):
        result += values[i]
    return result


@njit(cache=False)
def dot(a, b):
    result = 0.0
    for i in range(len(a)):
        result += a[i] * b[i]
    return result


@njit(cache=False, fastmath=True)
def dot_relaxed(a, b):
    result = 0.0
    for i in range(len(a)):
        result += a[i] * b[i]
    return result


@njit(cache=False)
def digest(values, modulus):
    result = 0
    for i in range(len(values)):
        result += values[i] % modulus
    return result


def timed(label, run, rounds):
    best = 1e9
    answer = None
    for _ in range(5):
        started = time.perf_counter()
        for _i in range(rounds):
            answer = run()
        best = min(best, (time.perf_counter() - started) / rounds)
    print(f"# {label}: {best * 1000:.4f} ms")
    return answer


def main():
    size = 8192
    x = np.arange(size, dtype=np.float64) * 0.001
    y = np.arange(size, dtype=np.float64) * 0.002
    counts = np.arange(size, dtype=np.int64) * 7919
    for _ in range(50):
        total(x)
        dot(x, y)
        dot_relaxed(x, y)
        digest(counts, 1000003)
    timed("total", lambda: total(x), 2000)
    timed("dot", lambda: dot(x, y), 2000)
    timed("dot_relaxed", lambda: dot_relaxed(x, y), 2000)
    timed("digest", lambda: digest(counts, 1000003), 2000)
    print(round(total(x), 6), round(dot(x, y), 6), round(dot_relaxed(x, y), 3), digest(counts, 1000003))

main()

kernels_numpy.py (Python)

"""The same kernels as NumPy array expressions: no loop to write."""

import time

import numpy as np


def timed(label, run, rounds):
    best = 1e9
    answer = None
    for _ in range(5):
        started = time.perf_counter()
        for _i in range(rounds):
            answer = run()
        best = min(best, (time.perf_counter() - started) / rounds)
    print(f"# {label}: {best * 1000:.4f} ms")
    return answer


def main():
    size = 8192
    x = np.arange(size, dtype=np.float64) * 0.001
    y = np.arange(size, dtype=np.float64) * 0.002
    counts = np.arange(size, dtype=np.int64) * 7919
    timed("total", lambda: float(np.sum(x)), 2000)
    timed("dot", lambda: float(np.dot(x, y)), 2000)
    timed("dot_relaxed", lambda: float(np.dot(x, y)), 2000)
    timed("digest", lambda: int((counts % 1000003).sum()), 2000)
    print(
        round(float(np.sum(x)), 6),
        round(float(np.dot(x, y)), 6),
        round(float(np.dot(x, y)), 3),
        int((counts % 1000003).sum()),
    )

main()

Source: examples/12_buffers_and_jit.