numpy compat: full ufunc coverage, sorting, safe reductions, tolerance comparisons for mp dtypes

core/CMakeLists.txt

@@ -570,6 +570,15 @@ target_precompile_headers(bertini2_exe PRIVATE
     <bertini2/parallel.hpp>
 )
 
+# clang validates a PCH by the mtime of its input; when the wheel build configures
+# twice in one job (test build + wheel build sharing bld/), cmake regenerates an
+# IDENTICAL cmake_pch.cxx with a fresh mtime and clang refuses the still-valid PCH
+# ("has been modified since the precompiled header was built: mtime changed").
+# -fno-pch-timestamp exists precisely for this (see the ccache docs on PCH).
+if(CMAKE_CXX_COMPILER_ID MATCHES "Clang")
+    target_compile_options(bertini2_exe PRIVATE "$<$<COMPILE_LANGUAGE:CXX>:SHELL:-Xclang -fno-pch-timestamp>")
+endif()
+
 # todo: this should be made a devmode thing
 #target_compile_options(bertini2 PRIVATE -Wall -Wextra)
 

docs/adr/0051-numpy-ufunc-coverage-and-owned-getitem.md

@@ -0,0 +1,161 @@
+# ADR-0051: Full numpy ufunc coverage for the mp dtypes; getitem returns owned copies
+
+**Status:** Accepted
+**Date:** 2026-07-08
+
+## Context
+
+The eigenpy-registered numpy dtypes for `real_mp` / `complex_mp` covered only the
+arithmetic core (add/subtract/multiply/divide, equality + real orderings,
+negative/square/sqrt, matmul, and the hardened `dotfunc`).  Everything else —
+`np.abs`, `np.conj`, the transcendental family, `power`, `sign`, `minimum/maximum`,
+rounding, the `isnan` predicates — raised `ufunc ... not supported`; `np.sort` /
+`np.argmax` failed on empty dtype slots ("type does not have compare function" /
+"data type not ordered"); and the docs declared the identity-seeded reductions
+(`np.sum`/`np.prod`/`np.mean`) permanently unsupported after a build-dependent
+`SystemError` (documented 2026-06-29, "not something Bertini can patch").
+
+Investigating the reductions on current numpy (2.3.2 and 2.4.6) showed the
+`SystemError` no longer reproduces — but exposed something worse hiding behind them:
+
+**eigenpy's `getitem` returns `boost::ref(slot)` — a Python scalar aliasing numpy
+array storage.**  Our heal-on-read `getitem` specializations (ADR-0006) had copied
+that behavior faithfully.  Consequences:
+
+- A scalar extracted from a *temporary* array dangles once the array is freed.
+  `s = np.sum(v)` is exactly that: the reduce result is a temporary 0-d array,
+  `s` aliased its buffer, and reading `s` later gave zeros, garbage precision, or an
+  MPFR assertion SIGABRT inside `str()`.
+- `np.mean` returned a **silently wrong `0`** through the same mechanism.
+- This is the root of the ADR-0031 / #259 hazard class ("indexing an eigenpy Vec
+  returns an aliasing view"), previously mitigated consumer-by-consumer with a
+  copy-at-extraction rule in Python code.
+
+## Decision
+
+1. **Register guarded loops for the full ufunc set** on both dtypes
+   (`python_bindings/include/eigenpy_interaction.hpp`, `registerGuardedUfunct`):
+
+   - both dtypes: `absolute` (complex → `real_mp` output), `conjugate`, `sign`
+     (numpy-2 semantics: complex sign is `z/|z|`), `positive`, `reciprocal`,
+     `power`, `exp`, `log`, `log10`, full trig/hyperbolic + inverses,
+     `isnan`/`isinf`/`isfinite` (→ bool);
+   - real only (ordering- or domain-dependent): `greater`/`less`/... , `fabs`,
+     `exp2`, `log2`, `expm1`, `log1p`, `cbrt`, `floor`, `ceil`, `trunc`, `rint`,
+     `signbit`, `arctan2`, `hypot`, `copysign`, `fmod`, `remainder`,
+     `floor_divide`, `minimum`, `maximum`, `fmin`, `fmax`.
+
+   Every loop reads through `value_or_zero` (the ADR-0006 uninitialized-slot
+   doctrine) and calls the same boost::multiprecision free function the
+   `bertini.multiprec` scalar function binds, so `np.f(a)[i] == mp.f(a[i])`
+   exactly.  Deliberate semantic choices, matching numpy's float64 behavior:
+   `rint` rounds half-to-even (direct `mpfr_rint` in `MPFR_RNDN`; boost's `rint`
+   rounds half away), `remainder`/`mod` takes the sign of the divisor (`fmod`
+   keeps C semantics), `minimum`/`maximum` propagate nan while `fmin`/`fmax`
+   ignore it.
+
+2. **Fill the `compare`/`argmax`/`argmin` dtype slots for `real_mp`**
+   (`HardenCompare`, `HardenArgMinMax`, same install pattern as `HardenDotfunc`).
+   Enables `np.sort`/`argsort`/`searchsorted`/`unique`/`median`/`argmax`/`argmin`.
+   Complex stays unordered on purpose — numpy's lexicographic complex ordering is
+   historical baggage we do not reproduce.
+
+3. **`getitem` returns an owned copy, never `boost::ref`.**  An indexed element is
+   a durable value, as numpy users expect.  This kills the ADR-0031/#259 hazard
+   class at the source (that ADR's copy-at-extraction rule in Python remains good
+   hygiene but is no longer load-bearing), and it is what makes the reductions
+   *actually* safe rather than accidentally readable.  Cost: one mp copy per
+   element read.
+
+4. **Reductions are supported and regression-tested**, on numpy ≥ 2.3 (verified
+   2.3.2 and 2.4.6; `python/test/classes/numpy_ufuncs_test.py::TestReductions`
+   pins them in CI on all three platforms).  The docs note the historical
+   `SystemError` and keep the `initial=` idiom as the fallback for older numpy.
+   `pyproject.toml` keeps `numpy` unpinned.
+
+5. **The float64 boundary stays closed for VALUES, open for tolerance
+   comparisons** (both decided 2026-07-08).  `double → mp` casts remain
+   registered *unsafe*, so float64 scalars/arrays do not silently promote into
+   mp arrays: the conversion itself is bit-exact, but a promoted float64 `0.1`
+   is not the decimal `0.1` the user typed — the user has to think.  However,
+   mixed mp-vs-float64 **ordering** loops (`<`, `<=`, `>`, `>=`, both operand
+   orders, real only) ARE registered: `np.abs(a - b) < 1e-10` is safe — the
+   result is a bool, no float flows into an mp value, the comparison is exact
+   (boost compares the number against the double directly), and it matches the
+   C++ solvers' double `ToleranceT` and the scalar `GreatLessVisitor<T,double>`
+   precedent.  Mixed EQUALITY stays unregistered (exact equality against a
+   float literal is the 0.1-intent trap; not bound at scalar level either), as
+   does mixed arithmetic; `np.isclose`/`np.allclose` still raise (they *compute*
+   with float64 tolerances internally).  `mp → complex128` casts are registered
+   unsafe alongside the pre-existing `mp → double`, so `arr.astype(complex)` /
+   `astype(float)` are the explicit, conscious truncations.  Registration-order
+   note: casts must be registered BEFORE the ufunc loops — registering the
+   mixed loops makes numpy query the mp↔double casts, and a cast first
+   registered after being queried is permanently ignored (numpy
+   RuntimeWarning).
+
+6. **Component access on complex arrays** goes through new array overloads of
+   `multiprec.real`/`imag`/`arg` (returning `real_mp` arrays).  numpy cannot know
+   a legacy user dtype is complex-like — `ndarray.real`/`.imag` are C getsets
+   gated on `PyArray_ISCOMPLEX`, a hardwired builtin-type-number check (verified
+   against numpy 2.4.x `getset.c`; neither the legacy user-dtype API nor the
+   NEP 42 new DType API offers a complex-like hook) — so on mp-complex arrays
+   `.real` returns the complex values and `.imag` returns zeros, silently.  Not
+   fixable at the attribute; defended everywhere reachable (decided 2026-07-08,
+   "a crash would be better than incorrect values"):
+
+   - `bertini.records.Solution` overrides `.real`/`.imag` at the subclass level —
+     solve results are simply correct;
+   - importing bertini wraps the module functions `np.real`/`np.imag`/`np.angle`
+     (`bertini._numpy_guard`): plain mp-complex arrays (and lists that would
+     convert to them) raise a `TypeError` naming the right tool; `np.angle`
+     raises for every mp-complex input since it branches on dtype and no subclass
+     property can reach it; everything else passes through untouched;
+   - the raw `.real`/`.imag` attributes on a plain self-built ndarray remain the
+     single lying spelling — documented with a warning and a pinning test.
+
+   (En route, fixed a copy-paste bug: the scalar `mp.imag` was bound to
+   `boost::multiprecision::real` and returned the real part.)
+
+7. **Loops write output slots only through `slot_write`** (added 2026-07-08, after
+   CI segfaults on numpy 2.5.1): numpy may hand a loop an output slot that
+   **bitwise-aliases** another slot's mpfr allocation — numpy 2.5 initializes the
+   accumulator of an *identityless* reduce (`np.min`/`np.max`) by `memcpy` of
+   element 0, so the accumulator and `v[0]` share one set of limbs.  A plain
+   BMP assignment move-frees the slot's old limbs (freeing `v[0]`'s storage:
+   use-after-free → double-free → corrupted allocator → SIGSEGV several calls
+   later — the CI crash landed in `np.median`, three tests after the damage)
+   and writes through shared storage (`np.max(v)` silently rewrote `v[0]`).
+   `slot_write` computes the value first, memsets the slot to BMP's
+   uninitialized sentinel, and move-assigns the fresh value in — no existing
+   allocation is ever freed or written through.  Same crash-into-bounded-leak
+   trade as `HardenSetitem`.  Diagnosed with valgrind (UAF pair between the
+   `maximum` and `minimum` reduce loops); regression test
+   `test_identityless_reduce_does_not_corrupt_input` runs the reduces
+   repeatedly and asserts the input array survives.  numpy < 2.5 never
+   aliased, which is why every pre-2.5 environment was green.
+
+## Consequences
+
+- The "Known gotchas" docs page shrinks to the real, permanent edges: the float64
+  boundary (by design), complex component access (numpy limitation), complex
+  ordering (by design), no mp→int casts.  The reductions section becomes a
+  historical note.
+- `python/test/classes/numpy_ufuncs_test.py` pins: element-wise agreement with the
+  scalar functions, the numpy-semantics corners above, precision preservation
+  through every loop shape (including `sign`/`reciprocal`, which cross the mixed
+  real/complex division path with the known boost precision-mis-tagging hazard —
+  loops re-tag via `at_precision_of`), unwritten-slot safety per loop shape,
+  sorting/arg-extrema (nan-wins semantics), reductions, and the dangling-scalar
+  regressions.
+- Not upstreamed to eigenpy (the guarded loops already diverge; see ADR-0006).
+  Upstreaming the owned-copy `getitem` would fix the aliasing class for all
+  eigenpy user types and may be worth an issue later.
+
+## Relation to prior ADRs
+
+- **ADR-0006** — the slot-guard doctrine these loops follow; its "known gaps" list
+  shrinks (compare/argmax slots now filled and guarded).
+- **ADR-0031 / #259** — root-caused and fixed at the binding level by the owned-copy
+  `getitem`; the Python-side copy rule is now belt-and-suspenders.
+- **ADR-0001/0008** — unrelated eigenpy hazards, unchanged.

