NumPy Source Code Architecture: The Internals

A deep dive into how NumPy 2.x actually works under the hood: the C structs behind every array, the dtype and ufunc engines, the iterator that strides through memory, the promotion rules that pick your loops, and the build pipeline that stitches it all together, written for engineers who use NumPy every day and have wondered what actually happens between import numpy as np and the result coming back.

Source: https://github.com/numpy/numpy

This guide describes the NumPy main branch (the 2.x series). Where behavior changed in 2.0 you will find a NumPy 1.x callout so you can port older code with confidence. The C snippets are simplified: the names are real, but bodies are trimmed to the lines that matter for understanding.

Table of Contents

1. Repository Structure
2. Package Layout
3. The ndarray: Anatomy of the Core Data Structure
4. The dtype System
5. Memory Layout: Strides, Views, and Copies
6. The UFunc Engine
7. Loop Selection, Casting, and Type Promotion
8. The nditer: NumPy’s Workhorse Iterator
9. Indexing: The Mapping Layer
10. Array Construction and Coercion
11. Scalar Types
12. The Build System and Code Generation
13. SIMD Vectorization
14. Random Number Generation
15. linalg and BLAS/LAPACK Integration
16. FFT: pocketfft
17. The .npy and .npz File Formats
18. Sorting, Searching, and Printing
19. Error Handling and Floating Point State
20. Duck Arrays: __array_ufunc__ and __array_function__
21. Complete Lifecycle: What Happens When You Call np.add
22. Key Algorithms Summary
23. Source Code Map

1. Repository Structure

The NumPy repository is laid out like a university campus: one central building houses the heavy machinery (the C core), satellite departments each own a subpackage (random, fft, linalg), the library keeps the reference material (docs and tests), and a utilities shed holds the build tooling. Once you know which building does what, you can navigate straight to the code you need.

Source: https://github.com/numpy/numpy

NumPy is a single installable package (pip install numpy), but the repository is a monorepo containing the Python package, the C/C++ source underneath it, the documentation site, and the developer tooling. The top level looks like this:

numpy/
├── numpy/                       # The Python package you import
│   ├── _core/                   # THE HEART: ndarray, dtype, ufunc (C + Python)
│   ├── fft/                     # pocketfft wrappers
│   ├── linalg/                  # BLAS/LAPACK wrappers
│   ├── ma/                      # Masked arrays
│   ├── matrixlib/               # Legacy np.mat
│   ├── polynomial/              # poly1d, polynomial classes
│   ├── random/                  # BitGenerators + Generator (Cython)
│   ├── testing/                 # Assertions, errstate tooling
│   ├── typing/                  # Runtime types + mypy plugin
│   └── lib/                     # npy format, format helpers, mixins
├── doc/                         # Sphinx documentation (numpy.org)
├── tools/                       # CI and packaging scripts
├── .spin/                       # `spin` developer commands
├── meson_options.txt            # Build options (CPU baseline, BLAS)
├── pyproject.toml               # meson-python build definition
└── LICENSE.txt                  # BSD 3-Clause

numpy/
├── numpy/                       # The Python package you import
│   ├── _core/                   # THE HEART: ndarray, dtype, ufunc (C + Python)
│   ├── fft/                     # pocketfft wrappers
│   ├── linalg/                  # BLAS/LAPACK wrappers
│   ├── ma/                      # Masked arrays
│   ├── matrixlib/               # Legacy np.mat
│   ├── polynomial/              # poly1d, polynomial classes
│   ├── random/                  # BitGenerators + Generator (Cython)
│   ├── testing/                 # Assertions, errstate tooling
│   ├── typing/                  # Runtime types + mypy plugin
│   └── lib/                     # npy format, format helpers, mixins
├── doc/                         # Sphinx documentation (numpy.org)
├── tools/                       # CI and packaging scripts
├── .spin/                       # `spin` developer commands
├── meson_options.txt            # Build options (CPU baseline, BLAS)
├── pyproject.toml               # meson-python build definition
└── LICENSE.txt                  # BSD 3-Clause

The layering between your code and the hardware is short and steep. Almost everything interesting happens inside one compiled extension module, _multiarray_umath, which both the array object code and the ufunc engine are compiled into:

graph TD
    U["Your code<br/>np.add(a, b)"]
    P["Python layer<br/>numpy/__init__.py, _core/*.py"]
    C["C extension: _multiarray_umath<br/>multiarray/ + umath/ + npysort/"]
    H["libm, BLAS/LAPACK,<br/>OS allocator, CPU"]

    U --> P
    P --> C
    C --> H

    style C fill:#ffe6cc,stroke:#333,stroke-width:2px

NumPy 1.x callout: In 1.x the heart of the package was numpy/core. The 2.0 release renamed it to numpy/_core (and numpy/lib was slimmed down) because downstream projects kept importing from numpy.core internals that were never public API. If you see ModuleNotFoundError: No module named 'numpy.core' the fix is to import from numpy itself; 2.x also provides numpy._core._multiarray_umath compatibility shims for a small set of legacy paths.

Two structural facts are worth internalizing before anything else:

1. There is one C extension module that matters. numpy/_core/_multiarray_umath.*.so contains the ndarray type, the dtype machinery, all ufunc loops, the nditer, and the sort code. Subpackages like random and linalg are separate Cython extensions that build on top of its C API.
2. The Python layer is deliberately thin. Files like _core/multiarray.py mostly re-export names and attach documentation strings; the implementations live in C. The real logic of the Python files you will actually enjoy reading is in code generation (section 12) and the module __getattr__ that lazily loads submodules.

2. Package Layout

Reading numpy/_core is like reading a professional kitchen’s station map: a prep station turns raw ingredients into standard portions (coercion), the stove room executes recipes in bulk (umath loops), the pass expedites orders to diners (the public Python API), and the pantry stores shared spices every station needs (common headers). Knowing which station owns which job keeps you from opening every drawer.

Source: numpy/_core/

numpy/_core/ The Core Package

_core/
├── __init__.py                  # Wires the module together, __getattr__ lazy loading
├── multiarray.py                # Re-exports from C: ndarray, flagsobject, ndarrayobject
├── umath.py                     # Re-exports ufuncs from C
├── numerictypes.py              # dtype hierarchy definitions in Python
├── fromnumeric.py               # np.sum, np.mean, ... delegating to methods
├── numeric.py                   # zeros, ones, arange-ish helpers, broadcasting helpers
├── shape_base.py                # stack, block, concatenate
├── _methods.py                  # Array method implementations (a.sum -> _methods._sum)
├── einsumfunc.py                # np.einsum pure-Python path
├── memmap.py                    # np.memmap
├── code_generators/             # Scripts that WRITE C source (section 12)
├── include/numpy/               # Public C API headers
│   ├── ndarraytypes.h           # PyArrayObject, PyArray_Descr definitions
│   ├── ndarrayobject.h          # PyArray_* function declarations
│   ├── ufuncobject.h            # PyUFuncObject + PyUFunc_* API
│   └── npy_2_compat.h           # Helpers for code building against 1.x and 2.x
└── src/                         # The C/C++ implementation
    ├── multiarray/              # The array object itself (biggest directory)
    ├── umath/                   # The ufunc engine and arithmetic loops
    ├── common/                  # Shared headers, SIMD abstraction (simd/)
    ├── npymath/                 # Portable math (npy_math, complex support)
    ├── npysort/                 # Sorting: quicksort, heapsort, radix, timsort
    └── umath/svml/              # Vendored Intel SVML (vectorized libm)

_core/
├── __init__.py                  # Wires the module together, __getattr__ lazy loading
├── multiarray.py                # Re-exports from C: ndarray, flagsobject, ndarrayobject
├── umath.py                     # Re-exports ufuncs from C
├── numerictypes.py              # dtype hierarchy definitions in Python
├── fromnumeric.py               # np.sum, np.mean, ... delegating to methods
├── numeric.py                   # zeros, ones, arange-ish helpers, broadcasting helpers
├── shape_base.py                # stack, block, concatenate
├── _methods.py                  # Array method implementations (a.sum -> _methods._sum)
├── einsumfunc.py                # np.einsum pure-Python path
├── memmap.py                    # np.memmap
├── code_generators/             # Scripts that WRITE C source (section 12)
├── include/numpy/               # Public C API headers
│   ├── ndarraytypes.h           # PyArrayObject, PyArray_Descr definitions
│   ├── ndarrayobject.h          # PyArray_* function declarations
│   ├── ufuncobject.h            # PyUFuncObject + PyUFunc_* API
│   └── npy_2_compat.h           # Helpers for code building against 1.x and 2.x
└── src/                         # The C/C++ implementation
    ├── multiarray/              # The array object itself (biggest directory)
    ├── umath/                   # The ufunc engine and arithmetic loops
    ├── common/                  # Shared headers, SIMD abstraction (simd/)
    ├── npymath/                 # Portable math (npy_math, complex support)
    ├── npysort/                 # Sorting: quicksort, heapsort, radix, timsort
    └── umath/svml/              # Vendored Intel SVML (vectorized libm)

