// HACKER NEWS — CYBERSECURITY
Tracing np.add, all the way down
The notes for this blog post have been sitting in my drafts folder
for half a year now. I’ve done a little work on NumPy itself in the past
year. Nothing notable, but enough to have to find my way around the
source. That gave me the idea to write this, but then, as other
obligations overshadowed my NumPy contributions, it just started rotting
quietly. One or two NumPy releases later I finally picked it up again,
retraced my steps, and here we are.
Here’s the premise: np.add(a, b) might well be among the
most executed lines of numerical Python in the world, and most of us
have a working mental model that’s equivalent to “it adds the arrays, in
C, quickly”. That model is correct, but there’s a lot of machinery
between the Python call and the loop that does the adding, and I think
it’s a fun machine to take apart. So today we’ll trace a single call,
np.add(a, b) with two float64 arrays, from the
Python entry point all the way down to the SIMD kernel, reading the
actual NumPy source as we go. We’ll learn a lot, I hope!
Everything below is pinned to NumPy 2.5.2, the
current release as I write this, and all links point into that tag. The
internals move around between versions1, so
if you’re spelunking along at home, check out the matching tag. I’ll
assume you’re at least somewhat comfortable reading C, but no NumPy
internals knowledge is required, that’s what we’re here for.
Before we dive in, here’s the treasure map, so you always know where
we are:
Each of these is a section below. Let’s start at the top.
The first thing to know is that np.add is not a normal
Python function. It’s an instance of numpy.ufunc, a
C-defined type2:
A ufunc is, at its core, a bundle of inner loops. One small C
function per supported type signature, plus metadata about how many
inputs and outputs there are. np.add ships 22 of them,
though I should say that types only lists the classic ones:
loops registered the modern way (more on that distinction later) live in
an internal mapping on the ufunc that Python never sees. The one we’re
chasing today is dd->d: double, double, to double. The
whole rest of this post is about how NumPy gets from your call to that
one entry, and what happens once it’s found. Other paths may vary, it’s
a big piece of kit!
When Python sees np.add(a, b), it calls the ufunc
object. The ufunc type implements the vectorcall protocol,
so the call lands in ufunc_generic_vectorcall,
which immediately forwards to the real workhorse, ufunc_generic_fastcall.
That function is long, but it reads like a checklist, and it is
the skeleton of the whole operation. Heavily abbreviated:
Parse, check for overrides, pick a loop, run it, wrap the result.
Looks like we have a plan! Now for the interesting parts.
Before NumPy commits to doing any work, it asks the arguments whether
they’d rather do it themselves (always a good modus operandi). PyUFunc_CheckOverride
walks all inputs and outputs and looks for a non-default
__array_ufunc__ method, the protocol defined in NEP
133. If any argument has one, NumPy
calls it and returns whatever it produces, and none of the machinery
we’re going to talk about below ever runs.