numpy ndarray .real/.imag return silently wrong values on complex_mp arrays — no user-dtype hook exists in numpy

[ofloveandhate]

The problem

On a numpy array of complex_mp, the ndarray attributes return silently wrong values:

w = np.array([complex_mp(1, 2), complex_mp(3, 4)])
w.real     # the complex values themselves — NOT the real parts
w.imag     # zeros — NOT the imaginary parts

No error, no warning. np.real, np.imag, and np.angle are thin wrappers over these attributes (or over the dtype check, for angle) and inherit the lie. The scalar properties are fine (w[0].imag is exact); it is only the array attributes.

Why it cannot be fixed in our bindings — receipts

ndarray.real / .imag are C-level getsets inside numpy, gated on PyArray_ISCOMPLEX, which is a hardwired check of the type number against numpy's three built-in complex types (complex64/128/256). From numpy 2.4.x numpy/_core/src/multiarray/getset.c:

/* array_real_get */
if (PyArray_ISCOMPLEX(self)) { ret = _get_part(self, 0); ... }
else { Py_INCREF(self); return (PyObject *)self; }        /* ← our array, returned whole */

/* array_imag_get */
if (PyArray_ISCOMPLEX(self)) { ret = _get_part(self, 1); }
else { /* ... _NPY_ARRAY_ZEROED ...; read-only zeros */ }

There is no slot, flag, or protocol a registered dtype can set to declare itself complex-like — not in the legacy user-dtype API (declared a frozen dead end by NEP 40), and not in the new DType API either (NEP 42 discusses complex-storage variants as motivation, but the hook was never built). Structurally, numpy's own .real/.imag are strided views into the interleaved (re, im) storage — for a user dtype numpy would need a "give me the dtype and offset of your real component" protocol that has never existed. Every complex-like user dtype ever written (quaternions included) has this hole.

What we did about it (PR #306, ADR-0051) — "a crash is better than incorrect values"

• bertini.records.Solution overrides .real/.imag at the subclass level (a Python property can shadow the C getset), so solve results are simply correct.
• bertini._numpy_guard, installed on import bertini: wraps np.real/np.imag/np.angle so that a plain mp-complex array (or a list that would convert to one) raises a TypeError naming the right tool instead of lying; np.angle raises for every mp-complex input (it branches on dtype and dies inside arctan2 — not even a subclass property can reach it). All other inputs pass through untouched.
• Sanctioned accessors: bertini.real/imag/arg and bertini.multiprec.real/imag/arg (array-capable, full precision, riding the native ufunc loops), plus the polymorphic bertini.operators namespace.
• Documented on the "Multiprecision numbers and NumPy" docs page (python/docs/source/numpy.rst), with a pinning test (numpy_ufuncs_test.py::TestComponentAccessors::test_ndarray_real_imag_attributes_are_untrustworthy) that will start failing the day numpy grows a hook.

What remains impossible (the reason this issue exists)

The raw .real / .imag attributes on a plain, self-built ndarray of complex_mp cannot be fixed, guarded, or made to raise — there is no point in numpy's attribute path where any bertini code runs. They will keep silently returning wrong values until numpy itself provides a complex-like protocol for user dtypes.

Possible upstream path

File a numpy issue / NEP discussion proposing a complex-component protocol for the new DType API — e.g. a DType slot declaring "complex-like, with component dtype D at offsets (0, itemsize/2)", which would let _get_part construct the same strided views it builds for the builtins, or a fallback get_component(which) hook for non-interleaved layouts. Until (and unless) that exists and we port real_mp/complex_mp from the legacy registration to new-style DTypes, the guard above is the ceiling of what is achievable.

Watch items for future numpy versions: re-run the pinning test on each numpy bump; if w.imag ever starts returning correct values (or raising), the guard and the docs warning can be relaxed.

🤖 Generated with Claude Code