Who Owns What

The two C directories interact through a clean seam: multiarray/ owns the container (memory, layout, indexing) while umath/ owns the computation (loops, dispatch, reductions). The nditer is the bridge: ufuncs never walk memory themselves, they always ask the iterator for (pointer, stride, count) blocks and hand those blocks to inner loops.

graph TD
    subgraph "Python surface"
        MI["multiarray.py, umath.py,<br/>fromnumeric.py, _methods.py"]
    end

    subgraph "_multiarray_umath extension"
        MA["multiarray/<br/>ndarray, dtype, nditer, indexing"]
        UM["umath/<br/>ufunc object, loops, reductions"]
    end

    subgraph "Sibling extensions"
        RD["random/"]
        LA["linalg/"]
        FT["fft/"]
    end

    MI --> MA
    MI --> UM
    RD -->|"C API"| MA
    LA -->|"C API"| MA
    FT -->|"C API"| MA
    UM -->|"uses"| MA

    style MA fill:#ffe6cc,stroke:#333
    style UM fill:#d4edda,stroke:#333

NumPy 1.x callout: 1.x eagerly imported every submodule on import numpy, which made startup noticeably slower. 2.x added lazy loading via module-level __getattr__, so heavyweight submodules are only imported when first touched.

3. The ndarray: Anatomy of the Core Data Structure

An ndarray is like a spreadsheet whose cells live on one long tape while the spreadsheet itself is just a header card describing how to walk the tape: how many rows and columns, how many bytes to skip per step in each direction, and how to decode the bytes at each stop. The card is cheap to copy; the tape is not, and nearly all of NumPy’s magic comes from that asymmetry.

Source: numpy/_core/src/multiarray/arrayobject.c

The single most important idea in NumPy, documented in the internal organization page of the developer guide, is that an array is a raw data buffer plus metadata that interprets it. The C struct behind every Python array you have ever made is this (simplified from numpy/_core/include/numpy/ndarraytypes.h):

typedef struct tagPyArrayObject_fields {
    PyObject_HEAD               /* Refcount + type pointer, like every Python object */
    char *data;                 /* Pointer to the first element (not the malloc block!) */
    int nd;                     /* Number of dimensions */
    npy_intp *dimensions;       /* Shape: elements along each dimension */
    npy_intp *strides;          /* Bytes to jump per step in each dimension */
    PyObject *base;             /* Owning object: who keeps the buffer alive? */
    PyArray_Descr *descr;       /* The dtype: how to interpret each element's bytes */
    int flags;                  /* Bitmask: C_CONTIGUOUS, ALIGNED, WRITEABLE, OWNDATA */
    PyObject *weakreflist;
} PyArrayObject_fields;

Each field maps onto something you can touch from Python:

Two allocations are involved in every array: a small, fixed-size header (the struct above, plus its shape and strides arrays) and a potentially huge data buffer on the heap (allocated by alloc.c, aligned so vectorized code can chew on it). The header is where shape, strides, and descr live. Views share the buffer and get their own header, which is the entire reason slicing is O(1):

>>> a = np.arange(12)            # 12 int64 elements, 96-byte buffer
>>> b = a[2:5]                   # New header, no new buffer
>>> b.base is a
True
>>> a[2] = 99
>>> b[0]
99

>>> a = np.arange(12)            # 12 int64 elements, 96-byte buffer
>>> b = a[2:5]                   # New header, no new buffer
>>> b.base is a
True
>>> a[2] = 99
>>> b[0]
99

graph LR
    HA["header a<br/>data -> start, shape (12,), strides (8,)"]
    HB["header b = a[2:5]<br/>data -> start + 16, shape (3,), strides (8,)<br/>base = a"]
    BUF[("96-byte heap buffer<br/>0 1 2 3 4 5 6 7 8 9 10 11")]

    HA -.points into.-> BUF
    HB -.points into.-> BUF

    style BUF fill:#ffe6cc,stroke:#333
    style HB fill:#d4edda,stroke:#333

The base chain is a directed link that keeps buffers alive: when a goes out of scope, the buffer survives because b.base holds a reference to a. An array that owns its buffer has the OWNDATA flag; slices and transposes never do. Following x.base.base.base... to the end tells you which allocation is really backing your data, and np.shares_memory(a, b) walks exactly those chains to answer “do these overlap?”.

Element access is pure arithmetic on the header fields. For any index tuple, the byte offset of the element is:

offset = i0 * strides[0] + i1 * strides[1] + ... + in * strides[n]

offset = i0 * strides[0] + i1 * strides[1] + ... + in * strides[n]

so a[3, 2] with a.strides == (32, 8) lands at byte 3*32 + 2*8 = 112 past a.data. No multiplication of sizes, no search, one add per dimension. Everything else in this guide, the iterator, the ufunc loops, indexing, is machinery for evaluating that formula fast, in bulk, for many arrays in lockstep.

4. The dtype System

A dtype is like the label on a jar of screws: not the screws themselves, but the description that tells you the size, how to count them, and how to read the markings. Every array carries exactly one such label, and nearly every decision NumPy makes (which loop runs, whether a cast is safe, how many bytes each element occupies) starts by reading it.

Source: numpy/_core/src/multiarray/descriptor.c

The runtime object behind a.dtype is PyArray_Descr. Simplified from ndarraytypes.h:

typedef struct _PyArray_Descr {
    PyObject_HEAD
    PyTypeObject *typeobj;          /* Python scalar type, e.g. np.int64 */
    char kind;                      /* Kind letter, see table below */
    char type;                      /* Legacy char code: 'l', 'd', 'V', ... */
    char byteorder;                 /* '<', '>', or '=' (native) */
    int type_num;                   /* Index into the builtin type table */
    npy_intp elsize;                /* Bytes per element */
    npy_uintp alignment;            /* Alignment requirement in bytes */
    PyArray_ArrayDescr *subarray;   /* For dtypes of the form (base, shape) */
    PyObject *fields;               /* Structured: name -> (descr, offset, title) */
    PyObject *names;                /* Ordered list of field names */
} PyArray_Descr;

typedef struct _PyArray_Descr {
    PyObject_HEAD
    PyTypeObject *typeobj;          /* Python scalar type, e.g. np.int64 */
    char kind;                      /* Kind letter, see table below */
    char type;                      /* Legacy char code: 'l', 'd', 'V', ... */
    char byteorder;                 /* '<', '>', or '=' (native) */
    int type_num;                   /* Index into the builtin type table */
    npy_intp elsize;                /* Bytes per element */
    npy_uintp alignment;            /* Alignment requirement in bytes */
    PyArray_ArrayDescr *subarray;   /* For dtypes of the form (base, shape) */
    PyObject *fields;               /* Structured: name -> (descr, offset, title) */
    PyObject *names;                /* Ordered list of field names */
} PyArray_Descr;

The eleven builtin kinds each answer one question: “how do I decode the bytes at data + offset?”

Three consequences of this design are worth pausing on:

1. Structured arrays are just V dtypes with a field map. descr->fields maps names to (sub-dtype, byte offset, optional title), so a['x'] is a view at a shifted pointer with the sub-dtype’s decoder. No data moves.
2. Byte order is part of the dtype, not the file. <i4 and >i4 are different PyArray_Descr instances, which is why big-endian files transparently become swapped dtypes rather than eagerly converted.
3. O dtype stores pointers, not values. Arithmetic on object arrays calls back into Python operators element by element, which is why it is 10 to 100 times slower than typed loops.

The 2.x DType Rework: Classes for dtypes

NumPy 1.x hardcoded all dtype behavior in C switch statements scattered across the codebase, which made user-defined dtypes impossible. The 2.x tree carries the result of NEPs 41, 42, and 43: dtypes are now real Python-level classes with pluggable behavior, and ufunc loops attach to them through a new abstraction called an ArrayMethod.

graph TD
    DT["PyArray_DTypeMeta<br/>type of a dtype, with slots:<br/>discover_descr_from_pyobject,<br/>common_dtype, common_instance"]
    BUILTIN["np.dtypes.Float64DType,<br/>Int32DType, ... one per builtin"]
    INST["PyArray_Descr instances<br/>np.dtype('float64') etc."]
    AM["ArrayMethod<br/>resolve_descriptors + strided inner loops<br/>(array_method.c)"]
    ARR["ndarray.descr"]

    DT --> BUILTIN
    BUILTIN -->|"singleton instances"| INST
    AM -->|"loop registered for"| BUILTIN
    ARR --> INST

    style DT fill:#ffe6cc,stroke:#333
    style AM fill:#d4edda,stroke:#333

The supporting cast in src/multiarray/:

