the fast path out of python
Sum ten million floating-point numbers, two ways. A plain Python loop, and numpy:
python for-loop total += x : 289.9 ms
numpy arr.sum() : 3.4 msThe same ten million input values, and numpy is eighty-five times faster. I first called them the same additions and the answer bit-for-bit identical. Both statements survive only with qualifications that became important later. Numpy changes the grouping of the additions, and this particular input happens to sum exactly in binary floating point. The easy explanation was still tempting: the Python loop spent nearly all its time on everything surrounding each add. Finding exactly what “surrounding” means led from object layout to reduction order, buffer ownership, and a lock I could release at the wrong moment. The earlier CPython investigation found the first piece I need here: a Python number is an object, not merely its numeric bits.
Every number is a box
The first cost is visible in the representation. In Python, the float 3.14 is not eight bytes of IEEE-754 sitting in a register. It’s a heap-allocated object:
one python float: 24 bytes (a PyObject: refcount + type pointer + the 8-byte double)Twenty-four bytes to hold eight bytes of actual number, because every Python value carries a reference count and a pointer to its type, the PyObject header from the CPython post, before it gets to the value itself. This is boxing: the real number is wrapped in an object box. And a Python list of these isn’t a block of numbers; it’s a block of pointers, each aimed at a separate float object somewhere on the heap:
list of 10,000,000 floats: 89 MB of pointers + 240 MB of float objects = ~329 MB, scatteredNow walk the loop total += x and count what each iteration actually does. Fetch the next pointer from the list. Follow it to a float object somewhere else in memory. Unbox the double out of the object. Add it to total. But total is itself a Python float object, and the result is a new number, so the addition allocates a new float object to hold it. And all of that happens through the interpreter: in the measured CPython 3.10.14 build, each += reaches an INPLACE_ADD bytecode dispatched by the eval loop, one trip around the interpreter per element. Ten million pointer fetches, ten million unbox-add-rebox cycles, ten million result allocations, and ten million bytecode dispatches surround the arithmetic.
I wanted to see the dispatch rather than assert it, so I disassembled the loop body with dis. total += x is not one instruction, it’s four:
LOAD_FAST total # push the accumulator object pointer
LOAD_FAST x # push the element object pointer
INPLACE_ADD # unbox both, add, box the result, push it
STORE_FAST total # pop it back into the accumulator slotFour bytecodes, ten million times, and the one that does the real work, INPLACE_ADD, also creates the result object. My first identity check was wrong in a useful way. I wrote 257 is 257, got True, and called it evidence that CPython cached 257. Both literals live in one code object, so the compiler is free to store one constant and load it twice. Constructing the values separately with int("257") produces two live objects and the identity test is False.
Two simultaneously live results of float(3) + 0.0 are distinct too. That establishes that the two additions did not return one shared result object. It does not establish that CPython never reuses float storage after an object dies. The allocator may reuse an address later. The claim I need is narrower: each loop iteration produces a new result before the old accumulator reference is released.
I also went looking for the “scattered” in that memory diagram, because I’d written it down as fact and wanted to see the scatter. In this CPython build, id() exposes the object’s address, so I printed it for the first eight float objects in a freshly built list:
xs[0] id = 4341692048
xs[1] id = 4341692304 gap +256
xs[2] id = 4341691216 gap -1088
xs[3] id = 4341692720 gap +1504
xs[4] id = 4354694704 gap +13,001,984
xs[5] id = 4354695824 gap +1120
xs[6] id = 4354695792 gap -32
xs[7] id = 4354696464 gap +672That is not what I expected. I’d pictured ten million objects flung uniformly across the heap, a cache miss on every access. What I actually got is mostly clustered: consecutive floats often sit a few hundred bytes apart, then one address jumps by megabytes. A fresh comprehension and the allocator’s pools can create local clusters. The earlier parenthetical “a cache miss, probably” was too confident. I did not churn and remeasure a long-lived list, so fragmentation over time remains a prediction rather than a result.
The 329 MB in that diagram bothered me too, because it is a getsizeof sum, an accounting number, not something the operating system has to agree with. My first attempt to settle it produced this:
baseline process RSS : 31 MB
after building the list of 10M floats: 441 MB (+410 MB)
after copying it into an np.array : 521 MB (+80 MB, exactly)
after `del` the list : 521 MB (unchanged)I called the last column RSS. It was not current RSS. The probe used resource.getrusage(...).ru_maxrss, and the Python documentation defines that field as the maximum resident set size. It is a high-water mark. A high-water mark cannot fall after del, even if the operating system reclaimed every page. The experiment was incapable of supporting the conclusion I drew from it.
I reran it in a fresh process and sampled both current RSS, through the process information interface, and the peak:
current RSS peak RSS
baseline 38.5 MB 39.0 MB
list of 5,000,000 floats 223.0 MB 223.0 MB
list plus float64 numpy copy 287.3 MB 287.3 MB
after del list, gc.collect(), 1 s wait 86.9 MB 287.3 MBDeleting the list reduced current RSS by 200.4 MB. The peak reduced by zero, as its definition requires. CPython’s allocator can retain arenas, and a process is not guaranteed to return every freed page, but this run did return most of them. The old measurement proved only that the process had once reached 521 MB.
The size arithmetic still teaches something, with narrower claims. Five million boxed floats plus the list of pointers raised current RSS by about 184.5 MB. The numpy allocation raised it by another 64.3 MB in this run, more than the forty million payload bytes because creating the array also needs temporary and allocator state. A list and an array store different representations, but one RSS delta is not a precise measurement of allocator overhead. To isolate payload layout, nbytes is the exact array payload, while getsizeof describes particular Python objects. RSS describes the process, including pages neither of those object-level counters assign neatly.
Even “current RSS” is not the same as bytes owned by these two variables. It counts resident pages mapped into the process. A page can contain objects from several allocations. An allocator can reserve virtual address space that has not become resident. The operating system can reclaim, compress, or share pages according to policies the Python object model does not expose. Sampling itself occurs after the allocation path has run and may catch temporary storage or unrelated interpreter activity.
The five-million-element probe reduces some of those ambiguities by running alone in a fresh process and forcing collection before its final sample. It cannot assign every resident page to one object. The useful facts are coarser: the boxed representation raised process residency far more than its forty-megabyte numeric payload, the packed array stored the same number of doubles in one forty-megabyte payload, and current residency fell materially after the boxed list became unreachable. The old high-water measurement could establish only the largest observed footprint.
There is a matching trap on allocation. Requesting eighty megabytes of virtual memory can return quickly before every page has physical backing. Filling an array writes the pages and makes that cost visible. np.empty and np.arange therefore answer different memory questions even when nbytes is identical. The probes here construct actual values before timing the sums, so the reduction is not being credited with the initial page population. A benchmark that allocates and immediately reads untouched pages would mix computation with first-touch faults.
This distinction matters later at the buffer boundary. Sharing a pointer avoids duplicating the logical payload. It does not promise that every page is resident, local to the executing CPU, warm in cache, or free of copy-on-write behavior. Zero copy is an ownership and data-movement claim at one interface. It is not a complete memory-performance model.
Python’s built-in sum(list) removed the bytecode loop and still took 17.7 ms in the original run, against 3.4 ms for numpy. It still follows PyObject* pointers and still converts Python floats to machine doubles. That comparison does not isolate boxing alone. The builtin and numpy use different iterator interfaces, type dispatch, accumulation rules, and summation algorithms. It does show that replacing Python bytecode with C while retaining the object-shaped input leaves substantial work on every element.
Then I found the result that made me sit up. I put the same ten million doubles into an array.array('d'), which stores them the way numpy does, raw eight-byte doubles packed contiguously, no per-element objects at all, and summed that with the built-in sum():
builtin sum(list of float objects) : 17.7 ms
builtin sum(array.array 'd' of raw doubles) : 41.5 msThe typed, contiguous, object-free array is more than twice as slow. I stared at that for a while, because it’s backwards from everything the rest of this post is about to argue. The reason is that sum() is a generic Python function: it pulls elements out one at a time as Python objects and adds them with the Python + protocol. A list of floats already has those objects sitting there, so sum() just borrows each pointer. But an array.array has no objects, it has raw bytes, so every time sum() asks for the next element, the array has to box a fresh float object on the spot for sum() to consume, then that object is thrown away. So the contiguous array, read through the object-shaped sum() API, pays for boxing on the fly, ten million allocations that the list had already paid for once at construction. Storing your data efficiently buys you nothing if the thing that reads it insists on objects at the door. It isn’t enough to have raw bytes; the consumer has to be willing to read raw bytes. That’s the whole reason arr.sum() and sum(arr) are different universes, and it’s the first hint of what numpy is really selling: not a better container, a consumer that never asks for an object.
The source makes the mismatch less mysterious. In CPython 3.10.14’s builtin_sum_impl, the builtin first obtains an iterator. Its specialized fast path handles exact Python integers. Once the input is floating point, the general loop repeatedly calls PyIter_Next, combines the returned object through PyNumber_Add, decrements references, and replaces the running result object.
The list iterator can return a reference to the float object already stored in its slot. The array.array implementation cannot return a pointer to a Python float because no such object exists in its element storage. Its float64 item getter reads a C double and creates a PyFloatObject for the iterator result. The efficient representation reaches an interface that requires an object, so the interface reconstructs one.
This also corrects my earlier claim that sum(list) isolated boxing and assigned the whole gap to it. It did not. The call still exercises generic iteration and Python numeric dispatch, and its accumulator follows CPython’s builtin implementation rather than numpy’s pairwise reducer. A clean experiment would need native loops with the same accumulation order and control structure, one reading boxed pointers and one reading raw doubles. I did not build that pair here. The measured comparison tells me that removing bytecode dispatch was insufficient. It does not apportion the remaining 14x among boxing, indirection, generic dispatch, and reduction design.
The array.array result exposes a more general boundary failure. Representation and interface must agree. A database column may hold packed integers and still lose the benefit if a row API materializes one language object per value. A network receive buffer may be contiguous and still get copied if the parser accepts only owned strings. A tensor view may share storage and still trigger a contiguous conversion when a kernel accepts one layout. The bytes can already be arranged well while the consumer asks the wrong question of them.
arr.sum() asks numpy for a reduction over an ndarray. That request carries dtype, dimensions, and strides into the native implementation in one call. sum(arr) asks Python for the next object until iteration ends. Both names contain “sum.” Their interfaces expose different information, so they make different execution paths possible.
The number moved, and chasing it found the accumulator
I need to come clean about something before going further, because it’s the kind of thing this blog is supposed to keep in rather than launder out. When I reran the headline loop to make the plots for this post, it didn’t come back at 289.9 ms. It came back at 106 ms. Same machine, same ten million floats, same total += x, and it was nearly three times faster than the number at the top of this page.
My first thought was thermal state or background load. I had not recorded enough about the original run to test either explanation. A 2.7x gap is too large to file under ordinary variation without evidence, so the old and new numbers have to remain two measurements whose full difference is unresolved. I did find one controlled difference. The original 289.9 ms loop summed into a variable at module level. My rerun summed into a variable inside a function. I put both forms in one process, alternated their order, and measured them again:
loop, accumulator is a function local (STORE_FAST) : 102.5 ms
loop, accumulator is a module global (STORE_GLOBAL): 179.6 msThe module form was 75 percent slower in this rerun. The dis output shows a mechanism that can explain that controlled difference. A function local is a numbered slot in the frame’s array of fast locals, so reading and writing it uses LOAD_FAST and STORE_FAST. A module global lives in the module dictionary, so the loop uses LOAD_GLOBAL and STORE_GLOBAL. The global path performs name lookup and a dictionary update on every iteration.
That does not explain the entire distance from 289.9 ms to 102.5 ms. Even the newly measured module version was 110 ms faster than the headline result. I no longer have the process state, run order, CPU state, or interpreter build record needed to reconstruct the first run. The headline remains the observation that started the investigation, not a reproducible baseline. The controlled comparison is the local against global pair because those two measurements shared the environment and changed one source property.
The fast path is a C loop over raw bytes
numpy’s answer is to stop boxing. An np.array of ten million floats is not ten million objects; it’s one contiguous block of eighty megabytes holding the raw doubles back to back, no headers, no pointers, no objects:
numpy float64 array: 80 MB, contiguous, no objects (4.1x less memory than the list, and cache-friendly)And arr.sum() doesn’t loop in Python. It makes one call into numpy and the reduction code walks the buffer in native code. The elements are adjacent, there is no PyFloatObject to unwrap for each one, and the running partial sums can stay in machine registers. There is one trip across the Python-to-C boundary instead of ten million trips around the eval loop.
I followed that call through the numpy 1.26.4 source because “a C loop” leaves too much unexplained. PyArray_Sum in calculation.c calls the generic reduction machinery with n_ops.add. The floating-point add reduction in loops_arithm_fp.dispatch.c.src selects numpy’s pairwise summation routine. That routine lives in loops_utils.h.src, where a 128-element block is accumulated through eight independent partial sums before larger inputs are recursively divided. The native path is not merely my Python loop translated to C. Its input representation and its reduction structure both changed.
That one trip isn’t free, though, and I wanted to know what it costs, because “touch Python once” is only a good deal if once is cheap. So I timed arr.sum() on arrays too small to do any real work, and found a floor:
np.ones(1).sum() : 1166 ns
np.ones(10).sum() : 1166 ns
np.ones(100).sum() : 1208 nsAbout 1.1 microseconds in the original run, and it barely moved whether the array held one element or a hundred. A July rerun on the same machine put that floor between 417 and 458 ns. The fixed work includes Python call dispatch, dtype and reduction selection, iterator setup, and construction of the returned scalar. I originally included releasing and reacquiring the GIL in that list. The numpy source contradicts me for these tiny inputs. ufunc_object.c uses NPY_BEGIN_THREADS_THRESHOLDED(count), and ndarraytypes.h defines that threshold as more than 500 elements. The one-element call does not release the GIL through that path. The floor is real, but my first attribution included work that does not occur at the sizes used to measure it.
That fixed cost means there is a size below which numpy loses to the plain Python loop, so I looked for the crossover rather than assuming it:
N= 1 loop= 208 ns numpy= 1125 ns python loop wins
N= 16 loop= 583 ns numpy= 1125 ns python loop wins
N= 64 loop= 1875 ns numpy= 1375 ns numpy wins
N=1024 loop=10708 ns numpy= 541 ns numpy winsFor an array of a handful of numbers, the boxed Python loop beat numpy in this probe because its short interpreter path cost less than numpy’s fixed call machinery. The July rerun kept the same ordering, with Python ahead through 16 elements and numpy ahead from 64, even though the absolute numpy floor was lower. This crossover is specific to the operation, numpy build, interpreter, and machine. A more expensive operation can repay the crossing sooner. A different dtype or strided input can repay it later.
I wanted to know what that native loop might be limited by, because “3.4 ms” is a number without a ceiling next to it. Eighty megabytes read in 3.4 ms is an effective payload rate of 23.5 GB/s. Warm reruns reached 59 GB/s, and the latest reached 67.8 GB/s. Those rates are derived from payload bytes divided by elapsed time. They are not hardware-counter measurements of DRAM traffic. Caches, prefetching, reduction dependencies, and loop overhead all remain possible contributors. Calling this memory-bound would require a stronger experiment, such as size scaling beyond the cache hierarchy together with counters. Changing the dtype is useful evidence, but it does not hold the executed instructions constant:
float64 80.0 MB 1.18 ms 67.8 GB/s
float32 40.0 MB 1.07 ms 37.4 GB/sHalf the payload bytes did not halve the time. That result rejects a model in which elapsed time is only payload bytes divided by one constant bandwidth. It does not identify the missing cost by itself. Then I tried float16, expecting its twenty-megabyte payload to finish first:
float16 20.0 MB 7.78 ms 2.6 GB/s # the answer is `inf`The float16 reduction was more than six times slower than float64, and its result was infinity. The overflow has a direct explanation: float16’s largest finite value is 65,504, so the requested sum does not fit in the output dtype. The performance result establishes that numpy 1.26.4 on this machine did not give this float16 reduction a fast path comparable to the measured float32 and float64 paths. I cannot infer from this timing that the M4 lacks native half-precision arithmetic, nor can I name the exact instruction sequence without disassembling or profiling numpy’s selected loop. Fewer payload bytes helped neither correctness nor speed in this case because dtype dispatch selected materially different work.
The other thing that quietly falls off the fast path is non-contiguity. Everything above assumed the eighty megabytes are packed back to back so the C loop can stride forward eight bytes at a time and the prefetcher can run ahead. Take that away and watch:
contiguous a.sum() : 1.39 ms
strided a[::2].sum() : 2.47 ms (same 10M elements, stride = 16 bytes)The strided version sums the same number of elements, but adjacent logical elements are sixteen bytes apart. It was 1.8x slower in the original run and 1.9x slower in the latest. The access pattern gives a consumer fewer adjacent values per region fetched and can require a different vector loop. I did not record cache-line events or actual memory traffic, so “moving twice the bytes across the bus” was too strong. The timing proves that this view’s layout changed the cost. It does not divide that cost between vectorization, caches, and DRAM.
I wrote the C loop myself, and it lost
Every claim so far leans on “numpy drops into a fast C loop,” and I’ve been saying it with the confidence of someone who hasn’t written the C loop. The last post in this series waved at ctypes, Cython, and pybind11 and admitted it hadn’t written one. So I wrote one. The entire fast path, the thing this whole post is worshipping, is five lines:
double csum(const double *x, long n) {
double total = 0.0;
for (long i = 0; i < n; i++) total += x[i];
return total;
}Compile that to a shared library with clang, load it with ctypes, and hand it the numpy array’s raw buffer pointer directly. First I checked that “hand it the pointer” really means zero copy and not a polite copy behind my back: arr.ctypes.data_as(...) gives back a pointer whose integer value is bit-identical to arr.__array_interface__['data'][0], the numpy buffer’s own address. Same address, no copy, my C function reads the exact eighty megabytes numpy is holding. So now I have the pure fast path, a compiled C loop over the raw buffer, no interpreter anywhere near it, and I expected it to at least match numpy. It didn’t:
ctypes csum, clang -O0 : 6.80 ms
ctypes csum, clang -O3 : 5.57 ms
numpy arr.sum() : 1.54 msMy hand-written C, at full optimization, is three and a half times slower than numpy. That stung, and I spent a while convinced I’d made a measurement error before I disassembled my own function and saw the problem in one line. Here’s the inner loop clang generated at -O3:
LBB0_5:
fadd d0, d0, d1 ; total = total + x[i], one double at a time
...fadd d0, d0, d1 is a scalar add. One double per instruction, and worse, every add reads d0 (the running total) and writes d0, so each iteration depends on the previous one. The chain is strictly serial: the adder can’t start element i+1 until element i finishes, and floating-point add has a few cycles of latency, so the loop runs at one add per several cycles and the vector units sit idle. -O3, the flag that’s supposed to make C fast, left it scalar. And it left it scalar on purpose, because it isn’t allowed to do otherwise: floating-point addition is not associative. (a + b) + c can differ from a + (b + c) in the last bit, so the compiler may not reorder or reassociate my sum without permission, and vectorizing a reduction requires reassociating it (you’d add lanes 0,4,8,… in one accumulator and 1,5,9,… in another, which changes the grouping). Strict IEEE semantics forbid that. So the correct, faithful C loop is the slow one.
The original fast-math build transformed:
ctypes csum, clang -Ofast : 1.17 ms (now faster than numpy)LBB0_5:
fadd.2d v0, v4, v0 ; four vector adds, each doing 2 doubles,
fadd.2d v1, v5, v1 ; so 8 doubles in flight across 4 independent
fadd.2d v2, v6, v2 ; accumulators, dependency chain broken
fadd.2d v3, v7, v3The generated instructions are direct evidence of what changed. The strict build contains scalar fadd; the fast-math build contains fadd.2d and multiple accumulators. Each vector add handles two doubles, and the independent accumulators shorten the dependency chain. The flag also relaxes floating-point rules beyond this one transformation, so it is not a harmless speed switch. NaNs, signed zero, infinities, rounding, and reproducibility can all matter to callers.
The speed ranking did not reproduce in July. I rebuilt the same five-line source in a temporary directory, alternated the calls, verified each result, and measured this:
ctypes csum, clang -O0 : 6.45 ms
ctypes csum, clang -O3 : 6.84 ms
ctypes csum, clang -O3 -ffast-math: 2.37 ms
numpy arr.sum() : 1.79 msFast math remained about 2.9 times faster than the strict C build, and its assembly remained vectorized. It no longer beat numpy. This is a useful contradiction, not a reason to choose the prettier run. The original probe did not preserve enough trial-level data to identify why 1.17 ms became 2.37 ms. Library caches, CPU state, compiler version, and measurement order are candidates, not answers. The stable finding across both runs is structural: the strict source compiled to a serial scalar reduction, while the relaxed source compiled to a vector reduction. The claim that my fast-math loop is faster than numpy belongs only to the historical run.
Numpy reached its own fast path by choosing the grouping in source. It does not sum left to right. Its pairwise routine handles small blocks through eight partial accumulators and recursively divides larger ranges. Within the same numpy build, dtype, shape, strides, and execution path, this gives a fixed grouping. “Bit-reproducible” across numpy versions, compiler builds, CPU features, or layouts would be a much broader promise, and the source does not make it.
One was too small to reach the accumulator
The need to choose a grouping comes from the way a double stores a number. Float64 has a sign, an exponent, and 53 bits of significand precision when the implicit leading bit is counted. It can represent every integer up to 2^53. Beyond that, representable values spread apart.
At 10^16, adjacent float64 values are two units apart:
>>> import math
>>> math.ulp(1e16)
2.0ulp here is the distance to the next representable value at that magnitude. Adding one to 10^16 produces a mathematical result halfway between representable floats. The rounding rule chooses a representable neighbor, and the one disappears:
>>> 1e16 + 1.0
1e16That makes three additions change with grouping:
>>> (1e16 + 1.0) + (-1e16)
0.0
>>> (1e16 + -1e16) + 1.0
1.0The real-number expression contains the same terms and has sum one. The first order loses the one when the large positive value is still in the accumulator. Subtracting the large value later cannot recover information that was rounded away. The second order cancels the large values first, leaving an accumulator small enough to represent the one.
This is why the compiler cannot generally turn one accumulator into four and promise the same bits. Four partial sums group elements differently. They may be closer to the real-number result or farther from it, depending on the data. Strict compilation preserves the source’s observable floating-point order. Fast-math gives the optimizer permission to apply algebraic transformations that real arithmetic permits but floating-point arithmetic does not.
Pairwise summation chooses a tree rather than leaving the grouping to whichever optimization happens to fire. For eight values, a simplified tree might first form four neighboring pairs, then two pair-of-pair sums, then one final result:
x0 x1 x2 x3 x4 x5 x6 x7
\ / \ / \ / \ /
s01 s23 s45 s67
\ / \ /
s0123 s4567
\ /
totalThe longest dependency path contains three additions rather than seven. Independent branches can be computed at the same time or mapped to vector lanes. The numerical advantage on ordinary same-scale data comes from combining partial sums of roughly comparable magnitude instead of repeatedly adding a small input into an ever-growing total.
It is not a proof that pairwise always returns the closest possible float. Adversarial signs and magnitudes can still produce substantial error. It is also not the same algorithm as a compensated sum, which tracks some lost low-order information in an extra value, or math.fsum, which maintains multiple partials. Those methods make different performance and accuracy trades.
Pairwise summation was also more accurate for the next input, and this is where I have to walk back a word in the opening paragraph. I said the loop and numpy give a “bit-for-bit identical” answer, and for the specific data in this post they do: I’m summing the integers 0 through 9,999,999 as floats, and their total, 49,999,995,000,000, is comfortably under 2^53, the point past which float64 can’t represent every integer, so no rounding happens in any order and every method agrees exactly. But that’s a property of the data, not of the methods. Sum ten million random floats in [0,1) instead, where rounding actually bites, and they diverge:
naive python loop : 4999991.227965012 error vs reference: 8.1e-07
numpy pairwise : 4999991.227965828 error vs reference: 6.5e-09
math.fsum reference: 4999991.2279658215I used math.fsum as the reference because it tracks partial sums to avoid much of the rounding loss of a naive reduction. Calling it mathematically exact was too broad. Its documentation describes an accurate floating-point sum and notes that some non-Windows builds can suffer a double-rounding difference in the least significant bit. Against that reference, numpy’s result in this run was about 120 times closer than the naive left-to-right loop. The long accumulator repeatedly combines a large partial total with a much smaller addend, while the pairwise tree tends to combine values of closer magnitude. The “identical answer” in the opening was a property of the integer-valued data, not a property shared by the two methods.
And the gap isn’t fixed, it grows with the array, which is the fingerprint of the two algorithms. I summed the same random floats at increasing sizes and watched numpy’s error stay flat near zero while the naive loop’s error climbed:
N= 128 numpy err=1.4e-14 naive err=4.3e-14
N= 1024 numpy err=1.1e-13 naive err=2.3e-13
N= 8192 numpy err=0.0 naive err=1.5e-11
N= 131072 numpy err=0.0 naive err=6.8e-10These four points are consistent with the expected advantage of a pairwise tree, but they do not prove an error-growth law. The numpy source states an error bound proportional to the logarithm of the element count for its pairwise routine, while the usual worst-case bound for naive accumulation grows with the count. The zeroes in the table mean agreement with math.fsum at the precision printed by the probe. They do not mean an exact real-number result.
The 128-element value was not inferred from this noisy table. It is PW_BLOCKSIZE in the numpy 1.26.4 source. Below that size, the routine uses eight accumulators; above it, it recursively divides the range. The measured values tell me that different grouping mattered for this sample. The source tells me which grouping the selected numpy version implemented.
I originally converted the timings to cycles per element by multiplying by “roughly four billion cycles per second”:
python loop : 10.6 ns/elem ~ 42 cycles to add one number
C -O3 (scalar) : 0.56 ns/elem ~ 2.2 cycles to add one number
numpy : 0.14 ns/elem ~ 0.5 cycles to add one number
C -Ofast (SIMD) : 0.12 ns/elem ~ 0.5 cycles to add one numberI did not measure the core frequency during those runs, so the cycle column is invented precision. Dynamic frequency makes the conversion especially weak. The elapsed values can be normalized without that assumption: the Python local loop took 10.6 ns per input, strict C took 0.56 ns, numpy took 0.15 ns, and the historical fast-math run took 0.12 ns. Assembly explains why the strict and relaxed C loops have different dependency structures. It does not let a wall-clock benchmark assign every remaining nanosecond to a specific instruction.
numpy is a C struct wearing a Python costume
To make that fast path usable, numpy needs your data to already be in that raw contiguous form, and the way it describes that form is worth opening up, because it explains numpy’s other superpower. An ndarray is, underneath, a small C struct: a pointer to the data buffer, a dtype (what kind of number, how many bytes), a shape (the dimensions), and, the clever part, strides: how many bytes to step to move one element along each axis. That last field means many operations that look like they rearrange data actually rearrange nothing:
a.reshape(3,4): shares memory? True strides=(32, 8)
a.T (transpose): shares memory? True strides=(8, 32) # swapped, zero bytes copied
a[::2] (slice): shares memory? True strides=(16,) # steps by two elementsReshaping a flat array to 3×4 copies nothing, it’s the same buffer, now read with a stride of 32 bytes per row and 8 per column. Transposing copies nothing, it swaps the two strides, so “move along a row” and “move down a column” trade places, and the same bytes are now read column-major. Slicing every other element copies nothing, it’s the same buffer with a stride of 16 bytes. These are views: different (shape, strides) descriptions pointing at one shared buffer, which is why numpy is both fast and memory-thrifty. Broadcasting, adding a vector to every row of a matrix, is the same trick with a stride of zero, so an axis of length one is read over and over without being copied. When you learn that a transpose is free and a reshape is free, you’ve learned that numpy’s speed and its memory efficiency are the same fact: it operates on descriptions of buffers, and moves data only when it truly must.
There’s a seam in “free,” though, and I stepped on it in the strided-sum measurement above. A transpose is free to create, but not always free to use. When I summed a transposed 2D array, a.T.sum() came back at 1.37 ms against 1.27 ms for the original. A full reduction is free to choose a traversal compatible with its reduction semantics, so the near-equal timing is consistent with numpy finding an efficient path for both layouts. I did not inspect that selected iterator path. a[::2].sum() was 1.8x slower, so that view’s layout materially changed execution. “Views are free” is exact for construction, not for every consumer.
Copying the slice into a contiguous array might recover the faster reduction, but the copy has to be repaid. This is not a stylistic choice between tidy and untidy arrays. It is an amortization question: pay once to rearrange the values, then save something on every later operation.
I tested that decision on 2,500,000 logical float64 values. The strided view occupied every other position in a forty-megabyte source allocation. Its contiguous copy held twenty megabytes. For each repetition count, one path reduced the strided view repeatedly. The other path made a fresh contiguous copy once and reduced that copy the same number of times. I randomized the two paths across twelve trials and report medians from the second complete run:
reductions keep stride copy once, then reduce
1 0.514 ms 1.946 ms
2 1.035 ms 2.183 ms
4 2.215 ms 3.207 ms
8 4.128 ms 3.977 ms
16 8.694 ms 6.218 msOne reduction should not pay a 1.4 ms copy to save less than one millisecond. At eight reductions, the copy path reached break-even in this run. At sixteen, it saved about 2.5 ms overall. The break-even count is not an array constant. It moves with view stride, element type, operation, allocation cost, cache state, and the amount of work each later operation performs.
The first complete run also exposed the variance hidden by that neat table. Its two-reduction strided median was 2.250 ms, more than twice the second run’s 1.035 ms, although its minimum was 0.948 ms. The other repetition counts kept the same broad ranking. This is why the table supports the decision at one, eight, and sixteen much more strongly than it supports a precise copy cost. A production implementation should not encode “copy after eight operations” from one laptop. It can carry layout through the operation graph, estimate the downstream uses, and measure representative shapes on the deployed hardware.
There is another option: change the consumer instead of the data. A loop that understands stride sixteen can unroll or gather values without allocating a copy. Numpy’s general iterator does this. Whether that path beats copying depends on how often the view is consumed and whether later consumers share the same layout capability. Copying is attractive when it converts several later operations to their preferred layout. It is wasteful when the next operation could have used the original view efficiently or when only one pass remains.
That is why the ndarray carries strides rather than a Boolean named contiguous. Layout is part of the value presented to the operation. A transpose can be cheap for a full reduction and expensive for a row-wise kernel. A zero-stride broadcast can avoid storage and create repeated reads. A negative stride can reverse an axis by placing the logical first element near the end of the allocation. The consumer decides which of those descriptions it can execute directly.
The pointer was not enough
If an array can be reduced by a pointer and a length, another library should be able to consume it without copying. My first version used ctypes and extracted a pointer through a numpy-specific interface. That proves one call can share one allocation. It does not define what the callee should do with a two-dimensional view, a byte-swapped dtype, a read-only mapping, or an exporter that dies while native code still holds the pointer.
CPython’s buffer protocol exists to carry that missing contract. An exporting object fills a Py_buffer. The structure contains a starting address, total byte length, element format, item size, dimension count, shape, and strides. It also contains obj, a reference to the exporting Python object. A consumer asks for the fields it needs by passing flags to PyObject_GetBuffer, then eventually calls PyBuffer_Release.
A Python memoryview exposes the same information without writing C. A three by four float64 array reported:
format = 'd'
itemsize = 8
ndim = 2
shape = (3, 4)
strides = (32, 8)
contiguous = True
readonly = Falseformat='d' identifies a native C double. itemsize=8 gives its width. Shape answers how many logical elements exist along each axis. Strides answer how many bytes the address changes when an index on that axis increases by one. These fields describe the view. They do not transfer ownership of the allocation and they do not promise that a consumer understands every layout they can express.
That last distinction was hidden by my ctypes example. It accepted a double * and a count, which silently assumes native-endian float64 values packed in one contiguous interval. Passing a[::2] through the same pointer would be wrong. The first value is at the advertised address, but the second logical value is sixteen bytes later. My loop advances eight bytes and reads the element the view meant to skip.
I replaced the bare ctypes call with a small CPython extension so the boundary itself could reject inputs it cannot interpret. The acquisition path is:
Py_buffer view = {0};
if (PyObject_GetBuffer(
exporter,
&view,
PyBUF_FORMAT | PyBUF_ND | PyBUF_STRIDES) < 0) {
return NULL;
}
if (view.format == NULL || strcmp(view.format, "d") != 0) {
PyBuffer_Release(&view);
PyErr_SetString(PyExc_TypeError,
"expected a native float64 buffer");
return NULL;
}
if (!PyBuffer_IsContiguous(&view, 'C')) {
PyBuffer_Release(&view);
PyErr_SetString(PyExc_BufferError,
"expected a C-contiguous buffer");
return NULL;
}The requested flags matter. PyBUF_FORMAT asks the exporter to supply the element format. PyBUF_ND asks for dimensions and shape. PyBUF_STRIDES asks for strides. A consumer that requests only a simple byte span cannot later assume it received all the metadata required to interpret an arbitrary numeric array.
The error branches matter too. Once PyObject_GetBuffer succeeds, the view is acquired. Every later path must release it exactly once. Returning after a failed dtype check without PyBuffer_Release would leak the export and its reference to the exporter. Releasing twice would violate the same ownership contract in the other direction.
The extension accepted one million contiguous float64 values and returned 499999500000.0, equal to the known arithmetic result and to numpy for this exactly representable input. The strided view was rejected:
native-sum elements=1000000 result=499999500000.0
strided-rejected stride=16 reason=expected a C-contiguous bufferI then passed values that differed along one contract dimension at a time:
accepted array.array('d') format=d strides=(8,) sum=6.0
accepted read-only numpy f64 format=d strides=(8,) sum=6.0
accepted empty numpy f64 format=d strides=(8,) sum=0.0
rejected numpy float32 format=f strides=(4,)
rejected bytes format=B strides=(1,)
rejected big-endian float64 format=>d strides=(8,)
rejected reversed float64 format=d strides=(-8,)
rejected Fortran-order 3 x 4 format=d strides=(8, 24)
rejected broadcast 3 x 4 format=d strides=(0, 8)The array.array result is important because it shows that this is not a numpy extension. The producer is a standard-library container. It exports the same native-double format and contiguous layout the consumer requested, so the native loop can read it without knowing its concrete Python type.
The read-only array is accepted because the loop only reads. If the extension requested PyBUF_WRITABLE, that exporter would reject acquisition. Checking view.readonly after requesting a read-only-capable view can also work, but asking for the required permission during acquisition produces the cleaner contract. A function that writes through view.buf after accepting a read-only exporter would have a correctness bug even if the operating system happened not to fault.
The empty array is a useful edge case. Its byte length is zero and the loop performs no dereference, so a null or otherwise non-dereferenceable buffer address is harmless. Code that reads values[0] before checking the element count would turn a mathematically ordinary empty reduction into an invalid memory access.
Float32 and bytes have contiguous storage but the wrong element representation. Reading four-byte floats with a double * would group pairs of elements into invented eight-byte values. Reading bytes as doubles would add bit patterns that were never numbers. Contiguity answers where the next byte lives. It says nothing about what those bytes mean.
The big-endian array also has eight-byte elements and stride eight. Its >d format says the byte order differs from the machine-native d expected by the C load. Ignoring that prefix would produce nonsense values on this little-endian host. A general consumer could swap the bytes or select a conversion loop. This one rejects the format because a silent conversion would violate its zero-copy premise and a raw load would be wrong.
The reversed view is compact in the sense that its elements occupy one interval, but its logical traversal begins at the high address and advances by negative eight. PyBuffer_IsContiguous(view, 'C') rejects it for the forward loop I wrote. The Fortran-order matrix stores its first axis contiguously and uses a different order for the second. It is a valid contiguous matrix in Fortran order, but not in the C order requested by the loop. Flattening it by increasing address would visit a different logical sequence than C-order flattening. Addition makes that distinction look harmless for finite values, yet order affects floating-point rounding and would be observable immediately for a non-commutative operation.
The broadcast view has a zero stride. Twelve logical positions refer to four stored numbers. A consumer that divides view.len by eight sees the physical span, not necessarily the logical element count implied by shape. My narrow consumer rejects it. Supporting it would require nested index traversal or a general iterator rather than one pointer increment.
One more assumption was absent from the first extension: alignment. Casting view.buf to double * asks the C implementation and the target architecture to treat the address as suitably aligned for a double. Shape, format, and contiguity do not independently prove that. I added a check that a nonempty address is divisible by _Alignof(double). Some processors can execute misaligned loads with a penalty. The C-level assumption still needs to be satisfied before dereferencing a double *; relying on a particular processor’s tolerance does not repair the language contract.
Rejecting the view is not a limitation to hide. It is the difference between a small consumer with a precise contract and a small consumer that sometimes reads the wrong values. I then made the conversion explicit:
strided = values[::2]
copied = np.ascontiguousarray(strided)
assert not np.shares_memory(strided, copied)
result = native.sum_buffer(copied)The copy contained four million bytes in this probe. The call site now exposes the cost and the ownership change instead of letting a wrapper make an undocumented temporary. A general reduction could implement arbitrary strides directly. That would require iterating by the supplied shape and stride arrays, dealing with negative and zero strides, checking format and byte order, and deciding how much vector performance to trade for generality. Numpy already has that machinery. The point of this extension is to make the narrow native boundary visible, not to replace it.
The smaller doorway reached C first
The extension is much narrower than numpy. That gave me a way to separate the cost of crossing into native code from the cost of everything a mature array library does after the call. I timed the extension and values.sum() on the same contiguous arrays. Each reported value is the median cost per call across nine batches. Small arrays used tens of thousands of calls per batch so the clock did not have to resolve one hundred-nanosecond intervals individually.
elements native extension numpy
0 110.5 ns 407.7 ns
1 109.9 ns 416.5 ns
16 116.5 ns 421.6 ns
64 136.0 ns 427.7 ns
1,024 652.3 ns 549.1 ns
1,000,000 569,864.6 ns 115,283.3 nsAn empty call into my extension cost about 110 ns in this environment. It still parsed a tuple, asked numpy for a buffer, checked the format, length, contiguity, and alignment, released the GIL around an empty loop, reacquired it, released the view, allocated a Python float result, and returned. “Crossing into C” is not a fixed one-microsecond law. The specific native API and wrapper decide how much happens at the doorway.
Numpy’s empty call cost about 408 ns, close to the 417 to 458 ns floor from the separate crossover rerun. It performs more dispatch because it supports axes, output arrays, dtype selection, where masks, arbitrary dimensions, and many layouts. Those capabilities have a fixed cost even when there are no elements. The extension rejects nearly all of that problem space and reaches its loop sooner.
At sixteen elements, the extension remained about 3.6 times faster than numpy. At 1,024, numpy moved ahead. At one million, numpy was about 4.9 times faster. The narrow wrapper wins the doorway and loses the road because its loop is the serial total += values[index] reduction compiled under strict floating-point rules. Numpy pays more before the first element and executes a better large reduction afterward.
The shape is the same fixed-plus-variable model that explained Python versus numpy, now inside the native side:
total time = fixed call work + per-input execution workFor two implementations A and B, the crossover appears when A’s fixed saving is consumed by its higher per-input cost. If the simple wrapper saves about 300 ns at the call and eventually spends roughly 0.45 ns more per input, a rough division puts the crossover in the hundreds of elements. The measured change in ordering between 64 and 1,024 is consistent with that estimate. It is not a fitted performance model because the per-element cost is not constant across cache sizes and numpy changes loop strategy.
This explains why replacing Python with a C extension can win a microbenchmark and still be the wrong optimization for the workload. For tiny calls, wrapper design dominates. For large arrays, the native algorithm, vectorization, layout, and memory access dominate. Batching changes both: it pays the fixed boundary once and gives the native implementation enough elements to justify a specialized loop.
The experiment also catches a mistake in the opening explanation. The 3.4 ms numpy result was not “roughly the actual floating-point work” in any universal sense. My extension performs the same count of source-level additions and takes much longer for a million elements. Source-level operation count does not describe dependency chains, grouping, vector width, or load behavior. The arithmetic has an execution structure, and numpy changed it.
The call returned through every layer it entered
ctypes let me point at a C symbol in a shared library. The extension has to participate in CPython’s object and error conventions. That adds code, but it also makes failures travel back through Python without inventing a second error channel.
The module publishes a method table:
static PyMethodDef methods[] = {
{"sum_buffer", sum_buffer, METH_VARARGS,
"Sum a contiguous float64 buffer."},
{NULL, NULL, 0, NULL}
};
static struct PyModuleDef module = {
.m_base = PyModuleDef_HEAD_INIT,
.m_name = "asquare_native_buffer",
.m_size = -1,
.m_methods = methods
};When Python imports the compiled file, the dynamic loader resolves PyInit_asquare_native_buffer. That function returns a new module object built from this definition. METH_VARARGS tells CPython to collect positional arguments into a tuple before entering sum_buffer. Faster calling conventions exist, but this ordinary one keeps the argument path visible.
The function parses "O|n" with PyArg_ParseTuple: one required Python object and one optional Py_ssize_t repetition count. The object pointer is borrowed for the duration of argument parsing. The function does not decrement it because it did not acquire a reference there. PyObject_GetBuffer later acquires a different resource, the view, whose obj field holds its own exporter reference.
These two facts look like reference-count trivia until an error occurs. Suppose the dtype is float32. Argument parsing succeeded, buffer acquisition succeeded, format validation failed. The return path must release the acquired view, set a Python exception, and return NULL:
if (validate_double_buffer(&view) < 0) {
PyBuffer_Release(&view);
return NULL;
}Returning NULL without an exception would make CPython raise a SystemError because the native function violated the error protocol. Returning a non-null object while an exception remained set would violate it in the other direction. The helper that reports the format problem calls PyErr_SetString, so the caller receives the intended TypeError.
The success path has its own ownership transfer. PyFloat_FromDouble(total) returns a new reference. Returning that pointer gives the result reference to the Python caller. Before creating it, the extension releases the buffer view. If allocating the result float fails and PyFloat_FromDouble returns NULL, no buffer lease remains to clean up.
The inner loop uses no Python API:
static double sum_repeated(
const Py_buffer *view,
Py_ssize_t repetitions) {
const double *values = (const double *)view->buf;
Py_ssize_t count = view->len / sizeof(double);
double total = 0.0;
for (Py_ssize_t r = 0; r < repetitions; ++r) {
for (Py_ssize_t i = 0; i < count; ++i) {
total += values[i];
}
}
return total;
}Py_ssize_t is CPython’s signed size type, chosen so Python container lengths fit in it. The validation checks that byte length is divisible by sizeof(double) before this division. A negative repetition count is rejected before the loop. Very large positive counts could keep the function busy for an impractical time, and the floating result can overflow to infinity. The extension treats both as caller-level behavior except that it rejects counts below one. A production API might impose a work limit or support cancellation. This probe needs repetition only to make the GIL experiment long enough to observe.
Building the extension produced a filename ending in cpython-310-darwin.so. That suffix records an implementation and ABI target. This binary was compiled against the CPython 3.10 headers in the active environment. It is not a portable library that can be copied unchanged into CPython 3.14, PyPy, Windows, or an x86-64 Linux container. The source can be rebuilt for those targets, subject to their compilers and Python development headers.
The earlier five-line csum library exposes a plain C ABI and is easier for other languages to call. It also trusts the caller to declare the signature correctly. In ctypes, forgetting to set a 64-bit return type can make the library interpret the returned bits as a C integer. Passing an incorrect pointer type can reach the loop without dtype or layout validation. The extension trades a Python-specific build for a richer checked boundary.
Neither interface is universally better. A stable C library with thin language bindings keeps the core reusable. A CPython extension can turn Python objects into validated native views with precise exceptions. Many systems use both: a language-neutral compute library, then a small binding that owns reference counts, buffer leases, and error translation. What matters here is that the thin layer is part of the execution path. It decides whether the native loop receives valid bytes and whether every owner survives long enough.
The owner disappeared before the pointer did
Layout validation prevents reading the wrong address sequence. It does not keep those addresses valid. I could acquire a pointer, return it to Python as an integer, delete the numpy array, and retain a number that still looks like an address. Nothing in the integer remembers the owner.
The buffer protocol’s view.obj solves that lifetime problem. The official contract says the field is a new reference to the exporting object and is owned by the consumer. PyBuffer_Release decrements it. I tested this by adding hold_buffer, a function that stores an acquired Py_buffer inside a PyCapsule. Its capsule destructor releases the view:
static void lease_destructor(PyObject *capsule) {
Py_buffer *view =
PyCapsule_GetPointer(capsule, "asquare.buffer_lease");
if (view != NULL) {
PyBuffer_Release(view);
PyMem_Free(view);
}
}The capsule is a lifetime handle. It is not the allocation owner, but it owns one acquired view, and that view keeps the exporter alive. I attached a weak-reference finalizer to a sixteen-element numpy array, acquired the lease, deleted the only Python name for the array, forced garbage collection, and summed through the lease:
lease after array name deleted:
exporter finalized = False
sum through pointer = 120.0
after lease deleted:
exporter finalizer events = ['exporter-released']Deleting the name did not free the array because the live Py_buffer held its reference. Deleting the lease ran the capsule destructor, released the view, and allowed the array finalizer to run. The address stayed valid for a reason I could name and test.
This is the part that “zero copy” usually omits. Sharing bytes replaces a copy with a relationship. The consumer must know how long the producer’s allocation remains alive. It must know whether writes are allowed. It must know which layouts and formats are understood. It must release its claim even when validation or later work fails. A copied array has an obvious independent lifetime. A view is cheaper because its lifetime is coupled to something else.
The same coupling appears inside numpy. A sliced array normally stores another object in its base field so the allocation owner outlives the slice. The public ndarray attributes show shape and strides, while the numpy 1.26.4 PyArrayObject_fields definition also includes data, base, descr, dimensions, and flags. A data pointer without base would be enough to address bytes and insufficient to keep them alive.
The contract does not solve concurrent mutation. A consumer may have a valid live pointer while another thread writes the same elements. Lifetime safety prevents use after free. It does not prevent a data race at the native level, produce an atomic snapshot, or make a multistep update coherent. My extension requests a readable buffer and performs an unsynchronized reduction. Its caller must ensure the values are not being mutated concurrently. “Same allocation” and “safe simultaneous access” are different properties.
I released the lock around the wrong code first
The native loop can run without touching Python objects, so I wanted it to release the GIL. The safe order is narrower than “release the GIL in the C function”:
Py_buffer view = {0};
if (PyObject_GetBuffer(exporter, &view, flags) < 0) {
return NULL;
}
if (validate_double_buffer(&view) < 0) {
PyBuffer_Release(&view);
return NULL;
}
double total;
Py_BEGIN_ALLOW_THREADS
total = sum_repeated(&view, repetitions);
Py_END_ALLOW_THREADS
PyBuffer_Release(&view);
return PyFloat_FromDouble(total);Argument parsing, buffer acquisition, validation, buffer release, and construction of the result all use the Python C API, so they run while the thread owns the GIL. Only sum_repeated, which reads the already validated native buffer and writes C locals, sits between Py_BEGIN_ALLOW_THREADS and Py_END_ALLOW_THREADS.
My first mental version released the lock as soon as the extension function began. That would allow calls into the Python C API without the required thread state. The lock is not a decoration around “slow work.” It protects access to interpreter-managed state. The correct boundary follows the operations, not the source indentation I happen to find convenient.
I also cannot release the GIL, read the buffer, and immediately release the Py_buffer before the loop finishes. The acquired view is what keeps the exporter alive. The loop and the lifetime lease cover the same interval. This order gives the function a compact invariant: while native code can dereference view.buf, one acquired view exists and the exporter remains referenced.
Timing four worker threads suggested that numpy releases the GIL for the large reduction:
pure-python loop 1 thread: 102.6 ms 4 threads: 407.2 ms
np.sum 1 thread: 1.18 ms 4 threads: 2.01 msFour Python sums performed four units of work in about four times the single-work duration. Four numpy sums performed four units in about 1.7 times the duration. That is compatible with native overlap, but scaling alone does not prove where the lock was released. It also does not prove that memory bandwidth caused the missing scaling. Scheduler placement, shared caches, CPU frequency, and the reduction implementation can affect the result.
The extension allowed a more direct negative control. I built two entry points over the same sum_repeated function. One held the GIL. One released it. A second Python thread incremented a counter and slept for one millisecond between increments. I raised Python’s thread switch interval to one second so the counter could not gain time through an ordinary interpreter handoff during the short native call.
native loop holding GIL:
elapsed = 65.9 ms
counter increments during call = 0
same native loop with GIL released:
elapsed = 45.5 ms
counter increments during call = 37The exact durations varied on repeated runs, but the control did not: zero increments while the extension held the GIL, and more than zero while it released the GIL. Both functions returned the same expected sum. The observation establishes that another Python thread executed during the released region.
Releasing the lock does not automatically make one reduction faster. In this run the released version happened to finish sooner, but the counter thread, scheduling, and repeated memory reads changed the environment. The reason to release it is that the native function no longer needs exclusive access to interpreter state. It permits other Python work to proceed. Whether concurrent work improves total throughput depends on the resources that work shares.
Numpy makes this choice by path and input size. The source threshold above explains why I should not say “arr.sum() releases the GIL” without qualification. The one-element reduction used to measure the fixed call floor stays below the 500-element threshold. The ten-million-element reduction qualifies for the threaded region in the examined numpy 1.26.4 source. Dtype-specific object loops may need the Python API and keep the lock. The GIL property belongs to an execution path, not to the spelling of a method call.
Releasing the GIL also removes an accidental scheduling barrier around the buffer. Another Python thread can now reach code that mutates the same numpy allocation while the C loop reads it. The acquired Py_buffer keeps the allocation alive, but it does not lock its contents.
I did not add a probe that writes values concurrently and treats the resulting sum as evidence. In C, an unsynchronized read racing with a write to the same non-atomic double is a data race. The language does not promise a meaningful mixture of old and new values. On one run it might appear to produce an intermediate sum. On another architecture a load might observe a torn representation, or the compiler might transform the reader under the assumption that no data race exists. A result printed by one machine would not repair the missing synchronization.
Marking a numpy array read-only controls writes attempted through array interfaces that honor its flag. It does not turn the physical pages into an immutable snapshot in every possible alias. Another existing writable view, a native pointer, or a library with different ownership rules may still reach the allocation. Read-only access is a permission granted to one handle, not proof that no writer exists.
A caller that needs a stable reduction has several designs. It can establish exclusive ownership until the call returns. It can protect mutations and reads with the same lock, acquiring that lock before entering the extension and keeping it across the released-GIL region. It can copy into private storage, paying bytes and time to remove the shared lifetime. It can use an application-level generation scheme in which writers publish a new immutable buffer and readers retain the old generation. Which design fits depends on write frequency, array size, latency, and whether consumers can share immutable storage.
The tempting lock design has one more trap. A Python lock acquired by the caller can remain held while the extension releases the GIL. Other threads can execute Python, but a thread that needs that application lock will block. That may be correct. If the native code later calls back into Python while still holding the application lock, lock ordering can produce a deadlock. This extension makes no callbacks, so its critical region is easy to reason about.
The GIL answers whether threads may manipulate CPython interpreter state concurrently. The buffer lease answers whether the exporter stays alive. An application lock or ownership rule answers whether the bytes stay logically stable. These are three separate questions. Removing one lock from the native loop does not answer the other two.
A tensor shared the bytes, but not through that contract
The original experiment also wrapped the numpy array with torch.from_numpy. The tensor and array reported the same address, and writes through either handle were visible through the other. That establishes storage sharing for that PyTorch interop function and that input. It does not establish that every PyTorch tensor exports CPython’s buffer protocol.
I had grouped memoryview, torch.from_numpy, Arrow, and DLPack under one protocol because their diagrams all contain pointers, shapes, and strides. Similar fields do not make the ownership rules interchangeable. The CPython buffer protocol is one in-process Python C API. torch.from_numpy is a library-specific conversion. Arrow defines data layouts for columnar arrays. DLPack passes tensor metadata and ownership between array libraries, including device identity and stream coordination.
That distinction becomes concrete as soon as the producer and consumer do not share CPython. A Py_buffer contains a PyObject * reference. That reference is meaningful only inside the interpreter that owns it. DLPack instead describes a managed tensor with a device, dtype, shape, optional strides, byte offset, and a deleter callback. The deleter is how the consumer eventually tells the producer that the shared tensor is no longer needed.
The Python DLPack specification gives the exchange a deliberately constrained lifetime. The producer returns a PyCapsule named dltensor. A successful consumer renames the capsule to used_dltensor, preventing another consumer from taking the same ownership token. When the consumer finishes, the managed tensor’s deleter releases the producer-side resource.
I tested that rule without needing PyTorch or a GPU. Numpy 1.26.4 can export and import DLPack. A source array returned device (1, 0), the DLPack code for CPU device zero. np.from_dlpack(source) created a view that shared memory:
DLPack device : (1, 0)
np.shares_memory : True
consumer writable : False
write source[0] = 42 : consumer[0] became 42Numpy documents its imported result as generally read-only, which is what this run produced. Zero copy did not imply equal mutation rights. The producer remained writable, so changing it changed the bytes observed by the consumer. Code that requires an immutable snapshot must copy or arrange stronger ownership and synchronization.
I attached a finalizer to the source, deleted its Python name, and kept only the DLPack consumer. The source did not finalize and the consumer still read the value. Deleting the consumer then triggered the source finalizer. This is the same lifetime question the buffer lease answered, implemented through a different protocol’s ownership chain.
I then bypassed numpy’s convenient producer object and saved one raw capsule:
capsule = np.arange(4, dtype=np.int64).__dlpack__()
first = np.from_dlpack(exporter_returning(capsule))
second = np.from_dlpack(exporter_returning(capsule))The first import succeeded. The second raised ValueError because the capsule had already been consumed and renamed. A DLPack capsule is not a reusable description that can be handed to arbitrary consumers. It is one transfer token. A producer method can create a fresh capsule for another consumer, subject to the producer’s rules, but reusing the same capsule would make final ownership ambiguous.
DLPack has to solve one problem that the CPU buffer probe avoided. Device work can be asynchronous. A GPU operation may have been enqueued but not completed when the producer returns a tensor handle. The pointer can be valid and the values can still be unfinished. Conversely, a consumer may enqueue work and return before it has finished reading, so the producer cannot immediately recycle the allocation.
The Python exchange lets the consumer pass a stream to __dlpack__. The producer must arrange for its prior work to become visible to that consumer stream. On a stream-based GPU API, this can be implemented by recording an event in the producer stream and making the consumer stream wait for it. That preserves dependency without forcing every stream on the device to stop. The exact mechanism belongs to the array library and device backend; the protocol carries enough context for them to coordinate it.
My numpy-to-numpy probe used CPU device zero, so it exercised none of that stream machinery. CPU operations in this path were synchronous from the probe’s point of view. The experiment verifies capsule consumption, shared storage, read-only import behavior, and lifetime on CPU. It cannot support a claim about CUDA stream ordering or asynchronous device copies.
The managed tensor also distinguishes its data address from a byte offset. A view may begin somewhere inside a larger allocation. Shape and strides describe logical indexing from that starting point. A null strides field can identify a compact row-major tensor; explicit strides can describe other layouts. As with Py_buffer, a consumer still has to reject dtypes, devices, or layouts it cannot execute. Protocol compatibility means it can inspect the description. It does not mean every consumer has a kernel for every description.
The deleter closes the asynchronous ownership loop only if it runs at the right time. A consumer that calls it before its queued device reads complete creates a use-after-free on the device. A consumer that never calls it leaks the producer’s resource. Libraries normally tie the managed tensor to an internal storage object whose destruction is delayed until the consuming work and references are finished. The one-consumer capsule rule makes responsibility unambiguous enough for that transfer.
DLPack is also not a general promise that CPU bytes can become GPU bytes without a transfer. The device tuple is part of the contract precisely because a consumer must know where the allocation lives. The specification permits a copy when requested through its copy semantics and addresses synchronization when streams are involved. If the source is on CPU device zero and the requested consumer needs CUDA device zero, a pointer description cannot erase the physical and API boundary. Some systems support mappings or unified memory that make particular sharing arrangements possible. Those are properties of the allocator, devices, and framework path, not consequences of writing __dlpack__.
The GPU result belongs to an older environment
The first version of this investigation included a PyTorch MPS run from the same M4:
.to('mps') for 40 MB float32 : 1.45 ms
.to('cpu') back : 1.46 ms
CPU tensor sum : 0.47 ms
MPS tensor sum : 0.82 msThe CPU tensor and MPS tensor reported different logical data pointers, and changing the CPU tensor did not change the MPS result. Those observations support a framework-level conclusion: that conversion produced distinct storage as exposed by PyTorch. They do not reveal whether the allocations occupied different physical DRAM on a unified-memory machine. A virtual address, device handle, or framework data pointer is not a physical-memory measurement.
The timings came from PyTorch 2.9.1 in the original environment. That installation is no longer available in the repository environment, so I could not rerun the MPS probe during this revision. The numbers remain historical measured results with that version label. I cannot verify whether a current PyTorch release chooses the same path.
I also cannot attribute the 0.82 ms sum to kernel launch overhead from those four stopwatch values. I did not capture a Metal trace, separate command submission from execution, or sweep input size. The MPS sum was slower than the CPU sum for that input in that run. The transfer plus MPS sum was slower still. Those are the measured claims. A profiler would be required to divide the time among dispatch, synchronization, reduction work, and result transfer.
This changes the decision the measurement can support. It does not say a GPU is generally slower for reductions, and it says nothing about matrix multiplication. It says moving one forty-megabyte tensor to MPS for one sum lost on that software and hardware combination. Reusing an already resident tensor across many operations would amortize the transfer. Increasing the problem size might change the reduction ranking. Changing the operation changes the arithmetic performed per byte. None of those alternatives were measured here.
The old stopwatch and the new one disagree
This investigation now contains two generations of evidence. The headline, the first RSS table, the original ctypes ranking, and the MPS numbers came from the earlier run. The July 2026 rerun used the same Apple M4 laptop with arm64 macOS 15.7.7, CPython 3.10.14, numpy 1.26.4, and Apple clang 15.0.0. PyTorch 2.9.1 was available for the historical MPS run and was not available in the rerun environment.
The first probes often reported one best time after warmup. The new native call-floor probe records nine randomized batches for each implementation and reports the median per call. The copy-amortization probe records twelve randomized trials per path and reports medians, while preserving the minimum so a disturbed batch is visible. The RSS probe runs in a fresh process and records current and peak memory separately. These methods do not make the laptop quiet, but they prevent one convenient sample from silently becoming the result.
Several contradictions survived:
Python local loop historical 106 ms rerun 102.5 ms
Python module loop historical 175 ms rerun 179.6 ms
numpy float64 sum historical 1.35 ms rerun 1.18 ms
fast-math C sum historical 1.17 ms rerun 2.37 ms
numpy call floor historical 1.1 us rerun about 0.42 usThe Python local and global comparison remained close. The numpy floor improved substantially. The fast-math loop became more than twice as slow and lost its lead over numpy. I have source and assembly for the native loops, but no profiler trace or preserved system-state record that explains the cross-run shifts. Pinning source code and library versions does not pin CPU frequency, thermal state, background activity, address layout, cache state, or every dependency selected when numpy was built.
The distinction between observation and mechanism helps when numbers move. The strict C assembly contains a serial scalar accumulation in both checks. The relaxed assembly contains vector fadd.2d instructions and multiple accumulators. That mechanism reproduced even when the absolute time did not. The buffer lifetime finalizer remained silent until the lease was released. The GIL-held negative control kept the counter at zero, while the released path let it advance. Those assertions are less sensitive to a few hundred microseconds of scheduling noise than the question of which optimized reduction wins by 0.3 ms.
The headline 289.9 ms is weaker. It did not reproduce, and the original environment record is incomplete. Keeping it is honest only if the article refuses to explain away the gap. It remains the strange behavior that caused the descent. The controlled local-versus-global pair explains part of it. The rest stays unknown.
The fast path acquired an owner
The pointer that made ctypes convenient was only the middle of the real path. Before the loop could safely read it, the extension had to acquire a view, validate the dtype and layout, and retain the exporter. After the loop, it had to release the view on every path. Only the region that touched no Python state could run without the GIL.
That work is not overhead accidentally left around the arithmetic. Each part closes a failure that a bare pointer creates. Format validation prevents interpreting integers as doubles. Stride validation prevents summing skipped elements as if they were adjacent. The exporter reference prevents use after free. The release pairs one ownership acquisition with one relinquishment. The GIL boundary separates interpreter operations from native operations. DLPack reconstructs comparable guarantees for a tensor exchange that cannot depend on a borrowed PyObject *.
The original ten-million-element loop crossed those boundaries ten million times because every input and every intermediate total remained a Python object. Numpy crossed once, then used its own buffer representation and its own reduction tree. My first C loop removed the objects and still lost because it retained a serial summation order. The fast-math loop changed that order and became vectorizable, then failed to reproduce its historical win over numpy. The native extension finally made the zero-copy claim honest by specifying which bytes it accepts and who keeps them alive.
The number left Python early. The allocation did not become ownerless when it left.