feat: back the statically sized scalars with a NumPy dtype

[oberbichler]

Why

Arrays of scalars were object arrays — a pointer per element and a Python object behind each one, so every NumPy operation was an object loop. The static variants have a fixed element size and, since #82, are exactly their data and trivially copyable, which is what a dtype needs.

Measured over 10 000 elements of DD3Scalar:

	dtype=object	dtype
a + b	349.5 ns	2.5 ns	138×
a * b	354.3 ns	3.9 ns	92×
np.sqrt(a)	349.2 ns	2.8 ns	123×
np.sin(a)	366.4 ns	7.0 ns	52×
np.sum(a)	353.1 ns	3.8 ns	94×
a @ b	673.0 ns	2.7 ns	246×

NumPy picks the dtype up on its own — no opt-in needed:

a = np.array(hj.variables([1.0, 2.0, 3.0]))
a.dtype     # DD3ScalarDType
a.nbytes    # 240, contiguous

What

One dtype per static type, generated from the same template as the bindings — which is exactly what #81 was the prerequisite for. No parametric dtype is needed: the size follows from the type. The dynamic variants hold a std::vector and stay object arrays.

27 ufunc loops per type: the four arithmetic ones, negative, positive, arctan2, hypot, matmul, vecdot and 19 unary mathematical functions. All of them copy through local values, which makes them safe for unaligned data at no cost — and NumPy requires that anyway (see below).

np.dot cannot be made to work

I tried, including implementing the legacy NPY_DT_PyArray_ArrFuncs_dotfunc slot, which is still active in dtype_api.h while cast and copyswap are commented out. It is unreachable:

TypeError: This function currently only supports native NumPy dtypes
           and old-style user dtypes

np.dot is not a ufunc but an _ArrayFunctionDispatcher. It dispatches via __array_function__, which looks at the array type rather than the dtype — our arrays are plain ndarray, so the protocol never fires — and then rejects anything that is not a native or old-style dtype before consulting any slot.

Registering as an old-style dtype via PyArray_RegisterDataType would work, and I decided against it: that is the API NumPy deprecated in 2.0 and intends to remove. Building the largest new feature of a library that just modernised to C++23 on an outgoing path, and facing the same migration later with users attached, is the worse trade.

@, np.matmul and np.vecdot compute the same thing and work on object arrays too, so the three tests that used np.dot now use @. np.linalg.norm also fails, inside NumPy's own Python code; np.sqrt(np.vecdot(a, a)) replaces it. np.array(a, dtype=object) converts back at any time and restores the object behaviour including np.dot. All documented in the README.

Undocumented requirements, each found only as an error

None of these are in the headers; each cost an iteration:

1. __repr__ and __str__ are mandatory
2. a cast between the DType's own instances is mandatory
3. that cast must handle unaligned data
4. …and declare NPY_METH_SUPPORTS_UNALIGNED
5. the DType needs its own tp_new; np.dtype.__new__ must not be inherited
6. a null type object is rejected — "Not giving a type object is currently not supported"
7. one type object cannot be shared between DTypes — "Can only map one python type to DType"

6 and 7 together are why the dtype has to hang off the scalar class, and why NumPy picking it up in np.array is not optional. That is the finding that ruled out the "arithmetic first, rest later" scope: a partial loop set would leave np.sin and friends failing hard where they used to fall back to the object path (33 failing tests, measured).

A mistake I made and caught

The first version of test_dtype.py compared ufunc(array) against ufunc(scalar). Both go through the very loop under test, so the tests confirmed themselves. With subtract deliberately computing plus and arccosh computing arcsinh, all 66 passed.

They now compare against the scalar operators (a - b, a.acosh()), and the same two mutations fail four tests. Without the mutation check this PR would have shipped worthless tests for its core behaviour.

Also fixed in my own code: register_loop kept the PyArrayMethod_Spec in function-local static buffers, which every call with the same <T, TArity> would have shared. NumPy may keep pointers into the spec, so they are now heap-allocated for the lifetime of the module.

Verification

• Python 239 passed, 66 of them new
• C++ 102 test cases, 1310 assertions — Release with -Werror and under clang with ASan+UBSan
• all 34 generated dtypes checked individually for element size, dtype identity, boxing and arithmetic — 0 deviations, from D0Scalar at 8 bytes to DD16Scalar at 1224
• clang-format --Werror and ruff format --check clean
• build time unchanged at 37.4 s; the module grows 2.8 → 3.4 MB

Notes

numpy>=2.0 moves into build-system.requires — the dtype is written against the NumPy C API, so its headers are needed at build time. It was already a runtime dependency.

Deliberately not included, both mechanical follow-ups through the same machinery:

• np.power — needs a mixed (T, float64) loop
• the comparison ufuncs — output dtype is bool, so a different loop shape