You do not need these files for everyday work, but they explain why 2.x can host things 1.x never could: pandas’ null-bearing dtypes, StringDType, and third-party dtypes all plug into the same ArrayMethod machinery that builtin int64 does.

Datetime Deserve a Special Mention

datetime64 packs a unit (ns, us, D, … and generic) into descriptor metadata, and datetime.c implements the tricky parts: unit conversion in casting (datetime64[s] to datetime64[ms] multiplies), business-day arithmetic (datetime_busday.c, datetime_busdaycal.c), and string parsing/formatting (datetime_strings.c). NaT (not a time) is a specific bit pattern, which is why it propagates through arithmetic without special loops.

5. Memory Layout: Strides, Views, and Copies

Strides are the walking instructions a museum visitor receives: to see the next painting within a room, walk 8 feet; to reach the same spot in the next room, take the stairs 32 steps. Two visitors can receive different instruction cards for the same gallery, and that is exactly what C order and Fortran order are.

Source: numpy/_core/src/multiarray/

The strides array is the single field that makes NumPy flexible. For a 3×4 array of int64 the same 96-byte buffer can be walked two ways:

>>> a = np.arange(12).reshape(3, 4)
>>> a.strides                       # C order: last axis fastest
(32, 8)
>>> f = np.asfortranarray(a)
>>> f.strides                       # Fortran order: first axis fastest
(8, 32)

graph TB
    subgraph "One 96-byte buffer, elements 0..11"
        ROW0["0  1  2  3"]
        ROW1["4  5  6  7"]
        ROW2["8  9  10 11"]
    end

    HC["header: C order<br/>shape (3, 4), strides (32, 8)"]
    HF["header: F order<br/>shape (4, 3)... as a.T<br/>strides (8, 32)"]

    HC -.walks rows.-> ROW0
    HC -.-> ROW1
    HF -.walks columns.-> ROW0

    style HC fill:#d4edda,stroke:#333
    style HF fill:#ffe6cc,stroke:#333

Both headers describe the same memory; only the walking instructions differ. A few consequences follow immediately from the offset formula from section 3:

• Transposes are free. a.T swaps two numbers in a new header.
• a[::-1] has negative strides ((-8,) for a 1-d int64 array). Walking backwards is as cheap as walking forwards; the loop just adds a negative step.
• Broadcasting needs no copies. When shapes (3, 1) and (1, 4) combine into (3, 4), the stretched dimensions get stride 0: every step re-reads the same element, so a scalar or row is repeated without copying anything.
• Degenerate arrays may have arbitrary strides. An array with a dimension of length 0 or 1 never steps anywhere, so NumPy does not bother normalizing its strides.

Views versus Copies

The flags bitmask on every array records the guarantees loops rely on: C_CONTIGUOUS, F_CONTIGUOUS, ALIGNED (data pointer suitably aligned for the dtype), WRITEABLE, OWNDATA, and WRITEBACKIFCOPY (a temporary that flushes into the original when released; the old UPDATEIFCOPY name was removed in 1.23). Checking a.flags['C_CONTIGUOUS'] is how library code decides whether a fast contiguous path is safe.

For allocation itself, alloc.c provides NumPy’s allocator: buffers are over-aligned (vectorized loops want aligned loads) and small blocks go through a modest cache to reduce malloc traffic. The dangerous np.lib.stride_tricks.as_strided function builds an arbitrary header over a buffer with no safety checks; sliding_window_view in the same module is the safe wrapper that computes correct strides for windowed views without copying.

6. The UFunc Engine

A ufunc is like a vending machine with many internal dispensers: you press one button (np.add), the machine inspects your coins (the operand dtypes), routes to one dispenser (an inner loop compiled for exactly those types), and processes every item in one mechanical sweep. Ninety percent of “NumPy is fast” lives inside this mechanism.

Source: numpy/_core/src/umath/ufunc_object.c

A ufunc is a Python object (so np.add can be passed around, called, and given methods like .reduce) whose C struct is essentially a table of inner loops:

typedef struct {
    PyObject_HEAD
    int nin, nout;                     /* 2 inputs, 1 output for np.add */
    int nargs;                         /* nin + nout */
    int identity;                      /* PyUFunc_Zero for add, One for multiply,
                                          None for logical_or's True identity */
    PyUFuncGenericFunction *functions; /* INNER LOOP TABLE, ntypes entries */
    char *types;                       /* nargs type chars per table row */
    void **data;                       /* per-loop context, often NULL */
    int ntypes;                        /* rows in the tables */
    char *core_signature;              /* gufuncs only, see below */
    const char *name;
} PyUFuncObject;

typedef struct {
    PyObject_HEAD
    int nin, nout;                     /* 2 inputs, 1 output for np.add */
    int nargs;                         /* nin + nout */
    int identity;                      /* PyUFunc_Zero for add, One for multiply,
                                          None for logical_or's True identity */
    PyUFuncGenericFunction *functions; /* INNER LOOP TABLE, ntypes entries */
    char *types;                       /* nargs type chars per table row */
    void **data;                       /* per-loop context, often NULL */
    int ntypes;                        /* rows in the tables */
    char *core_signature;              /* gufuncs only, see below */
    const char *name;
} PyUFuncObject;

For np.add, one row of the table says: for types (double, double, double) use the function DOUBLE_add; the next row handles (float, float, float); another (int64, int64, int64); and so on for roughly 20 type combinations, with more added by the SIMD layer of section 13.

The Inner Loop Contract

Every entry in the table is a C function with the same tiny signature, the one shown to loop authors since the 1990s:

typedef void (PyUFuncGenericFunction)(
    char **args,                /* one data pointer per operand */
    npy_intp const *dimensions, /* size of the inner dimension */
    npy_intp const *steps,      /* byte stride per operand */
    void *func                  /* elementary function, used by scalar wrappers */
);

typedef void (PyUFuncGenericFunction)(
    char **args,                /* one data pointer per operand */
    npy_intp const *dimensions, /* size of the inner dimension */
    npy_intp const *steps,      /* byte stride per operand */
    void *func                  /* elementary function, used by scalar wrappers */
);

A loop does not know what an ndarray is. It receives raw pointers, a count, and byte steps, and it processes dimensions elements in a strided walk. The generated add loop for doubles, sketched from the templates in loops.h.src (the BINARY_LOOP macro comes from fast_loop_macros.h):

static void
DOUBLE_add(char **args, npy_intp const *dimensions,
           npy_intp const *steps, void *NPY_UNUSED(func))
{
    BINARY_LOOP {
        const npy_double in1 = *(const npy_double *)ip1;
        const npy_double in2 = *(const npy_double *)ip2;
        const npy_double out = in1 + in2;
        *(npy_double *)op1 = out;
    }
}

static void
DOUBLE_add(char **args, npy_intp const *dimensions,
           npy_intp const *steps, void *NPY_UNUSED(func))
{
    BINARY_LOOP {
        const npy_double in1 = *(const npy_double *)ip1;
        const npy_double in2 = *(const npy_double *)ip2;
        const npy_double out = in1 + in2;
        *(npy_double *)op1 = out;
    }
}

The loop advances ip1 += is1, ip2 += is2, op1 += os1 each iteration. When the caller (the nditer) has arranged for a contiguous block, those steps equal 8 and the loop body reduces to a tight vectorizable loop that compilers turn into AVX2/AVX-512 instructions.

What Happens on a Call

ufunc_generic_call in ufunc_object.c is the entry point for every ufunc invocation. The choreography, with each step expanded in its own section of this guide:

flowchart TD
    A["np.add(a, b)"] --> B{"Any operand defines<br/>__array_ufunc__ ?<br/>(override.c)"}
    B -->|"yes"| C["Delegate to the other library<br/>(section 20)"]
    B -->|"no"| D["Type resolution:<br/>common dtype, pick loop row<br/>(ufunc_type_resolution.c,<br/>section 7)"]
    D --> E["Build nditer over operands<br/>with cast buffers if needed<br/>(section 8)"]
    E --> F["For each block:<br/>run the selected inner loop"]
    F --> G["Check floating point flags<br/>raise/warn per errstate<br/>(extobj.c, section 19)"]
    G --> H["Allocate/return output ndarray"]

    style C fill:#e0f0ff,stroke:#333
    style F fill:#d4edda,stroke:#333,stroke-width:2px

Where Loops Come From

You will almost never write an inner loop, but knowing where they live helps when reading stack traces or benchmarking:

The last row deserves a sentence: reductions are just ufunc methods driven by the same nditer with a zero-stride output, so np.sum(a, axis=0) runs the ordinary add loop while the iterator points the output pointer at the same accumulator element for the whole inner pass.

Generalized UFuncs (gufuncs)

core_signature is what makes a ufunc generalized. A signature like (i,j),(j,k)->(i,k) says: my loop processes whole small matrices at once, and the iterator should iterate over all leading (stacking) dimensions while the core dimensions i, j, k are guaranteed contiguous inside the loop. matmul, np.linalg‘s operations, and pocketfft‘s transforms (sections 15 and 16) are gufuncs, which is why they broadcast over stacks of matrices for free.