docs/adr/README.md

@@ -67,3 +67,4 @@ Each ADR follows the template:
 | [0047](0047-casual-records-surface.md) | The casual records surface: bertini.solve/save/load, Solution = points that remember, CLI records-on-by-default beside b1 files (no flag) | Python + CLI / records |
 | [0048](0048-cauchy-endgame-security-and-operating-zone.md) | Cauchy divergence handling: security check watches the ENDPOINT, truncates only in the operating zone, no pole-growth truncation (the acceptance gate alone cures junk-success); restores cyclic-6's 156 solutions (refines #70) | Core / endgames |
 | [0050](0050-docs-deploy-from-branch-not-deploy-pages.md) | Docs publish by serving the docs-store branch directly (Pages branch-source), NOT actions/deploy-pages — which failed structurally with BlobNotFound on the versioned store; custom domain via a /CNAME file; /stable/ is a redirect, .doctrees dropped | CI / docs |
+| [0051](0051-numpy-ufunc-coverage-and-owned-getitem.md) | Full numpy ufunc coverage for the mp dtypes (abs/conj/transcendentals/min-max/rounding/predicates, guarded loops), sort/argmax dtype slots (real only), getitem returns owned copies (kills the ADR-0031/#259 aliasing class; reductions genuinely safe); float64 boundary stays closed | Python bindings / numpy |

python/bertini/__init__.py

@@ -122,9 +122,6 @@
 from . import _slice_ops as _slice_ops
 _slice_ops.install(nag_algorithm.Slice)
 
-from . import operators                          # `from bertini.operators import *` -> just the math ops
-
-
 # --- sympy interop (#295) ----------------------------------------------------------------------
 # Make every function-tree node auto-convert to sympy (the `_sympy_` protocol), so sympy.sympify(node),
 # sympy.Matrix(array_of_nodes), and sympy.det(J) work directly.  sympy is an optional dependency; the
@@ -146,21 +143,39 @@ def __getattr__(name):
     raise AttributeError("module {!r} has no attribute {!r}".format(__name__, name))
 
 
-# --- numpy-compatible elementwise helpers for the mp types (#298, #301) ------------------------
-# Vectorized real/imag/abs/conj/round/sum/norm/is_real that stay mp-native (numpy's ufuncs/attrs on
-# the custom dtypes are unreliable -- see docs/source/known_gotchas.rst).  These are attributes of the
+# --- numpy-friendly elementwise helpers for the mp types (#298, #301) --------------------------
+# Vectorized real/imag/abs/conj/round/sum/norm/is_real that stay mp-native, riding the native numpy
+# ufunc loops on mp-dtype arrays and working element-wise on lists/mixed input (see the
+# "Multiprecision numbers and NumPy" docs page, docs/source/numpy.rst).  These are attributes of the
 # top-level module (bertini.abs, bertini.real, ...).  The builtin-shadowing names (abs, round, sum) are
 # deliberately kept OUT of __all__, so `from bertini import *` never clobbers the Python builtins.
 from . import _numpy_helpers as _numpy_helpers
+
+# make np.real/np.imag/np.angle raise (instead of silently returning wrong values)
+# on plain mp-complex arrays -- numpy has no user-dtype hook for component access,
+# and a crash is better than incorrect values.  See _numpy_guard for the story.
+from . import _numpy_guard as _numpy_guard
+_numpy_guard.install()
+
 real = _numpy_helpers.real
 imag = _numpy_helpers.imag
 conj = _numpy_helpers.conj
+arg = _numpy_helpers.arg
 norm = _numpy_helpers.norm
 is_real = _numpy_helpers.is_real
 abs = _numpy_helpers.abs        # noqa: A001  (bertini.abs; not exported via *)
 round = _numpy_helpers.round    # noqa: A001
 sum = _numpy_helpers.sum        # noqa: A001
 
+# --- the one-stop math vocabulary ---------------------------------------------------------------
+# `from bertini.operators import *` gives sin/cos/.../abs/arg/real/imag/... that work on symbolic
+# expressions AND numbers AND numpy containers alike, dispatching per argument.  The top-level
+# elementary functions are rebound to the polymorphic versions (a strict superset of the symbolic
+# ones bound above): bertini.sin(x) works for a Variable, a real_mp, or an array.
+from . import operators
+from .operators import (sin, cos, tan, asin, acos, atan, exp, log, sqrt,   # noqa: F811
+                        sinh, cosh, tanh, asinh, acosh, atanh)
+
 
 
 # https://stackoverflow.com/questions/44834/what-does-all-mean-in-python
@@ -170,7 +185,7 @@ def __getattr__(name):
            'jacobian','random_matrix','random_vector','random_real','random_complex','coefficient','coefficients',
            'complex_mp','real_mp','int_mp','rational_mp',
            'nag_algorithm','default_precision','is_distinct_up_to',
-           'real','imag','conj','norm','is_real',
+           'real','imag','conj','arg','norm','is_real',
            'tracking','endgame','logging','symbolics','parse','multiprec','random','parallel',
            'operators',
            # everyday classes hoisted to the top level
@@ -181,6 +196,7 @@ def __getattr__(name):
            # symbolic constants
            'E','Pi','I',
            'sin','cos','tan','asin','acos','atan','exp','log','sqrt',
+           'sinh','cosh','tanh','asinh','acosh','atanh',
            'canonicalize','monomial_order']