7. Loop Selection, Casting, and Type Promotion

Picking which compiled loop runs is like ordering at a bubble tea shop that brews many cup sizes: the staff first agree with you on one size everyone at the table can drink from, then pour. Type promotion decides the size; casting rules decide whether pouring between cups is allowed at all.

Source: numpy/_core/src/umath/ufunc_type_resolution.c

When np.add(int32_array, int64_array) executes, the engine cannot just call a loop. The two operands have different dtypes, and the ufunc only has loops registered for exact type combinations. Before any arithmetic happens, three decisions must be made:

1. What common dtype can all inputs be viewed or cast into?
2. Is a loop registered for that common dtype?
3. Is each input-to-common conversion permitted by the casting rules?

The pipeline looks like this:

flowchart TD
    A["np.add(a, b)"] --> B["Collect operand dtypes<br/>e.g. int32, int64"]
    B --> C["Find common dtype<br/>common_dtype.c: int64"]
    C --> D{"Loop registered<br/>for (int64, int64, int64)?"}
    D -->|"yes"| E["Pack casts into nditer buffers<br/>(section 8)"]
    D -->|"no"| F["TypeError:<br/>ufunc did not contain a loop<br/>with signature matching types"]
    E --> G["Execute inner loop"]
    G --> H{"Caller passed dtype= ?"}
    H -->|"yes"| I["Cast result to requested dtype"]
    H -->|"no"| J["Return result"]
    I --> J

    style F fill:#f8d7da,stroke:#333
    style G fill:#d4edda,stroke:#333,stroke-width:2px

The Common dtype Algorithm

common_dtype.c implements pairwise promotion, the same logic behind np.result_type. Within the numeric kinds, each dtype sits on a ladder per kind and sign:

bool < int8 < int16 < int32 < int64
       < uint8 < uint16 < uint32 < uint64
float16 < float32 < float64 < longdouble

bool < int8 < int16 < int32 < int64
       < uint8 < uint16 < uint32 < uint64
float16 < float32 < float64 < longdouble

Promotion walks both operands up their ladders until they meet. Two classic rules fall out of the ladder walk:

• A signed and an unsigned integer of the same bit width never meet on their own rungs, so they promote to the next larger signed type: int16 + uint16 gives int32.
• Any integer promotes into a floating dtype only if that dtype can hold its range exactly enough for the rules: int64 + float32 gives float64 because float32 cannot represent all int64 values.

>>> np.result_type(np.int16, np.uint16)
dtype('int32')
>>> np.result_type(np.int64, np.float32)
dtype('float64')

>>> np.result_type(np.int16, np.uint16)
dtype('int32')
>>> np.result_type(np.int64, np.float32)
dtype('float64')

NumPy 1.x callout: In 1.x, promotion from int64 + float32 followed the same ladder, but promotion involving Python scalars was value based, which is the single biggest behavioral change in 2.0 and deserves its own subsection below.

NEP 50: Weak Scalars

In NumPy 2.0, NEP 50 replaced value-based casting with a simple rule: Python scalars are weak, array dtypes are strong. A Python int or float never dictates the result dtype; it participates only if the array side can hold it:

>>> np.int16(3) + 3          # Python int is weak, int16 wins
np.int16(6)
>>> np.float32(1.0) + 3.0    # Python float is weak, float32 wins
np.float32(4.0)
>>> np.int16(3) + 3000       # 3000 does not fit in int16
OverflowError: Python integers 3000 out of bounds for int16

>>> np.int16(3) + 3          # Python int is weak, int16 wins
np.int16(6)
>>> np.float32(1.0) + 3.0    # Python float is weak, float32 wins
np.float32(4.0)
>>> np.int16(3) + 3000       # 3000 does not fit in int16
OverflowError: Python integers 3000 out of bounds for int16

NumPy scalars follow array rules, so a 0-d NumPy value is strong and promotes normally:

>>> np.float32(1.0) + np.float64(2.0)
np.float64(3.0)

>>> np.float32(1.0) + np.float64(2.0)
np.float64(3.0)

NumPy 1.x callout: Under 1.x value-based casting, np.uint8(200) + 200 promoted to uint16 because the values overflowed uint8, while np.uint8(1) + 1 stayed uint8. Code that silently relied on value-dependent promotion must be revisited when porting to 2.0: overflow now raises instead of quietly widening.

Casting Kinds

Whether a conversion may happen at all is governed by casting levels. can_cast_table.h holds the lookup tables, and convert_datatype.c exposes the logic through PyArray_CanCastTypeTo:

These levels gate astype, np.array(..., dtype=...), assignment, and the ufunc dtype= keyword (which replaced the old sig= spelling). The same machinery decides whether the nditer’s cast buffers in the next section are legal to build.

Where Each Piece Lives

The last row is the quiet workhorse: every astype and every buffered ufunc call eventually runs a cast loop selected from dtype_transfer.c, which registers transfer functions for each (from-dtype, to-dtype, stride pattern) triple, including SIMD variants for the common cases.

8. The nditer: NumPy’s Workhorse Iterator

The nditer is like a professional tour guide who plans one route through several warehouses simultaneously, walking you past matching shelves in every building so nobody ever loses their place. Underneath every ufunc call, reduction, and broadcast assignment, this iterator computes the single lockstep walk over all operands.

Source: numpy/_core/src/multiarray/nditer_constr.c

Arithmetic over strided arrays could naively be written as nested loops per operand, but that breaks down when operands have different shapes (broadcasting), different dtypes (casting), different memory layouts (C order, Fortran order, negative strides), and when vectorized loops want aligned contiguous memory. The nditer, implemented across nditer_constr.c (construction), nditer_api.c (C accessors), and nditer_pywrap.c (the np.nditer Python wrapper), solves all of it at once.

What It Computes at Construction Time

Given the operands, the constructor computes four things up front:

1. A broadcast shape. All operands are virtually stretched to a common shape. Stretched axes get stride 0, so repeated values cost no memory and no copy.
2. An iteration order. Axes are sorted so the innermost axis is the one where the most memory-coherent operand is contiguous, letting the inner loop process long runs at full speed.
3. Axis collapsing. When several axes share a compatible stride pattern (for example both axes of a C-contiguous operand), they merge into one longer virtual axis, so the inner loop runs once over the whole span.
4. Optional cast buffers. With buffering enabled, each operand flows through a small aligned temporary buffer, refilled by cast callbacks selected from dtype_transfer.c.

Buffered Iteration in Action

The Python-level view of the same machinery:

it = np.nditer([a, b], flags=['external_loop', 'buffered'],
               op_dtypes=['float64', 'float32'],   # cast b on the fly
               op_flags=[['readonly'], ['readonly']])
for a_block, b_block in it:
    ...   # each block is a contiguous, aligned 1-d view

it = np.nditer([a, b], flags=['external_loop', 'buffered'],
               op_dtypes=['float64', 'float32'],   # cast b on the fly
               op_flags=[['readonly'], ['readonly']])
for a_block, b_block in it:
    ...   # each block is a contiguous, aligned 1-d view

Two details make this design powerful:

• Alignment for SIMD. The buffers are aligned, so even an unaligned or strided input is handed to the vectorized loop as a contiguous aligned block. This is how strided arrays still benefit from the .dispatch.c.src loops of section 13.
• external_loop mode. The iterator hands back (pointer, stride, count) triples rather than single elements, so the same loop signature described in section 6 drives both the ufunc engine and user code.

Reductions Ride the Same Machinery

umath/reduction.c implements .reduce(), .accumulate(), .reduceat(), and .outer() on top of the nditer with the NPY_ITER_REDUCE_OK flag. The reduce operand keeps stride 0 in the inner loop, so each pass accumulates into the same output element; add.reduce over an axis is literally the ordinary add loop with a zero-stride output pointer.

9. Indexing: The Mapping Layer

Indexing code is like a hotel front desk that handles three kinds of requests: show one room (an integer), walk a wing and note every nth door (a slice), or fetch a scattered guest list (arrays of room numbers). The desk either hands you a pass that points at existing rooms (a view) or compiles a new binder of photocopies (a copy).

Source: numpy/_core/src/multiarray/mapping.c

All of a[...] flows through mapping.c, which first classifies the index object:

flowchart TD
    IDX["a[index]"] --> CLS{"What is index?"}

    CLS -->|"int / slice / Ellipsis / None"| BASIC["Basic indexing"]
    CLS -->|"ndarray / list of ints"| FANCY["Advanced (fancy) indexing"]
    CLS -->|"bool array / mask"| MASK["Convert via nonzero()<br/>item_selection.c"]
    CLS -->|"string field name"| FIELD["Descriptor field lookup<br/>a['x'] or a.x"]

    BASIC --> VIEW["Pure pointer and stride math<br/>returns a VIEW, zero copy"]
    FANCY --> MAPITER["MapIter: broadcast index arrays<br/>gather/scatter, returns a COPY"]
    MASK --> MAPITER
    FIELD --> VIEW

    style VIEW fill:#d4edda,stroke:#333
    style MAPITER fill:#ffe6cc,stroke:#333

Basic Indexing: A View by Arithmetic

Basic indexing never touches data. It adjusts three things: the data pointer, the shape, and the strides:

>>> a = np.arange(12, dtype=np.int32).reshape(3, 4)
>>> b = a[1:, ::-2]
>>> b.shape, b.strides
((2, 2), (16, -8))

>>> a = np.arange(12, dtype=np.int32).reshape(3, 4)
>>> b = a[1:, ::-2]
>>> b.shape, b.strides
((2, 2), (16, -8))

For the slice 1:, the new data pointer becomes a.data + 1 * 16 (one row of four int32 values), the first dimension shrinks to 2, and the row stride stays 16. For ::-2, the column dimension keeps 2 of its 4 positions stepping backwards, so the column stride becomes -8. Nothing was copied: b is a pair of glasses looking at the same 48 bytes from a different angle. This is why transposes, reshapes (when strides allow), and slices are effectively free, and why mutating b mutates a.

Advanced Indexing: The MapIter

When any index component is an array of integers, mapping.c switches to its MapIter machinery:

1. All index arrays are broadcast against each other.
2. The iterator walks the broadcast index shape in C order.
3. At each step it gathers the selected element from the source array into a freshly allocated result.

>>> a = np.arange(12).reshape(3, 4)
>>> rows = np.array([0, 0, 2, 2])
>>> cols = np.array([0, 3, 0, 3])
>>> a[rows, cols]
array([ 0,  3,  8, 11])
>>> a[a > 9]
array([10, 11])

>>> a = np.arange(12).reshape(3, 4)
>>> rows = np.array([0, 0, 2, 2])
>>> cols = np.array([0, 3, 0, 3])
>>> a[rows, cols]
array([ 0,  3,  8, 11])
>>> a[a > 9]
array([10, 11])

A boolean mask first goes through nonzero() (implemented in item_selection.c) to become an integer index tuple, then follows the same path. This is why a[mask] always returns a copy shaped like the number of true values, never a view.

One subtlety the MapIter implements: when multiple advanced indices are separated by a slice, the broadcast dimensions move to the front of the result (a[rows, :, cols] has shape (4, 3)); when they are adjacent, they stay in place (a[rows[:, None], cols] broadcasts the two index arrays to (4, 4) first, then appends the sliced dimensions).

Assignment Goes Through Its Own Files

Reading gathers; writing scatters. arr[idx] = value is implemented in array_assign_array.c (value is an array) and array_assign_scalar.c (value is a scalar), which reuse the MapIter walk, apply broadcasting rules to the right-hand side, and apply casting rules (the same can_cast machinery from section 7) before each store.

The Neighborhood in item_selection.c

item_selection.c hosts the index-adjacent workhorses: take, put, choose, repeat, nonzero, compress, unique, and friends. Many are thin conveniences over the same gather and scatter primitives that MapIter uses, which is why a.take(idx) is often interchangeable with a[idx] for integer indices, differing mainly in mode flags such as out-of-bound handling.

10. Array Construction and Coercion

Coercing Python objects into an array is like customs inspection at an airport: the officer opens every bag, decides one category that fits all of them, repacks what can be repacked, and sends inconsistent luggage back where it came from. ctors.c is that customs hall, and most of its rules exist to make sure you never get a bag full of mixed contents without noticing.

Source: numpy/_core/src/multiarray/ctors.c

np.array(obj) looks like one function but is a dispatcher over wildly different inputs:

flowchart TD
    IN["np.array(obj, dtype=?)"] --> Q{"What is obj?"}

    Q -->|"ndarray"| A1["Copy, view, or borrow per copy= rule"]
    Q -->|"nested sequences<br/>lists, tuples"| A2["Discover shape recursively<br/>(array_coercion.c)"]
    Q -->|"Python scalars"| A3["Find builtin dtype,<br/>build 0-d array"]
    Q -->|"object with __array__"| A4["Call obj.__array__()<br/>(array_converter.c)"]
    Q -->|"bytes / buffer"| A5["frombuffer path:<br/>zero-copy view over bytes"]
    Q -->|"object with __dlpack__"| A6["Consume DLPack capsule<br/>zero-copy (dlpack.c)"]

    A2 --> D["Find common dtype across elements<br/>(find_descriptor)"]
    D --> ALLOC["Allocate empty result"]
    ALLOC --> FILL["Convert element by element<br/>via descr setitem"]
    FILL --> OUT["ndarray"]

    style A6 fill:#d4edda,stroke:#333
    style FILL fill:#ffe6cc,stroke:#333

The interesting branch is nested sequences, handled by array_coercion.c:

1. Shape discovery walks the nesting depth and each level’s length. If two sibling lists have different lengths, the input is ragged.
2. Dtype discovery peeks at the elements (a bounded number, not the whole list) to find a common dtype: ints promote toward uint if needed, floats win over ints, strings win over numbers, and anything unrecognizable falls back to object dtype.
3. Element conversion allocates the result and runs each element through the dtype’s setitem, where a non-fitting value raises its familiar ValueError.

>>> np.array([[1, 2], [3, 4, 5]])
ValueError: setting an array element with a sequence.
The requested array has an inhomogeneous shape...

>>> np.array([[1, 2], [3, 4, 5]])
ValueError: setting an array element with a sequence.
The requested array has an inhomogeneous shape...

NumPy 1.x callout: In 1.x the ragged case created an object-dtype array after a VisibleDeprecationWarning. 2.0 turned it into a hard error, which is the single most common porting friction besides promotion changes.

The copy= Contract

2.0 made copy semantics explicit:

NumPy 1.x callout: In 1.x, copy=False meant “avoid copying if possible but copy when needed”, so the same code could silently copy on one platform and not on another. The 2.0 semantics make performance claims checkable.

Fast Paths and Special Constructors

The constructors you use daily have dedicated, hand-written C implementations in multiarraymodule.c rather than going through generic coercion: array_empty, array_zeros, array_ones, array_full, array_arange, array_frombuffer, array_fromfile, array_fromiter, and the array_copy behind .copy(). np.arange fills a contiguous buffer with start + i*step arithmetic; np.linspace (in _core/function_base.py) computes the same kind of ramp in Python before one vectorized multiply.

Two neighbors worth knowing:

• src/multiarray/textreading/ is the 2.0 rewrite of np.loadtxt: a chunked parser that tokenizes CSV-like text directly into typed arrays without creating intermediate Python objects per cell.
• dlpack.c implements the DLPack protocol, the zero-copy handshake NumPy uses with PyTorch, JAX, CuPy, and others. __dlpack__ exports a memory capsule; np.from_dlpack(obj) imports one, and if the producer lives on a GPU or uses an unsupported dtype you get a clean error instead of silent copying.

11. Scalar Types

NumPy scalars are like single-serve cereal boxes: the same cereal as the big ndarray boxes, printed with the same nutrition label (the dtype), but packaged for exactly one serving. They exist so that pulling one element out of an array keeps the type information instead of collapsing to whatever Python feels like.

Source: numpy/_core/src/multiarray/scalartypes.c.src

The scalar hierarchy (generated by the templates in scalartypes.c.src and exposed through scalarapi.c) mirrors the dtype tree: np.generic at the root, then bool_, number (integer / inexact), flexible, and so on, ending in concrete leaves like np.int64 and np.complex128.

>>> x = np.arange(3)[1]         # pulling one element yields a scalar, not a 0-d array
>>> x
np.int64(1)
>>> x.dtype, x.shape, x.strides, x.ndim
(dtype('int64'), (), (), 0)
>>> x.item()                    # escape hatch to a plain Python int
1

>>> x = np.arange(3)[1]         # pulling one element yields a scalar, not a 0-d array
>>> x
np.int64(1)
>>> x.dtype, x.shape, x.strides, x.ndim
(dtype('int64'), (), (), 0)
>>> x.item()                    # escape hatch to a plain Python int
1

Scalars behave like 0-d arrays for most purposes: they carry .dtype, .itemsize, .shape, and support .astype, .tobytes, and friends. But they are distinct C types, implemented by scalarapi.c, with conversion functions in both directions (PyArray_FromScalar boxes a scalar into an array; setitem/getitem unbox them during coercion).

Two files carry the interesting behavior:

Two details routinely surprise people and are worth internalizing:

1. np.float64 subclasses Python float (and np.complex128 subclasses complex). This is why NumPy float64 scalars flow into math.sqrt and format strings without complaint, while np.float32 does not.
2. Python scalars are weak, NumPy scalars are strong, per NEP 50 in section 7. float32_arr + 1.0 stays float32; float32_arr + np.float64(1.0) promotes to float64.

NumPy 1.x callout: NEP 51 changed scalar reprs in 2.0: repr(np.float64(3.0)) is now 'np.float64(3.0)' while str remains '3.0'. The old bare-number repr is available with np.set_printoptions(legacy='1.25') while porting test suites that asserted on printed output.

12. The Build System and Code Generation

Building NumPy is like a factory with a pre-assembly shop: before any welding happens, robots stamp out racks of repeated parts from templates (the .c.src generators), and only then does the main assembly line (Meson, Ninja, and the C compilers) weld everything into extension modules. Understanding the pre-assembly step is the key to reading the tree, because some files you need to read do not exist until build time.

Source: numpy/_core/code_generators/

NumPy 2.x builds with Meson driven by meson-python (pyproject.toml), with project-specific switches in meson_options.txt:

The developer workflow is the spin command (spin build, spin test, spin gdb), which replaced the old script.

NumPy 1.x callout: 1.x used setup.py / numpy.distutils and runtests.py. numpy.distutils was removed in 1.26 deprecation cycle and is gone in 2.0; downstream packages had to migrate to Meson or setuptools. The spin tool is the supported entry point for working on NumPy itself.

The Generators

Several enormous C files in the tree are outputs, produced at build time by Python scripts in code_generators/:

The template syntax (convtemplate) handles two kinds of repetition. Substitution with @name@, and repeat blocks that stamp a block once per list item:

/**begin repeat
 * #type = BYTE,  SHORT,  INT,  LONG#
 * #TYPE = npy_byte, npy_short, npy_int, npy_long#
 */
static void
@TYPE@_multiply(char **args, npy_intp const *dimensions, ...)
{
    /* one compiled function per type in the list */
}
/**end repeat**/

/**begin repeat
 * #type = BYTE,  SHORT,  INT,  LONG#
 * #TYPE = npy_byte, npy_short, npy_int, npy_long#
 */
static void
@TYPE@_multiply(char **args, npy_intp const *dimensions, ...)
{
    /* one compiled function per type in the list */
}
/**end repeat**/

When you open loops.c.src and it looks like one /**begin repeat**/ block covering dozens of dtypes, remember the built file is an order of magnitude larger. The API versioning guards are also generated: import_array() (which every C extension using NumPy calls at module init) compares the compiled table hash against cversions.txt; the dreaded “module compiled against API version 0x…” error is that check firing.

The Pipeline

flowchart LR
    P["pyproject.toml<br/>meson.build, meson_options.txt"]
    G["code_generators/<br/>expand .c.src templates"]
    C["Compile<br/>baseline + one object per<br/>SIMD dispatch target"]
    L["Link extensions:<br/>_multiarray_umath, random,<br/>linalg, fft"]
    W["wheel / editable install<br/>via meson-python"]

    P --> G --> C --> L --> W

    style G fill:#ffe6cc,stroke:#333
    style C fill:#d4edda,stroke:#333

The Cython subpackages (numpy/random/*.pyx, linalg/_umath_linalg.pyx, fft/_pocketfft_umath.pyx, _core/_scaled_float_dtype.pyx for the float8/float16 formats) go through the same Meson graph with Cython as an extra translation step.

13. SIMD Vectorization

SIMD dispatch is like a multilingual help desk: one sign at the front, several clerks behind the counter, and the one who answers is chosen by which languages your passport shows at that exact machine, checked at runtime, every time a process starts.

Source: numpy/_core/src/common/simd/

Modern CPUs process 4, 8, or 16 numbers per instruction (SSE2, AVX2, AVX-512, NEON on ARM). NumPy exploits this with the .dispatch.c.src suffix convention you met in section 6:

1. At build time, Meson compiles each dispatch file once per configured SIMD target, with each variant fenced to a specific instruction set.
2. At import time, src/common/cpu_dispatch.c probes the actual CPU (via CPUID and friends) and fills a feature bitmask, exposed as numpy._core._multiarray_umath.__cpu_features__ and printable with np.show_runtime().
3. When a ufunc first needs a loop, the dispatch table (dispatching.cpp) picks the most advanced compiled variant whose feature bits are all present, and caches the choice.

flowchart TD
    CALL["np.exp(a)"] --> DT["Dispatch table for exp loops"]
    PROBE["Runtime CPU features<br/>(cpu_dispatch.c at import)"]
    DT --> SEL{"Best match?"}
    PROBE --> SEL
    SEL -->|"AVX512F + SVML"| L1["AVX-512 + SVML loop"]
    SEL -->|"AVX2"| L2["AVX2 loop"]
    SEL -->|"baseline"| L3["SSE2 loop<br/>(always present)"]
    L1 --> CACHE["Function pointer cached<br/>for subsequent calls"]
    L2 --> CACHE
    L3 --> CACHE

    style L1 fill:#d4edda,stroke:#333
    style L3 fill:#e0f0ff,stroke:#333

The npyv Portability Layer

Writing nine variants of every loop by hand would be madness, so src/common/simd/ provides npyv_*: a set of header-only macros and inline functions that map one spelling onto every architecture (avx512/, avx2/, sse/, neon/, vsx/ for POWER, and WASM SIMD). A dispatch loop sketch, labeled as such because the real ones are macro-heavy:

/* Sketch of the shape of a dispatch loop, not literal source */
for (i = 0; i + vstep <= n; i += vstep) {
    const npyv_f32 va = npyv_load_f32(ip1 + i);
    const npyv_f32 vb = npyv_load_f32(ip2 + i);
    npyv_store_f32(op + i, npyv_add_f32(va, vb));
}
/* scalar tail for the remainder */

/* Sketch of the shape of a dispatch loop, not literal source */
for (i = 0; i + vstep <= n; i += vstep) {
    const npyv_f32 va = npyv_load_f32(ip1 + i);
    const npyv_f32 vb = npyv_load_f32(ip2 + i);
    npyv_store_f32(op + i, npyv_add_f32(va, vb));
}
/* scalar tail for the remainder */

The same source compiles against AVX2, NEON, or VSX with the compiler flag Meson chose for that variant.

Tuning and Debugging the Knobs

One famous resident of this layer: umath/svml/ vendors Intel’s Short Vector Math Library routines (Apache-2.0 licensed files) so exp, log, sin, and friends get properly vectorized float32/float64 loops, a change that made some transcendental workloads several times faster. Meanwhile loops_autovec.dispatch.c.src marks simple loops for compiler auto-vectorization, a lighter-weight path than hand-written npyv code that still beats the generic baseline.

14. Random Number Generation

The random subpackage is like a casino with two strictly separated roles: the dice factory (BitGenerators) manufactures long streams of raw 64-bit rolls, and the dealer (Generator) turns those rolls into game-specific outcomes like bell curves or die throws. Because the roles are separated, you can swap dice brands without retraining any dealers.

Source: numpy/random/

The package layout separates the two stages cleanly:

numpy/random/
├── __init__.py            # default_rng, Generator, bit generator exports
├── _pickle.py             # __reduce__ helpers for pickling generators
├── bit_generator.pyx      # BitGenerator base: next_uint64 / next_double
├── _pcg64.pyx             # PCG64 and PCG64DXSM (the default)
├── _philox.pyx            # Philox4x64, counter based, parallel friendly
├── _sfc64.pyx             # SFC64, small fast counter
├── _mt19937.pyx           # Mersenne Twister (legacy stream)
├── _generator.pyx         # Generator: normal, integers, choice, permutation
├── mtrand.pyx             # Legacy RandomState, frozen for compatibility
└── src/
    ├── distributions/     # C: ziggurat tables, bounded ints, binomial, ...
    └── legacy/            # C routines used only by RandomState

numpy/random/
├── __init__.py            # default_rng, Generator, bit generator exports
├── _pickle.py             # __reduce__ helpers for pickling generators
├── bit_generator.pyx      # BitGenerator base: next_uint64 / next_double
├── _pcg64.pyx             # PCG64 and PCG64DXSM (the default)
├── _philox.pyx            # Philox4x64, counter based, parallel friendly
├── _sfc64.pyx             # SFC64, small fast counter
├── _mt19937.pyx           # Mersenne Twister (legacy stream)
├── _generator.pyx         # Generator: normal, integers, choice, permutation
├── mtrand.pyx             # Legacy RandomState, frozen for compatibility
└── src/
    ├── distributions/     # C: ziggurat tables, bounded ints, binomial, ...
    └── legacy/            # C routines used only by RandomState

The division of labor at runtime:

flowchart LR
    SEED["seed / entropy"] --> SS["SeedSequence<br/>hashes entropy pool into<br/>initialization words"]
    SS --> BG["BitGenerator (PCG64)<br/>next_uint64(): raw 64 bits"]
    BG --> DIST["distributions.c<br/>ziggurat / Lemire / inversion"]
    DIST --> G["Generator.normal(size=n)<br/>typed output array"]

    style BG fill:#ffe6cc,stroke:#333
    style DIST fill:#d4edda,stroke:#333

• Stage 1, the bit stream. A BitGenerator exposes one primitive, next_uint64 (and a derived next_double that folds 53 bits into [0, 1)). PCG64 is the default; PCG64DXSM is the tweaked variant recommended when you generate from enormous numbers of parallel streams; Philox is counter based, so it jumps to arbitrary offsets instantly, which makes it the choice for SPMD and distributed work.
• Stage 2, the distribution. Generator methods consume raw bits and produce typed arrays. standard_normal walks a ziggurat table (precomputed staircase boundaries in src/distributions/), so a normal sample costs one comparison in the common case instead of transcendental math. integers uses masked rejection sampling (Lemire’s method) for unbiased bounds without modulo bias. Binomial, Poisson, hypergeometric, and friends each get a dedicated algorithm in the same C file, which the Cython _generator.pyx calls with the GIL released.

SeedSequence sits in front of stage 1: it hashes arbitrary entropy (integers, bytes, OS randomness) into a reproducible initialization state, supports .spawn() for independent child streams, and is what np.random.default_rng(seed) passes to PCG64 under the hood.

NumPy 1.x callout: The legacy interface still exists: np.random.seed(0) plus np.random.normal(...) is a global RandomState implemented in mtrand.pyx, whose bit-for-bit output stream is frozen forever for reproducibility. Its algorithms (Mersenne Twister, polar-method normals, modulo-based bounded integers) are kept in src/legacy/. New code should use default_rng(); the old stream will never improve or change.

15. linalg and BLAS/LAPACK Integration

np.linalg is like an interpreter’s booth: NumPy speaks array and iteration, the Fortran-descended BLAS/LAPACK world speaks column-major matrices and workspace arrays, and the booth translates in both directions so neither side ever notices the other’s grammar.

Source: numpy/linalg/

The package is a thin Python layer (linalg.py, mostly argument checking and docstrings) over one Cython module, _umath_linalg.pyx, that exposes LAPACK operations as gufuncs (section 6): solve, inv, det, slogdet, eig, eigvals, svd, cholesky, qr, lstsq, and matrix_power‘s building blocks. The gufunc machinery is what gives linalg its signature ergonomics:

>>> A = np.random.default_rng(0).normal(size=(4, 3, 3))
>>> b = np.ones((4, 3))
>>> np.linalg.solve(A, b).shape    # batched over the leading axis, for free
(4, 3)

>>> A = np.random.default_rng(0).normal(size=(4, 3, 3))
>>> b = np.ones((4, 3))
>>> np.linalg.solve(A, b).shape    # batched over the leading axis, for free
(4, 3)

flowchart TD
    CALL["np.linalg.solve(A, b)<br/>shape (4, 3, 3) and (4, 3)"] --> GUF["gufunc core signature<br/>(m,m),(m,n)->(m,n)"]
    GUF --> NDIT["nditer iterates the stack<br/>4 outer steps"]
    NDIT --> LAP["LAPACK dgesv on each<br/>contiguous 3x3 core"]
    LAP --> OUT["output (4, 3)"]

    style GUF fill:#ffe6cc,stroke:#333
    style LAP fill:#d4edda,stroke:#333

Under the gufunc, each core iteration calls into the linked LAPACK through thin C shims. Which LAPACK depends on the build:

Errors surface as np.linalg.LinAlgError: when a LAPACK routine reports a nonzero info code (singular matrix in solve, no convergence in svd), the Cython wrapper raises instead of returning garbage. np.linalg.norm, matrix_rank, and the pinv/pinv-adjacent helpers are pure Python composing the gufuncs above, which is also why they inherit batching automatically.

16. FFT: pocketfft

pocketfft is like a foldable travel telescope: a classical optical design (the FFTPACK lineage) reengineered in modern C++ to fold into a tiny package that ships inside the NumPy wheel, with no external mounts required. You get solid accuracy and decent speed without a planning framework.

Source: numpy/fft/

Since NumPy 1.17 the FFT implementation has been pocketfft, a C++ rewrite of the old FFTPACK Fortran algorithms, vendored under numpy/fft/_pocketfft_umath/ and exposed through _pocketfft_umath.pyx. The design mirrors linalg: each transform is a gufunc whose core dimension is the transformed axis, so batching over every other axis comes from the nditer, not from hand-written loops.

flowchart LR
    A["np.fft.fft(a, n=128, axis=-1)"] --> B["Python wrapper: shape/type checks<br/>numpy/fft/__init__.py"]
    B --> C["gufunc c2c core<br/>(n) -> (n)"]
    C --> D["pocketfft C++ kernel<br/>on contiguous core"]
    D --> E["complex128 result array"]

    style C fill:#ffe6cc,stroke:#333
    style D fill:#d4edda,stroke:#333

The public surface (fft, ifft, rfft, irfft, hfft, ihfft, fftn, rfftn, and the real transform helpers) maps onto a small set of kernels: complex-to-complex, real-to-complex, and the Hermitian variants, each with its own gufunc signature. Real transforms exploit symmetry to do roughly half the work of complex ones.

Two practical notes from the source layout. First, pocketfft computes plans on the fly per call (it has no FFTW-style plan cache), which keeps the code small and the first call cheap. Second, the NumPy docs themselves describe np.fft as the convenience FFT: SciPy vendors the same algorithms with more transforms, better workers/plan options, and scipy.fft is the recommended production interface when it is available.

17. The .npy and .npz File Formats

The .npy format is like a museum crate: a small provenance label glued inside the lid (the header dict describing dtype, shape, and memory order), followed by the artifacts themselves laid out in one strict arrangement (the raw buffer, byte for byte). Anyone can unpack the crate without the museum’s help, which is the point: the format is intentionally so simple that other languages implement it independently.

Source: numpy/lib/format.py

np.save(a) writes this layout:

magic:  \x93NUMPY        6 bytes
version: 1 byte major, 1 byte minor   (v1.0, v2.0 big headers, v3.0 utf-8)
hlen:   2 bytes (v1) or 4 bytes (v2+)  header length
header: a Python-literal dictionary, padded with spaces to a 64-byte boundary, ending in \n
        {'descr': '<i8', 'fortran_order': False, 'shape': (3, 4), }
data:   the raw buffer, exactly as laid out in memory

magic:  \x93NUMPY        6 bytes
version: 1 byte major, 1 byte minor   (v1.0, v2.0 big headers, v3.0 utf-8)
hlen:   2 bytes (v1) or 4 bytes (v2+)  header length
header: a Python-literal dictionary, padded with spaces to a 64-byte boundary, ending in \n
        {'descr': '<i8', 'fortran_order': False, 'shape': (3, 4), }
data:   the raw buffer, exactly as laid out in memory

Everything an array is, minus the Python object, lives in that one header line: descr encodes dtype and byte order, fortran_order records the stride convention, shape gives the dimensions. np.load parses the literal, reconstructs the dtype, and either memmaps (np.load(..., mmap_mode='r'), lib.format.open_memmap) or reads the buffer. Because the data section is verbatim memory, loading an int64 array of a billion elements is essentially one file read plus one header parse.

.npz is one level up: a ZIP archive of .npy members. np.savez writes uncompressed entries, np.savez_compressed deflates them, and np.load returns a lazy NpzFile whose members are only materialized when you access them by name, data['x']. That laziness matters for the common “save twelve arrays, reload one” workflow.

Two sharp edges are deliberate design decisions documented in format.py:

• Object dtype cannot round-trip without pickle. The header describes fixed-size elements only, so arrays of Python objects are pickled into the data section, and np.load refuses them unless you pass allow_pickle=True (the default flipped to False after pickle-file exploits circulated; treat untrusted .npy files as executable).
• Structured dtypes are serialized as the header’s descr list, which is why exotic field layouts can stretch the header past the original 64 KB v1.0 limit; v2.0 widened the length field, and v3.0 allowed non-ASCII field titles. Writers pick the smallest version that fits, and readers accept all three.

18. Sorting, Searching, and Printing

npysort is like a tailor with several needle sets: quick invisible basting (introsort) for general work, flat-felled seams (timsort) when original order among equals must survive, and an industrial jig (radix sort) that handles bulk integer fabric faster than any needle could. NumPy ships all three and picks per dtype and per requested stability.

Source: numpy/_core/src/npysort/

The sort directory holds one file per algorithm, each compiled generically over dtypes: quicksort.cpp (which is actually introsort: quicksort that falls back to heapsort when recursion goes deep and to insertion sort for short runs), heapsort.cpp, timsort.cpp, and radixsort.cpp. np.sort‘s kind= argument maps onto them:

np.argsort runs the same algorithms with an index permutation alongside the data (the arg variants are generated from the same source with index-tracking macros), and argfunc.dispatch.c.src in the multiarray directory provides SIMD scans for argmin/argmax (single pass, one running best value and index, several elements in flight).

The searching neighborhood splits across two files:

Printing: dragon4

The numbers you see in a REPL are also a source file’s work. multiarray/dragon4.c implements shortest-round-trip floating point printing (the “Dragon4” family of algorithms, written directly in NumPy rather than delegated to the C library): given a float64, emit the fewest digits that parse back to the same bits, which is why np.float64(0.1) prints as 0.1 and not a 17-digit approximation. _core/arrayprint.py builds on it with the summarization machinery (edgeitems, threshold, linewidth), dtype-specific formatters, and the legacy= switch used during the 2.0 repr migration described in section 11.

19. Error Handling and Floating Point State

NumPy’s error system is like a smoke detector with a configurable response dial: it can chirp (warn), trip the alarm (raise), or be muted, and every room (thread) keeps its own dial settings. The detector watches the CPU’s own floating point status flags, which is how one add loop can produce a thousand errors and you still only hear about them once.

Source: numpy/_core/src/umath/extobj.c

Floating point hardware defines five sticky status flags per thread: invalid, divide-by-zero, overflow, underflow, and inexact. Vectorized loops set them as a side effect of regular arithmetic (0.0/0.0 raises the invalid flag in hardware without branching). After each inner loop completes, the ufunc machinery in extobj.c polls the flags, compares them against the current error policy, and acts:

flowchart TD
    L["Inner loop finished<br/>(float64 divide)"] --> P["Poll FP status flags<br/>npy_get_floatstatus()"]
    P --> C{"Flag raised and<br/>policy for that class?"}
    C -->|"ignore"| N["Nothing"]
    C -->|"warn"| W["warnings.warn RuntimeWarning<br/>'divide by zero encountered in divide'"]
    C -->|"raise"| E["raise FloatingPointError"]
    C -->|"call / print / log"| O["User callback / stderr / logger"]

    style N fill:#e0f0ff,stroke:#333
    style W fill:#fff3cd,stroke:#333
    style E fill:#f8d7da,stroke:#333

The policy lives in thread-local storage managed by extobj.c, and np.errstate is its context manager:

>>> with np.errstate(divide='ignore'):
...     np.float64(1.0) / np.float64(0.0)
inf
>>> np.geterr()
{'divide': 'warn', 'over': 'warn', 'under': 'ignore', 'invalid': 'warn'}

>>> with np.errstate(divide='ignore'):
...     np.float64(1.0) / np.float64(0.0)
inf
>>> np.geterr()
{'divide': 'warn', 'over': 'warn', 'under': 'ignore', 'invalid': 'warn'}

Three practical consequences of this design:

1. NumPy follows IEEE 754, Python scalars do not. 1/0 in Python raises ZeroDivisionError; np.float64(1)/np.float64(0) computes inf and emits a warning. 0/0 gives nan plus the invalid warning. If you port scalar Python code to arrays, the division behavior changes deliberately.
2. Warnings appear once per location, not per element. The action goes through Python’s warnings module, whose default filter deduplicates by message and origin, so dividing a million elements by zero prints one line.
3. under defaults to ignore. Underflow to zero is routine in numerical code and almost never a bug, so it stays silent unless you ask.

NumPy 1.x callout: The mechanism above has been stable across both series; 2.0 mostly improved messages and made some integer overflow paths (NEP 50) raise instead of wrap silently, as section 7 covered.

20. Duck Arrays: __array_ufunc__ and __array_function__

Duck typing protocols are like standard electrical sockets: any appliance that fits the socket can draw power from NumPy’s grid without being manufactured by NumPy. CuPy arrays, sparse matrices, unit-bearing quantities, and xarray labels all work because NumPy explicitly steps aside when it sees a compatible plug.

Source: numpy/_core/src/umath/override.c

NEP 13: __array_ufunc__

Before a ufunc runs its machinery, override.c checks every operand for a different __array_ufunc__ implementation (comparing against ndarray‘s own). If one exists, NumPy does not compute; it delegates:

class TrackedArray(np.ndarray):
    pass  # inherits ndarray.__array_ufunc__, so no override happens

result = np.add(tracked, plain)   # ndarray path, sections 6-8

# Conceptually what override.c executes when a foreign operand appears:
cupy_result = cupy_array.__array_ufunc__(np.add, '__call__', ndarray_a, cp_b)

class TrackedArray(np.ndarray):
    pass  # inherits ndarray.__array_ufunc__, so no override happens

result = np.add(tracked, plain)   # ndarray path, sections 6-8

# Conceptually what override.c executes when a foreign operand appears:
cupy_result = cupy_array.__array_ufunc__(np.add, '__call__', ndarray_a, cp_b)

If the foreign object returns NotImplemented for a combination, and no other operand can handle it, NumPy raises TypeError rather than silently coercing: this is the rule that makes ndarray + cupy an error while cupy + ndarray works, because each library owns its own concatenation semantics.

For implementers, numpy.lib.mixins.NDArrayOperatorsMixin removes the boilerplate: it defines all 40-odd dunders (__add__, __matmul__, comparisons) as thin calls into the corresponding np.* ufunc, so a subclass that implements only __array_ufunc__ gets the full operator surface.

NEP 18: __array_function__

Ufunc delegation covers operators, but what about plain functions like np.mean or np.linalg.solve? Public NumPy functions are decorated (wired through arrayfunction_override.c and the dispatch tables in numpy/lib/overrides.py) to check their arguments first:

flowchart TD
    CALL["np.mean(x, axis=0)"] --> CHK{"Any arg implements<br/>__array_function__ ?"}
    CHK -->|"no"| FAST["Normal NumPy path<br/>(_methods._mean -> add.reduce)"]
    CHK -->|"cupy array"| DEL["x.__array_function__(<br/>np.mean, (x,), {'axis': 0})"]
    DEL --> GPU["CuPy executes its own mean<br/>on the GPU, no copying"]

    style FAST fill:#d4edda,stroke:#333
    style GPU fill:#ffe6cc,stroke:#333

The pattern is what lets libraries interoperate without explicit conversion calls: Dask arrays stay chunked, sparse arrays stay sparse, pint quantities carry units through np.sin. The dispatch check is deliberately cheap (a type-slot lookup per argument), but it is why the pure-Python wrappers around hot functions stay thin: every added Python-level work multiplies across every NumPy call in the ecosystem.

21. Complete Lifecycle: What Happens When You Call np.add

This chapter replays the whole movie in one sitting: from the moment your finger presses the button to the result landing in your variable, with every backstage crew member credited in order. If you remember one diagram from this guide, make it this one, because every earlier section is a scene in it.

Setup: a is an int32 array of shape (1000000,), b is int64 of the same shape, and you call c = np.add(a, b).

The whole trip costs a few microseconds of setup (override check, type resolution, iterator construction) and then runs at memory bandwidth. Every performance story in NumPy is about shrinking the setup or widening the inner block, and every correctness story (promotion, casting, errors, overrides) lives in the setup steps.

Architecture at a Glance

graph TB
    subgraph "Python layer"
        PY["numpy/__init__.py, _core/*.py,<br/>lib/format.py, arrayprint.py"]
    end

    subgraph "_multiarray_umath (C core)"
        MA["multiarray/<br/>ndarray, dtype, nditer,<br/>mapping, coercion"]
        UM["umath/<br/>ufunc object, type resolution,<br/>loops + dispatch, reductions"]
        NS["npysort/ + compiled_base<br/>sorts, partitions, searches"]
    end

    subgraph "Sibling extensions"
        RD["random/<br/>bit generators + distributions"]
        LA["linalg/<br/>LAPACK gufuncs"]
        FT["fft/<br/>pocketfft gufuncs"]
    end

    subgraph "System"
        SYS["libm / SVML, BLAS/LAPACK,<br/>malloc, CPU SIMD units"]
    end

    PY --> MA
    PY --> UM
    PY --> NS
    RD --> MA
    LA --> UM
    FT --> UM
    UM --> MA
    MA --> SYS
    LA --> SYS

    style UM fill:#ffe6cc,stroke:#333,stroke-width:2px
    style MA fill:#d4edda,stroke:#333,stroke-width:2px

22. Key Algorithms Summary

This table is the pocket card for the whole guide: one row per load-bearing algorithm, where it lives, and the trick that makes it work. Skim it when you need to find a file fast, or before an interview.

23. Source Code Map

The bibliography is a map of the territory: every row is a file worth opening when you want the ground truth for its topic. All links point at the main branch, and every path in this guide can be reached from one of these starting points.

This guide is based on the NumPy main branch (2.x series). The biggest shifts from 1.x worth remembering: the numpy.core to numpy._core rename, NEP 50’s weak-scalar promotion replacing value-based casting, NEP 51’s scalar reprs, ragged input becoming an error, Meson replacing numpy.distutils, and the DType/ArrayMethod rework (NEPs 41-43) that opened the door to user-defined dtypes like StringDType.