From 39899184975a56337c77ce8dfd70f08e40852823 Mon Sep 17 00:00:00 2001 From: silviana amethyst Date: Wed, 8 Jul 2026 15:38:56 +0000 Subject: [PATCH 1/8] feat(bindings): register the full numpy ufunc set for mp dtypes; sort/argmax slots np.abs, np.conj, the transcendental family (exp/log/trig/hyperbolic + inverses), power, sign, reciprocal, minimum/maximum/fmin/fmax, the rounding family (floor/ceil/trunc/rint half-to-even via mpfr_rint), remainder/fmod/floor_divide (numpy sign semantics), arctan2/hypot/ copysign, and the isnan/isinf/isfinite/signbit predicates now have guarded loops on the real_mp/complex_mp numpy dtypes -- previously "ufunc not supported" TypeErrors. All loops read through value_or_zero per the uninitialized-slot doctrine (ADR-0006) and call the same boost::multiprecision free functions the multiprec scalar functions bind. absolute on complex outputs real_mp. New HardenCompare/HardenArgMinMax fill the compare/argmax/argmin PyArray_ArrFuncs slots for real_mp: np.sort/argsort/searchsorted/ median/argmax/argmin now work (complex stays unordered by design). Also fixes a copy-paste bug where the scalar free function mp.imag() returned the REAL part, and adds array overloads of multiprec real/imag/arg: numpy's ndarray .real/.imag attributes silently return wrong values for legacy user dtypes (numpy cannot know complex_mp is complex-like), so these are the sanctioned array component accessors and the np.angle replacement. Co-Authored-By: Claude Fable 5 --- .../include/eigenpy_interaction.hpp | 451 ++++++++++++++++-- python_bindings/src/mpfr_export.cpp | 32 +- 2 files changed, 437 insertions(+), 46 deletions(-) diff --git a/python_bindings/include/eigenpy_interaction.hpp b/python_bindings/include/eigenpy_interaction.hpp index bcd0455b1..cf00b022f 100644 --- a/python_bindings/include/eigenpy_interaction.hpp +++ b/python_bindings/include/eigenpy_interaction.hpp @@ -8,6 +8,7 @@ #include "python_common.hpp" #include +#include #include #include @@ -182,6 +183,222 @@ namespace eigenpy } }; + // trait: is this scalar the complex mp type? several ops (conjugate, sign, + // the isnan/isinf/isfinite predicates) need a different body for complex. + template struct is_complex_mp : std::false_type {}; + template <> struct is_complex_mp : std::true_type {}; + + // re-tag `val` to carry `ref`'s precision. guards against boost's mixed + // real/complex arithmetic occasionally mis-tagging the result's precision + // (division is the known offender); .precision(n) preserves the value. + template + inline T at_precision_of(T val, T const& ref) + { + if (val.precision() != ref.precision()) + val.precision(ref.precision()); + return val; + } + + // ----- unary ops, same-type output ------------------------------------ + // bodies call the same boost::multiprecision free functions the multiprec + // module binds as scalar functions, so np.exp(arr)[i] == mp.exp(arr[i]). + + struct op_positive { template static T apply(T const& x) { return x; } }; + struct op_reciprocal + { + template static T apply(T const& x) + { + T one(1); + one.precision(x.precision()); + return at_precision_of(T(one / x), x); + } + }; + struct op_conjugate + { + template static T apply(T const& x) + { + if constexpr (is_complex_mp::value) + return T(conj(x)); + else + return x; + } + }; + // numpy-2 sign semantics: real -> -1/0/+1; complex -> z/|z| (0 at 0). + struct op_sign + { + template static T apply(T const& x) + { + if constexpr (is_complex_mp::value) + { + if (x == 0) + return at_precision_of(T(0), x); + bertini::real_mp const mag(abs(x)); + return at_precision_of(T(bertini::real_mp(x.real() / mag), + bertini::real_mp(x.imag() / mag)), x); + } + else + { + T res(x > 0 ? 1 : (x < 0 ? -1 : 0)); + return at_precision_of(std::move(res), x); + } + } + }; + + struct op_exp { template static T apply(T const& x) { return T(exp(x)); } }; + struct op_log { template static T apply(T const& x) { return T(log(x)); } }; + struct op_log10 { template static T apply(T const& x) { return T(log10(x)); } }; + struct op_exp2 { template static T apply(T const& x) { return T(exp2(x)); } }; + struct op_log2 { template static T apply(T const& x) { return T(log2(x)); } }; + struct op_expm1 { template static T apply(T const& x) { return T(expm1(x)); } }; + struct op_log1p { template static T apply(T const& x) { return T(log1p(x)); } }; + struct op_cbrt { template static T apply(T const& x) { return T(cbrt(x)); } }; + + struct op_sin { template static T apply(T const& x) { return T(sin(x)); } }; + struct op_cos { template static T apply(T const& x) { return T(cos(x)); } }; + struct op_tan { template static T apply(T const& x) { return T(tan(x)); } }; + struct op_arcsin { template static T apply(T const& x) { return T(asin(x)); } }; + struct op_arccos { template static T apply(T const& x) { return T(acos(x)); } }; + struct op_arctan { template static T apply(T const& x) { return T(atan(x)); } }; + struct op_sinh { template static T apply(T const& x) { return T(sinh(x)); } }; + struct op_cosh { template static T apply(T const& x) { return T(cosh(x)); } }; + struct op_tanh { template static T apply(T const& x) { return T(tanh(x)); } }; + struct op_arcsinh { template static T apply(T const& x) { return T(asinh(x)); } }; + struct op_arccosh { template static T apply(T const& x) { return T(acosh(x)); } }; + struct op_arctanh { template static T apply(T const& x) { return T(atanh(x)); } }; + + // real-only rounding family + struct op_floor { template static T apply(T const& x) { return T(floor(x)); } }; + struct op_ceil { template static T apply(T const& x) { return T(ceil(x)); } }; + struct op_trunc { template static T apply(T const& x) { return T(trunc(x)); } }; + // numpy's rint is round-half-to-EVEN; boost's rint rounds half away from + // zero, so call mpfr directly in MPFR_RNDN (nearest, ties to even). + struct op_rint + { + template static T apply(T const& x) + { + T out(0); + out.precision(x.precision()); + mpfr_rint(out.backend().data(), x.backend().data(), MPFR_RNDN); + return out; + } + }; + + // ----- unary ops, cross-type output ----------------------------------- + + // absolute: real -> real, complex -> real (the magnitude). + struct op_absolute + { + template static bertini::real_mp apply(T const& x) + { + return bertini::real_mp(abs(x)); + } + }; + struct op_fabs { template static T apply(T const& x) { return T(fabs(x)); } }; + + // predicates -> bool. the isnan/isinf/isfinite family are function-like + // macros in C , so call the boost versions qualified. + struct op_isnan + { + template static bool apply(T const& x) + { + if constexpr (is_complex_mp::value) + return boost::multiprecision::isnan(x.real()) || boost::multiprecision::isnan(x.imag()); + else + return boost::multiprecision::isnan(x); + } + }; + struct op_isinf + { + template static bool apply(T const& x) + { + if constexpr (is_complex_mp::value) + return boost::multiprecision::isinf(x.real()) || boost::multiprecision::isinf(x.imag()); + else + return boost::multiprecision::isinf(x); + } + }; + struct op_isfinite + { + template static bool apply(T const& x) + { + if constexpr (is_complex_mp::value) + return boost::multiprecision::isfinite(x.real()) && boost::multiprecision::isfinite(x.imag()); + else + return boost::multiprecision::isfinite(x); + } + }; + struct op_signbit + { + template static bool apply(T const& x) + { + return boost::multiprecision::signbit(x); + } + }; + + // ----- binary ops ------------------------------------------------------ + + struct op_power { template static T apply(T const& x, T const& y) { return T(pow(x, y)); } }; + struct op_arctan2 { template static T apply(T const& x, T const& y) { return T(atan2(x, y)); } }; + struct op_hypot { template static T apply(T const& x, T const& y) { return T(hypot(x, y)); } }; + struct op_copysign { template static T apply(T const& x, T const& y) { return T(copysign(x, y)); } }; + // numpy fmod keeps C fmod's sign-of-dividend semantics + struct op_fmod { template static T apply(T const& x, T const& y) { return T(fmod(x, y)); } }; + // numpy remainder/mod takes the sign of the DIVISOR (python % semantics) + struct op_remainder + { + template static T apply(T const& x, T const& y) + { + T r(fmod(x, y)); + if (r != 0 && ((r < 0) != (y < 0))) + r += y; + return r; + } + }; + struct op_floor_divide + { + template static T apply(T const& x, T const& y) + { + return T(floor(x / y)); + } + }; + // minimum/maximum propagate nan (numpy semantics); fmin/fmax ignore it + struct op_minimum + { + template static T apply(T const& x, T const& y) + { + if (boost::multiprecision::isnan(x)) return x; + if (boost::multiprecision::isnan(y)) return y; + return y < x ? y : x; + } + }; + struct op_maximum + { + template static T apply(T const& x, T const& y) + { + if (boost::multiprecision::isnan(x)) return x; + if (boost::multiprecision::isnan(y)) return y; + return x < y ? y : x; + } + }; + struct op_fmin + { + template static T apply(T const& x, T const& y) + { + if (boost::multiprecision::isnan(x)) return y; + if (boost::multiprecision::isnan(y)) return x; + return y < x ? y : x; + } + }; + struct op_fmax + { + template static T apply(T const& x, T const& y) + { + if (boost::multiprecision::isnan(x)) return y; + if (boost::multiprecision::isnan(y)) return x; + return x < y ? y : x; + } + }; + template void guarded_binary_op( char **args, EIGENPY_NPY_CONST_UFUNC_ARG npy_intp *dimensions, @@ -240,6 +457,28 @@ namespace eigenpy } } + // unary loop with an output type different from the input type + // (absolute: complex -> real; the isnan family: T -> bool). Writes into + // mp-typed output slots go through BMP operator=, which initializes a + // zeroed destination itself; bool slots are plain bytes. + template + void guarded_unary_op_out( + char **args, EIGENPY_NPY_CONST_UFUNC_ARG npy_intp *dimensions, + EIGENPY_NPY_CONST_UFUNC_ARG npy_intp *steps, void * /*data*/) + { + npy_intp is = steps[0], os = steps[1], n = *dimensions; + char *i = args[0], *o = args[1]; + const T zero(0); + for (npy_intp k = 0; k < n; ++k) + { + T const& x = value_or_zero(*reinterpret_cast(i), zero); + OutT& res = *reinterpret_cast(o); + res = Op::apply(x); + i += is; + o += os; + } + } + // guarded matmul: mirrors eigenpy::internal::{matrix_multiply,gufunc_matrix_multiply} // stride logic, with the inner dot product reading through value_or_zero. template @@ -326,6 +565,53 @@ namespace eigenpy *reinterpret_cast(op) = acc; } + // guarded element comparison for the PyArray_ArrFuncs `compare` slot + // (np.sort / argsort / searchsorted / unique). eigenpy leaves this slot + // empty for user dtypes ("type does not have compare function"). Only + // installed for the real type — complex has no ordering. nan compares + // false both ways (weak-ordering violation, same as C doubles): sorting + // arrays containing nan gives an unspecified nan position, not a crash. + template + int guarded_compare(const void *a, const void *b, void * /*arr*/) + { + const T zero(0); + T const& x = value_or_zero(*static_cast(a), zero); + T const& y = value_or_zero(*static_cast(b), zero); + if (x < y) return -1; + if (y < x) return 1; + return 0; + } + + // guarded argmax/argmin for the PyArray_ArrFuncs slots (np.argmax / + // np.argmin / np.max / np.min dispatch through these for user dtypes on + // some numpy paths). Mirrors numpy's float semantics: a nan wins + // immediately (first nan is the arg-extremum). numpy hands these a + // contiguous buffer. + template + int guarded_argminmax(void *data, npy_intp n, npy_intp *extremum_ind, void * /*arr*/) + { + const T zero(0); + T const* p = static_cast(data); + *extremum_ind = 0; + if (n == 0) + return 0; + T best = value_or_zero(p[0], zero); + if (boost::multiprecision::isnan(best)) + return 0; + for (npy_intp k = 1; k < n; ++k) + { + T const& v = value_or_zero(p[k], zero); + if (boost::multiprecision::isnan(v) || (Max ? best < v : v < best)) + { + *extremum_ind = k; + if (boost::multiprecision::isnan(v)) + return 0; + best = v; + } + } + return 0; + } + } // namespace internal @@ -378,6 +664,29 @@ namespace eigenpy funcs->dotfunc = reinterpret_cast(&internal::guarded_dotfunc); } + // Fill the element-comparison slot (empty in eigenpy's registration), enabling + // np.sort / np.argsort / np.searchsorted / np.unique. Real type only — + // complex has no ordering. Call immediately after eigenpy::registerNewType. + template + void HardenCompare() + { + PyArray_Descr *descr = Register::getPyArrayDescr(); + PyArray_ArrFuncs *funcs = PyDataType_GetArrFuncs(descr); + funcs->compare = reinterpret_cast(&internal::guarded_compare); + } + + // Fill the argmax/argmin slots (empty in eigenpy's registration), enabling + // np.argmax / np.argmin ("data type not ordered" otherwise). Real type only. + // Call immediately after eigenpy::registerNewType. + template + void HardenArgMinMax() + { + PyArray_Descr *descr = Register::getPyArrayDescr(); + PyArray_ArrFuncs *funcs = PyDataType_GetArrFuncs(descr); + funcs->argmax = reinterpret_cast(&internal::guarded_argminmax); + funcs->argmin = reinterpret_cast(&internal::guarded_argminmax); + } + // register a single guarded loop on the named numpy ufunc, mirroring the // error handling of eigenpy's EIGENPY_REGISTER_*_UFUNC macros. inline void registerGuardedLoop(PyObject *numpy, char const *ufunc_name, @@ -414,13 +723,20 @@ namespace eigenpy // i lifted this from EigenPy and adapted it: all loops are the guarded // versions from internal:: above (eigenpy's read input slots unguarded — - // see the header comment), and the ordering comparitors are a compile-time - // option because they are NOT defined for complex types (instantiating - // them for complex_mp would be a hard error). + // see the header comment), and the ordering-dependent set is a compile-time + // option because ordering is NOT defined for complex types (instantiating + // those functors for complex_mp would be a hard error). Coverage beyond + // eigenpy's arithmetic core (absolute, conjugate, the transcendental family, + // rounding, min/max, the isnan predicates) closes the documented + // "ufunc not supported" gotchas — every loop body calls the same + // boost::multiprecision free function the multiprec module binds as the + // scalar function of the same name. template void registerGuardedUfunct() { const int type_code = Register::getTypeCode(); + const int bool_code = Register::getTypeCode(); + const int real_code = Register::getTypeCode(); PyObject *numpy_str; #if PY_MAJOR_VERSION >= 3 @@ -434,6 +750,19 @@ namespace eigenpy import_ufunc(); + // registration helpers: (in...) -> out signatures. the types array is + // copied by PyUFunc_RegisterLoopForType, so stack storage is fine. + auto unary = [&](char const* name, PyUFuncGenericFunction loop, int out_code) + { + int types[2] = {type_code, out_code}; + registerGuardedLoop(numpy, name, type_code, loop, types, 2); + }; + auto binary = [&](char const* name, PyUFuncGenericFunction loop, int out_code) + { + int types[3] = {type_code, type_code, out_code}; + registerGuardedLoop(numpy, name, type_code, loop, types, 3); + }; + // Matrix multiply { int types[3] = {type_code, type_code, type_code}; @@ -442,49 +771,81 @@ namespace eigenpy types, 3); } - // Binary operators - { - int types[3] = {type_code, type_code, type_code}; - registerGuardedLoop(numpy, "add", type_code, - &internal::guarded_binary_op, types, 3); - registerGuardedLoop(numpy, "subtract", type_code, - &internal::guarded_binary_op, types, 3); - registerGuardedLoop(numpy, "multiply", type_code, - &internal::guarded_binary_op, types, 3); - registerGuardedLoop(numpy, "divide", type_code, - &internal::guarded_binary_op, types, 3); - } - - // Comparison operators - { - int types[3] = {type_code, type_code, Register::getTypeCode()}; - registerGuardedLoop(numpy, "equal", type_code, - &internal::guarded_compare_op, types, 3); - registerGuardedLoop(numpy, "not_equal", type_code, - &internal::guarded_compare_op, types, 3); - - if constexpr (WithOrderingComparitors) // NOT defined for complex types - { - registerGuardedLoop(numpy, "greater", type_code, - &internal::guarded_compare_op, types, 3); - registerGuardedLoop(numpy, "less", type_code, - &internal::guarded_compare_op, types, 3); - registerGuardedLoop(numpy, "greater_equal", type_code, - &internal::guarded_compare_op, types, 3); - registerGuardedLoop(numpy, "less_equal", type_code, - &internal::guarded_compare_op, types, 3); - } - } - - // Unary operators + // Binary arithmetic + binary("add", &internal::guarded_binary_op, type_code); + binary("subtract", &internal::guarded_binary_op, type_code); + binary("multiply", &internal::guarded_binary_op, type_code); + binary("divide", &internal::guarded_binary_op, type_code); + binary("power", &internal::guarded_binary_op, type_code); + + // Equality comparisons (defined for real and complex alike) + binary("equal", &internal::guarded_compare_op, bool_code); + binary("not_equal", &internal::guarded_compare_op, bool_code); + + // Unary, same-type output + unary("negative", &internal::guarded_unary_op, type_code); + unary("positive", &internal::guarded_unary_op, type_code); + unary("square", &internal::guarded_unary_op, type_code); + unary("sqrt", &internal::guarded_unary_op, type_code); + unary("reciprocal", &internal::guarded_unary_op, type_code); + unary("conjugate", &internal::guarded_unary_op, type_code); + unary("sign", &internal::guarded_unary_op, type_code); + unary("exp", &internal::guarded_unary_op, type_code); + unary("log", &internal::guarded_unary_op, type_code); + unary("log10", &internal::guarded_unary_op, type_code); + unary("sin", &internal::guarded_unary_op, type_code); + unary("cos", &internal::guarded_unary_op, type_code); + unary("tan", &internal::guarded_unary_op, type_code); + unary("arcsin", &internal::guarded_unary_op, type_code); + unary("arccos", &internal::guarded_unary_op, type_code); + unary("arctan", &internal::guarded_unary_op, type_code); + unary("sinh", &internal::guarded_unary_op, type_code); + unary("cosh", &internal::guarded_unary_op, type_code); + unary("tanh", &internal::guarded_unary_op, type_code); + unary("arcsinh", &internal::guarded_unary_op, type_code); + unary("arccosh", &internal::guarded_unary_op, type_code); + unary("arctanh", &internal::guarded_unary_op, type_code); + + // absolute: real -> real, complex -> real (magnitude) + unary("absolute", &internal::guarded_unary_op_out, real_code); + + // predicates -> bool + unary("isnan", &internal::guarded_unary_op_out, bool_code); + unary("isinf", &internal::guarded_unary_op_out, bool_code); + unary("isfinite", &internal::guarded_unary_op_out, bool_code); + + if constexpr (WithOrderingComparitors) // the ordering-dependent set; NOT defined for complex types { - int types[2] = {type_code, type_code}; - registerGuardedLoop(numpy, "negative", type_code, - &internal::guarded_unary_op, types, 2); - registerGuardedLoop(numpy, "square", type_code, - &internal::guarded_unary_op, types, 2); - registerGuardedLoop(numpy, "sqrt", type_code, - &internal::guarded_unary_op, types, 2); + binary("greater", &internal::guarded_compare_op, bool_code); + binary("less", &internal::guarded_compare_op, bool_code); + binary("greater_equal", &internal::guarded_compare_op, bool_code); + binary("less_equal", &internal::guarded_compare_op, bool_code); + + // real-only unary: rounding family, fabs, real-only transcendentals + unary("floor", &internal::guarded_unary_op, type_code); + unary("ceil", &internal::guarded_unary_op, type_code); + unary("trunc", &internal::guarded_unary_op, type_code); + unary("rint", &internal::guarded_unary_op, type_code); + unary("fabs", &internal::guarded_unary_op, type_code); + unary("exp2", &internal::guarded_unary_op, type_code); + unary("log2", &internal::guarded_unary_op, type_code); + unary("expm1", &internal::guarded_unary_op, type_code); + unary("log1p", &internal::guarded_unary_op, type_code); + unary("cbrt", &internal::guarded_unary_op, type_code); + + unary("signbit", &internal::guarded_unary_op_out, bool_code); + + // real-only binary + binary("arctan2", &internal::guarded_binary_op, type_code); + binary("hypot", &internal::guarded_binary_op, type_code); + binary("copysign", &internal::guarded_binary_op, type_code); + binary("fmod", &internal::guarded_binary_op, type_code); + binary("remainder", &internal::guarded_binary_op, type_code); + binary("floor_divide", &internal::guarded_binary_op, type_code); + binary("minimum", &internal::guarded_binary_op, type_code); + binary("maximum", &internal::guarded_binary_op, type_code); + binary("fmin", &internal::guarded_binary_op, type_code); + binary("fmax", &internal::guarded_binary_op, type_code); } Py_DECREF(numpy); diff --git a/python_bindings/src/mpfr_export.cpp b/python_bindings/src/mpfr_export.cpp index d6ed10e65..ef3447ce5 100644 --- a/python_bindings/src/mpfr_export.cpp +++ b/python_bindings/src/mpfr_export.cpp @@ -259,10 +259,30 @@ namespace bertini{ using boost::multiprecision::imag; real_mp (*reeeal)(const T&) = &boost::multiprecision::real; - real_mp (*imaaag)(const T&) = &boost::multiprecision::real; + // regression note: this was `&boost::multiprecision::real` (copy-paste), + // so mp.imag(z) returned the REAL part. + real_mp (*imaaag)(const T&) = &boost::multiprecision::imag; def("real",reeeal, (arg("val")), "get the real part"); //,return_value_policy() def("imag",imaaag, (arg("val")), "get the imaginary part"); //,return_value_policy() + // array overloads: numpy's .real/.imag ndarray attributes (and np.real / + // np.imag / np.angle) return silently WRONG values for legacy user + // dtypes -- numpy does not know complex_mp is complex-like, so .real + // returns the complex values themselves and .imag returns zeros. These + // element-wise overloads are the sanctioned array component accessors. + def("real", +[](Vec const& v) { + Vec r(v.size()); + for (Eigen::Index i = 0; i < v.size(); ++i) + r(i) = v(i).real(); + return r; + }, (arg("val")), "get the real parts of an array of complex numbers, as an array of real_mp"); + def("imag", +[](Vec const& v) { + Vec r(v.size()); + for (Eigen::Index i = 0; i < v.size(); ++i) + r(i) = v(i).imag(); + return r; + }, (arg("val")), "get the imaginary parts of an array of complex numbers, as an array of real_mp"); + // and then a few more free functions // def("abs2",&T::abs2); @@ -276,6 +296,14 @@ namespace bertini{ real_mp (*aaaarg)(const T&) = &boost::multiprecision::arg; def("arg",aaaarg, "the argument, or the angle from 0. beware the branch cut."); + // array overload: the np.angle replacement (np.angle goes through the + // broken .imag/.real ndarray attributes -- see the note at real/imag). + def("arg", +[](Vec const& v) { + Vec r(v.size()); + for (Eigen::Index i = 0; i < v.size(); ++i) + r(i) = boost::multiprecision::arg(v(i)); + return r; + }, (arg("val")), "the arguments (angles from 0) of an array of complex numbers, as an array of real_mp. beware the branch cut."); // def("square",&square); // def("cube",&cube); @@ -448,6 +476,8 @@ namespace bertini{ eigenpy::registerNewType(); eigenpy::HardenSetitem(); // zero slots before assignment — see eigenpy_interaction.hpp & ADR-0003 eigenpy::HardenDotfunc(); // np.dot/np.inner guard — see eigenpy_interaction.hpp + eigenpy::HardenCompare(); // element compare slot: np.sort/argsort/searchsorted + eigenpy::HardenArgMinMax(); // argmax/argmin slots // guarded loops (real type — orderings included); eigenpy's registerCommonUfunc // loops read input slots unguarded and crash on never-written np.zeros/np.empty slots. eigenpy::registerGuardedUfunct(); From b122956e48185ce73eaa6906a4fd64b071784080 Mon Sep 17 00:00:00 2001 From: silviana amethyst Date: Wed, 8 Jul 2026 15:47:34 +0000 Subject: [PATCH 2/8] fix(bindings): getitem returns an owned copy, not a boost::ref into the buffer Stock eigenpy's getitem (and our heal-on-read specializations, which copied it) returns boost::ref(slot) -- a Python scalar ALIASING numpy array storage. Scalars extracted from temporary arrays dangled once the array was freed: np.sum(v) kept past the statement read freed memory (zeros, garbage precision, MPFR assertion SIGABRT on str()), and np.mean returned a silently-wrong 0 through the same mechanism. This is the root of the ADR-0031 / #259 hazard class, previously worked around consumer-by-consumer with copy-at-extraction. Returning an owned copy kills the class at the source: an indexed element is a durable value, as numpy users expect, and the identity-seeded reductions (np.sum/np.prod/np.mean) are now genuinely safe rather than accidentally readable. Adds python/test/classes/numpy_ufuncs_test.py: element-wise agreement of every new ufunc with the multiprec scalar functions, numpy semantics corners (rint half-to-even, remainder sign-of-divisor, nan propagation in minimum/maximum vs fmin/fmax), sort/argmax/searchsorted/ median, reductions regression, precision preservation through every loop shape, unwritten-slot safety, the mp.imag regression, the array real/imag/arg accessors, and named regressions for the dangling-scalar fix. Co-Authored-By: Claude Fable 5 --- python/test/classes/numpy_ufuncs_test.py | 436 ++++++++++++++++++ .../include/eigenpy_interaction.hpp | 17 +- 2 files changed, 449 insertions(+), 4 deletions(-) create mode 100644 python/test/classes/numpy_ufuncs_test.py diff --git a/python/test/classes/numpy_ufuncs_test.py b/python/test/classes/numpy_ufuncs_test.py new file mode 100644 index 000000000..18dad87fd --- /dev/null +++ b/python/test/classes/numpy_ufuncs_test.py @@ -0,0 +1,436 @@ +# This file is part of Bertini 2. +# +# python/test/classes/numpy_ufuncs_test.py is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# python/test/classes/numpy_ufuncs_test.py is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with python/test/classes/numpy_ufuncs_test.py. If not, see . +# +# Copyright(C) Bertini2 Development Team +# +# See for a copy of the license, +# as well as COPYING. Bertini2 is provided with permitted +# additional terms in the b2/licenses/ directory. + +# individual authors of this file include: +# +# silviana amethyst +# summer 2026 + +""" +Tests for the numpy ufunc coverage of the multiprecision dtypes. + +Every registered loop calls the same boost::multiprecision free function that +the ``bertini.multiprec`` scalar function of the same name binds, so a ufunc +applied to an array must agree element-wise with the scalar function -- +exactly, not approximately, since both run identical code at identical +precision. + +Also pins: +- the identity-seeded reductions (np.sum/np.prod/np.mean) that older docs + declared unsupported (they work on current numpy; regression-guard them), +- the sort/argsort/searchsorted/argmax family (dtype compare/arg slots), +- numpy semantics corners: rint half-to-even, remainder sign-of-divisor, + nan propagation in minimum/maximum vs fmin/fmax, +- precision preservation through every loop shape, +- unwritten-slot safety (np.empty) per the ADR-0006 guard doctrine, +- the mp.imag copy-paste regression (returned the real part), and the array + overloads of multiprec real/imag/arg that replace the silently-wrong + ndarray .real/.imag attributes and np.angle. +""" + +import numpy as np +import pytest + +import bertini.multiprec as mp +from bertini.multiprec import complex_mp, real_mp + + +@pytest.fixture(params=[real_mp, complex_mp], ids=["real_mp", "complex_mp"]) +def dtype(request): + return request.param + + +def _sample(dtype): + """A small array of values safe for every domain-restricted function.""" + if dtype is real_mp: + return np.array([real_mp('0.25'), real_mp('0.5'), real_mp('0.75')]) + return np.array([complex_mp('0.25', '0.125'), + complex_mp('0.5', '-0.25'), + complex_mp('-0.75', '0.375')]) + + +# ufuncs defined for both dtypes, with the matching multiprec scalar function +UFUNCS_BOTH = [ + (np.exp, mp.exp), (np.log, mp.log), (np.sqrt, mp.sqrt), + (np.sin, mp.sin), (np.cos, mp.cos), (np.tan, mp.tan), + (np.arcsin, mp.asin), (np.arccos, mp.acos), (np.arctan, mp.atan), + (np.sinh, mp.sinh), (np.cosh, mp.cosh), (np.tanh, mp.tanh), + (np.arcsinh, mp.asinh), (np.arccosh, None), (np.arctanh, mp.atanh), +] + + +class TestElementwiseAgreesWithScalarFunctions: + """np.f(arr)[i] must equal mp.f(arr[i]) exactly (identical code path).""" + + @pytest.mark.parametrize( + "ufunc,scalar", UFUNCS_BOTH, + ids=[u.__name__ for u, _ in UFUNCS_BOTH]) + def test_transcendental(self, dtype, ufunc, scalar): + if ufunc is np.arccosh: + # acosh needs |x| >= 1 on the real line + v = (np.array([real_mp('1.5'), real_mp(2)]) if dtype is real_mp + else np.array([complex_mp('1.5', '0.5')])) + scalar = mp.acosh + else: + v = _sample(dtype) + out = ufunc(v) + assert out.dtype == np.dtype(dtype) + for got, x in zip(out, v): + assert got == scalar(x) + + def test_absolute(self, dtype): + v = _sample(dtype) + out = np.abs(v) + # output is ALWAYS real, also for complex input (the magnitude) + assert out.dtype == np.dtype(real_mp) + for got, x in zip(out, v): + assert got == mp.abs(x) + + def test_conjugate(self, dtype): + v = _sample(dtype) + out = np.conj(v) + assert out.dtype == np.dtype(dtype) + for got, x in zip(out, v): + assert got == (mp.conj(x) if dtype is complex_mp else x) + + def test_power(self, dtype): + v = _sample(dtype) + out = np.power(v, v) + for got, x in zip(out, v): + assert got == x ** x + + def test_reciprocal(self, dtype): + v = _sample(dtype) + out = np.reciprocal(v) + for got, x in zip(out, v): + assert got == dtype(1) / x + + def test_square_negative_positive(self, dtype): + v = _sample(dtype) + assert all(np.square(v)[i] == v[i] * v[i] for i in range(len(v))) + assert all(np.negative(v)[i] == -v[i] for i in range(len(v))) + assert all(np.positive(v)[i] == v[i] for i in range(len(v))) + + def test_real_only_transcendentals(self): + v = _sample(real_mp) + pairs = [(np.log10, mp.log), (np.exp2, None), (np.log2, None), + (np.expm1, None), (np.log1p, None), (np.cbrt, None)] + # spot-check values against the double versions loosely; the exact + # contract (same boost call as a scalar) has no bound scalar twin for + # these, so compare against float64 at double precision. + for ufunc, _ in pairs: + got = ufunc(v) + want = ufunc(np.array([float(x) for x in v])) + for g, w in zip(got, want): + assert abs(float(g) - w) < 1e-14, ufunc.__name__ + + def test_sign_real(self): + v = np.array([real_mp(-3), real_mp(0), real_mp(2)]) + assert [str(s) for s in np.sign(v)] == ['-1', '0', '1'] + + def test_sign_complex_is_unit_modulus(self): + # numpy-2 semantics: sign(z) = z/|z|, 0 at 0 + z = complex_mp(3, 4) + s = np.sign(np.array([z, complex_mp(0)])) + assert s[0] == complex_mp('0.6', '0.8') + assert s[1] == complex_mp(0) + + def test_arctan2_hypot_copysign(self): + # atan2/hypot use dedicated algorithms, so agreement with the composed + # formulas is to the last ulp, not bit-exact + tol = real_mp('1e-25') + a = np.array([real_mp(1), real_mp(-2)]) + b = np.array([real_mp(3), real_mp(4)]) + assert mp.abs(np.arctan2(a, b)[0] - mp.atan(a[0] / b[0])) < tol + assert mp.abs(np.hypot(a, b)[1] - mp.sqrt(a[1] * a[1] + b[1] * b[1])) < tol + got = np.copysign(a, np.array([real_mp(-1), real_mp(1)])) + assert [str(x) for x in got] == ['-1', '2'] + + +class TestNumpySemanticsCorners: + """Corners where numpy's semantics differ from naive C/boost calls.""" + + def test_rint_rounds_half_to_even(self): + # regression: boost's rint rounds half AWAY from zero; numpy (and the + # loop, via mpfr_rint in MPFR_RNDN) rounds half to even + v = np.array([real_mp('0.5'), real_mp('1.5'), real_mp('2.5'), + real_mp('-0.5'), real_mp('-2.5')]) + assert [str(x) for x in np.rint(v)] == ['0', '2', '2', '-0', '-2'] + + def test_floor_ceil_trunc(self): + v = np.array([real_mp('1.7'), real_mp('-1.7')]) + assert [str(x) for x in np.floor(v)] == ['1', '-2'] + assert [str(x) for x in np.ceil(v)] == ['2', '-1'] + assert [str(x) for x in np.trunc(v)] == ['1', '-1'] + + def test_remainder_takes_sign_of_divisor(self): + a = np.array([real_mp(7), real_mp(-7), real_mp(7), real_mp(-7)]) + b = np.array([real_mp(3), real_mp(3), real_mp(-3), real_mp(-3)]) + assert [str(x) for x in np.mod(a, b)] == ['1', '2', '-2', '-1'] + # fmod keeps C semantics (sign of dividend) + assert [str(x) for x in np.fmod(a, b)] == ['1', '-1', '1', '-1'] + assert [str(x) for x in np.floor_divide(a, b)] == ['2', '-3', '-3', '2'] + + def test_minimum_maximum_propagate_nan_fmin_fmax_ignore_it(self): + nan = real_mp('nan') + one = np.array([real_mp(1)]) + nans = np.array([nan]) + assert mp.abs(np.fmax(nans, one)[0] - real_mp(1)) == 0 + assert mp.abs(np.fmin(nans, one)[0] - real_mp(1)) == 0 + assert np.isnan(np.maximum(nans, one))[0] + assert np.isnan(np.minimum(nans, one))[0] + + def test_minimum_maximum_values(self): + a = np.array([real_mp(1), real_mp(5)]) + b = np.array([real_mp(3), real_mp(2)]) + assert [str(x) for x in np.minimum(a, b)] == ['1', '2'] + assert [str(x) for x in np.maximum(a, b)] == ['3', '5'] + + def test_predicates(self, dtype): + good = _sample(dtype) + assert not np.isnan(good).any() + assert np.isfinite(good).all() + assert not np.isinf(good).any() + assert np.isnan(good).dtype == np.dtype(bool) + if dtype is real_mp: + bad = np.array([real_mp('nan'), real_mp('inf'), real_mp(1)]) + else: + bad = np.array([complex_mp('nan', '0'), complex_mp('0', 'inf'), + complex_mp(1, 1)]) + assert list(np.isnan(bad)) == [True, False, False] + assert list(np.isinf(bad)) == [False, True, False] + assert list(np.isfinite(bad)) == [False, False, True] + + def test_signbit(self): + v = np.array([real_mp(-2), real_mp(0), real_mp(3)]) + assert list(np.signbit(v)) == [True, False, False] + + +class TestSortingAndArgExtrema: + """The dtype compare/argmax/argmin slots (real only -- complex is unordered).""" + + def test_sort_and_argsort(self): + v = np.array([real_mp(3), real_mp(1), real_mp(2)]) + assert [str(x) for x in np.sort(v)] == ['1', '2', '3'] + assert list(np.argsort(v)) == [1, 2, 0] + + def test_argmax_argmin_max_min(self): + v = np.array([real_mp(3), real_mp(1), real_mp(7), real_mp(2)]) + assert np.argmax(v) == 2 + assert np.argmin(v) == 1 + assert str(np.max(v)) == '7' + assert str(np.min(v)) == '1' + + def test_argmax_nan_wins(self): + # numpy float semantics: the first nan is the arg-extremum + v = np.array([real_mp(1), real_mp('nan'), real_mp(3)]) + assert np.argmax(v) == 1 + assert np.argmin(v) == 1 + + def test_searchsorted_and_median(self): + v = np.array([real_mp(1), real_mp(2), real_mp(4)]) + assert np.searchsorted(v, real_mp(3)) == 2 + assert str(np.median(v)) == '2' + + def test_complex_stays_unordered(self): + w = np.array([complex_mp(1, 2), complex_mp(0, 1)]) + with pytest.raises(TypeError): + np.sort(w) + + +class TestReductions: + """Identity-seeded reductions -- previously documented as crashing. + + They work on current numpy; these tests exist so any numpy/eigenpy + combination that breaks them again fails loudly here instead of in + user code. + """ + + def test_sum_prod_mean_real(self): + v = np.array([real_mp(1), real_mp(2), real_mp(3)]) + assert np.sum(v) == real_mp(6) + assert np.prod(v) == real_mp(6) + assert np.mean(v) == real_mp(2) + + def test_sum_prod_mean_complex(self): + w = np.array([complex_mp(1, 2), complex_mp(3, 4)]) + assert np.sum(w) == complex_mp(4, 6) + assert np.prod(w) == complex_mp(-5, 10) + assert np.mean(w) == complex_mp(2, 3) + + def test_bare_reduce_and_cumsum(self, dtype): + v = np.array([dtype(1), dtype(2), dtype(3)]) + assert np.add.reduce(v) == dtype(6) + assert list(np.cumsum(v)) == [dtype(1), dtype(3), dtype(6)] + + def test_reduce_with_explicit_initial_still_works(self, dtype): + # the old workaround from the gotchas page must keep working + v = np.array([dtype(1), dtype(2)]) + assert np.add.reduce(v, initial=dtype(10)) == dtype(13) + + +class TestCloseness: + """np.isclose/np.allclose need float64 promotion, which stays disabled by + design (double -> mp casts are deliberately unsafe: the user has to think + about float-literal intent). This is the sanctioned all-mp idiom.""" + + def test_allclose_idiom(self, dtype): + v = _sample(dtype) + w = v + dtype('1e-20') + assert np.all(np.abs(v - w) <= real_mp('1e-10')) + assert not np.all(np.abs(v - (w + dtype(1))) <= real_mp('1e-10')) + + def test_isclose_itself_still_raises(self, dtype): + # if this ever starts passing, numpy grew user-dtype promotion -- + # revisit the known-gotchas page + v = _sample(dtype) + with pytest.raises(TypeError): + np.isclose(v, v) + + +class TestPrecisionPreservation: + """Outputs carry the operands' precision, not the ambient default -- + including through the mixed real/complex division inside sign and + reciprocal (the known boost precision-mis-tagging hazard).""" + + HIGH = 50 + + def _high_precision_sample(self, dtype): + mp.default_precision(self.HIGH) + v = (np.array([real_mp('1.5')]) if dtype is real_mp + else np.array([complex_mp('1.5', '2.5')])) + mp.default_precision(30) + assert v[0].precision == self.HIGH + return v + + @pytest.mark.parametrize("ufunc", [np.exp, np.sqrt, np.abs, np.sign, + np.reciprocal, np.conj, np.rint], + ids=lambda u: u.__name__) + def test_unary_output_precision(self, dtype, ufunc): + if dtype is complex_mp and ufunc is np.rint: + pytest.skip("rint is real-only") + v = self._high_precision_sample(dtype) + assert ufunc(v)[0].precision == self.HIGH + + def test_binary_output_precision(self, dtype): + v = self._high_precision_sample(dtype) + assert np.power(v, v)[0].precision == self.HIGH + + +class TestUnwrittenSlotSafety: + """Every new loop shape must survive never-written np.empty slots + (which hold the all-zero BMP sentinel) -- the ADR-0006 doctrine.""" + + def test_unary_loops_on_empty(self, dtype): + e = np.empty(3, dtype=dtype) + for ufunc in (np.exp, np.sin, np.conj, np.sign, np.abs, + np.isnan, np.isfinite): + ufunc(e) # must not crash + + def test_binary_loops_on_empty(self, dtype): + e = np.empty(3, dtype=dtype) + np.power(e, e) + if dtype is real_mp: + np.minimum(e, e) + np.arctan2(e, e) + np.mod(e, e) + + def test_sort_and_argmax_on_empty(self): + e = np.empty(4, dtype=real_mp) + np.sort(e) + np.argmax(e) + np.argmin(e) + + def test_reductions_on_zeros(self, dtype): + z = np.zeros(3, dtype=dtype) + assert np.sum(z) == dtype(0) + + +class TestScalarsAreOwnedCopies: + """Regression tests for the getitem aliasing fix. + + getitem used to return boost::ref into the numpy buffer (as stock eigenpy + does), so a scalar extracted from a temporary array -- most visibly the + result of np.sum/np.mean -- dangled once the array was freed: reading it + later gave zeros/garbage or SIGABRT inside mpfr (the ADR-0031 / #259 + hazard class). getitem now returns an owned copy. + """ + + def test_reduce_scalar_survives_its_array(self, dtype): + s = np.sum(np.array([dtype(1), dtype(2), dtype(3)])) + # the source (temporary) array is gone; s must still be intact + assert s == dtype(6) + assert str(s) is not None # printing used to MPFR-assert on the corpse + assert s / dtype(3) == dtype(2) + + def test_mean_is_correct_not_silently_zero(self, dtype): + # np.mean's internal divide ran on a dangling extraction and returned 0 + v = np.array([dtype(1), dtype(2), dtype(3)]) + assert np.mean(v) == dtype(2) + + def test_stored_indexed_elements_stay_distinct(self): + # the ADR-0031 shape, now safe at the binding level (still copy in + # Python code by convention) + pts = np.array([complex_mp(1, 1), complex_mp(2, 2), complex_mp(3, 3)]) + kept = [pts[i] for i in range(3)] + del pts + assert [str(k) for k in kept] == ['(1,1)', '(2,2)', '(3,3)'] + + def test_mutating_an_extracted_scalar_leaves_the_array_alone(self): + v = np.array([real_mp(1), real_mp(2)]) + x = v[0] + x += real_mp(10) + assert v[0] == real_mp(1) + + +class TestComponentAccessors: + """The multiprec real/imag/arg array overloads, and the scalar imag + regression.""" + + def test_scalar_imag_returns_imaginary_part(self): + # regression: mp.imag was bound to boost::multiprecision::real by a + # copy-paste error, so it returned the REAL part + z = complex_mp(1, 2) + assert mp.real(z) == real_mp(1) + assert mp.imag(z) == real_mp(2) + + def test_array_real_imag(self): + w = np.array([complex_mp(1, 2), complex_mp(3, 4)]) + r, i = mp.real(w), mp.imag(w) + assert r.dtype == np.dtype(real_mp) and i.dtype == np.dtype(real_mp) + assert [str(x) for x in r] == ['1', '3'] + assert [str(x) for x in i] == ['2', '4'] + + def test_array_arg_replaces_np_angle(self): + w = np.array([complex_mp(1, 1), complex_mp(-1, 0)]) + a = mp.arg(w) + assert a.dtype == np.dtype(real_mp) + assert a[0] == mp.arg(w[0]) + assert a[1] == mp.arg(w[1]) + + def test_ndarray_real_imag_attributes_are_untrustworthy(self): + # documenting-by-test: numpy cannot know a legacy user dtype is + # complex-like, so ndarray .real returns the complex values themselves + # and .imag returns zeros. If numpy ever fixes this, the accessors + # above become optional and the gotchas page should be updated. + w = np.array([complex_mp(1, 2)]) + assert w.real.dtype == np.dtype(complex_mp) # not real_mp! + assert w.imag[0] == complex_mp(0) # wrong value, by numpy diff --git a/python_bindings/include/eigenpy_interaction.hpp b/python_bindings/include/eigenpy_interaction.hpp index cf00b022f..b50b8eba6 100644 --- a/python_bindings/include/eigenpy_interaction.hpp +++ b/python_bindings/include/eigenpy_interaction.hpp @@ -77,7 +77,15 @@ namespace eigenpy } - // template specialization for real numbers + // template specialization for real numbers. + // + // NB: eigenpy's stock getitem (and an earlier version of this one) returns + // boost::ref(slot) — a Python object ALIASING the numpy buffer. That is + // the root of the ADR-0031 hazard family (#259: stored elements silently + // collapse or SIGABRT once the buffer is reused/freed), and it makes + // scalars extracted from temporary arrays — np.sum/np.mean results, most + // visibly — dangle outright. Returning an owned COPY kills the whole + // class: an indexed element is a durable value, as numpy users expect. template <> struct getitem { @@ -91,13 +99,14 @@ namespace eigenpy { mpfr_scalar = NumT(0); } - boost::python::object m(boost::ref(mpfr_scalar)); + boost::python::object m(mpfr_scalar); // owned copy — never boost::ref (see above) Py_INCREF(m.ptr()); return m.ptr(); } }; - // a template specialization for complex numbers + // a template specialization for complex numbers; see the real one for the + // copy-not-ref rationale. template <> struct getitem { @@ -111,7 +120,7 @@ namespace eigenpy { mpfr_scalar = NumT(0); } - boost::python::object m(boost::ref(mpfr_scalar)); + boost::python::object m(mpfr_scalar); // owned copy — never boost::ref (see above) Py_INCREF(m.ptr()); return m.ptr(); } From 4b97cf82b2862e2820865a2a0e59bbf8b3cf3dc1 Mon Sep 17 00:00:00 2001 From: silviana amethyst Date: Wed, 8 Jul 2026 15:55:42 +0000 Subject: [PATCH 3/8] docs: rewrite known-gotchas for the new numpy coverage; ADR-0051 The gotchas page now documents what IS supported (the full ufunc set, sorting, reductions -- verified on numpy 2.3/2.4 and regression-tested) and shrinks the gotchas to the real, permanent edges: the float64 boundary (by design -- double->mp casts stay unsafe, reaffirmed 2026-07-08: the user has to think about float-literal intent; the sanctioned closeness idiom is np.all(np.abs(a-b) <= real_mp('1e-10'))), complex component access (.real/.imag/np.angle silently wrong for legacy user dtypes -- use the new mp.real/imag/arg array overloads), complex ordering (deliberately unimplemented), and no mp->int casts. The reductions section becomes a version-qualified historical note with the initial= idiom kept as the old-numpy fallback. ADR-0051 records the whole decision set, including the owned-copy getitem that root-causes the ADR-0031/#259 aliasing class. Also updates the stale endgame_test comment, the multiprec docstring (and its pre-rename type names), and a stale dtype name in the precision-models tutorial. All 180 sphinx doctests pass. Co-Authored-By: Claude Fable 5 --- ...-numpy-ufunc-coverage-and-owned-getitem.md | 117 +++++++++++ docs/adr/README.md | 1 + python/bertini/multiprec/__init__.py | 16 +- python/docs/source/known_gotchas.rst | 189 +++++++++++++----- .../precision_models/index.rst | 4 +- python/test/tracking/endgame_test.py | 9 +- 6 files changed, 273 insertions(+), 63 deletions(-) create mode 100644 docs/adr/0051-numpy-ufunc-coverage-and-owned-getitem.md diff --git a/docs/adr/0051-numpy-ufunc-coverage-and-owned-getitem.md b/docs/adr/0051-numpy-ufunc-coverage-and-owned-getitem.md new file mode 100644 index 000000000..ca95383b5 --- /dev/null +++ b/docs/adr/0051-numpy-ufunc-coverage-and-owned-getitem.md @@ -0,0 +1,117 @@ +# 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** (reaffirmed 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. + Visible consequences, documented as gotchas rather than fixed: + `np.isclose`/`np.allclose` raise `DTypePromotionError` (sanctioned idiom: + `np.all(np.abs(a - b) <= real_mp('1e-10'))`), and `np.round(a, decimals≠0)` + fails (it scales by a float internally). + +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, so the ndarray `.real`/`.imag` attributes + (and `np.real`/`np.imag`/`np.angle`) return silently wrong values — `.real` + gives the complex values, `.imag` gives zeros. Not hookable from the bindings; + documented, plus 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.) + +## 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. diff --git a/docs/adr/README.md b/docs/adr/README.md index c50fe2370..878340bc5 100644 --- a/docs/adr/README.md +++ b/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 | diff --git a/python/bertini/multiprec/__init__.py b/python/bertini/multiprec/__init__.py index 4504c670b..eaf463fbf 100644 --- a/python/bertini/multiprec/__init__.py +++ b/python/bertini/multiprec/__init__.py @@ -33,12 +33,12 @@ """ Multiprecision types, and functions that operate on them. -Numeric types exposed are +Numeric types exposed are -* Complex (Boost.Multiprecision mpc) -* Float (Boost.Multiprecision mpfr) -* Int (Boost.Multiprecision mpz) -* Rational (Boost.Multiprecision.mpq) +* complex_mp (Boost.Multiprecision mpc) +* real_mp (Boost.Multiprecision mpfr) +* int_mp (Boost.Multiprecision mpz) +* rational_mp (Boost.Multiprecision.mpq) This namespace also includes the mathematical operators, like `cos`, etc. """ @@ -48,7 +48,11 @@ from bertini._pybertini.multiprec import * # (no Vector helper: eigenpy makes the mp number types work as numpy dtypes directly, so a plain -# numpy array -- e.g. np.zeros(n, dtype=bertini.complex_mp) -- is the vector.) +# numpy array -- e.g. np.zeros(n, dtype=bertini.complex_mp) -- is the vector. numpy ufuncs +# (np.abs, np.exp, np.sum, ...) work on such arrays; the exceptions are documented on the +# "Known gotchas" docs page. For the real/imaginary parts or argument of a COMPLEX ARRAY use +# this module's real()/imag()/arg() -- the ndarray .real/.imag attributes and np.angle return +# silently wrong values for user-defined dtypes, a numpy limitation.) __all__ = dir(_pybmp) diff --git a/python/docs/source/known_gotchas.rst b/python/docs/source/known_gotchas.rst index ce13f3acb..20c870ffe 100644 --- a/python/docs/source/known_gotchas.rst +++ b/python/docs/source/known_gotchas.rst @@ -4,78 +4,167 @@ A few sharp edges fall out of how Bertini 2's multiprecision numbers are exposed to NumPy. They are collected here with the idiom that works and the idiom that bites. -.. _gotcha-numpy-reductions: +Bertini 2 exposes :class:`~bertini.real_mp` (variable-precision real) and +:class:`~bertini.complex_mp` (variable-precision complex) as **custom NumPy dtypes** +(via eigenpy), and nearly everything works on arrays of them: + +* element-wise arithmetic, comparisons, ``@`` / :func:`numpy.dot` / + :func:`numpy.linalg.norm`, +* the whole ufunc family: :func:`numpy.abs`, :func:`numpy.conj`, ``exp`` / ``log`` / + the trigonometric and hyperbolic functions and their inverses, ``power``, ``sign``, + ``sqrt``, the rounding family (``floor`` / ``ceil`` / ``trunc`` / ``rint``), + ``minimum`` / ``maximum``, ``isnan`` / ``isinf`` / ``isfinite``, and friends, +* sorting and order statistics on the real type: :func:`numpy.sort`, + :func:`numpy.argsort`, :func:`numpy.searchsorted`, :func:`numpy.argmax` / + :func:`numpy.argmin`, :func:`numpy.median`, +* the identity-seeded reductions: :func:`numpy.sum`, :func:`numpy.prod`, + :func:`numpy.mean`, :func:`numpy.cumsum`, bare ``ufunc.reduce``. + +Every ufunc loop calls the same multiprecision function that the scalar +:mod:`bertini.multiprec` function of the same name binds, at the operands' precision -- +``np.exp(a)[i]`` is exactly ``mp.exp(a[i])``. + +The sharp edges that remain are below. They are *by design* (the float64 boundary) or +*by NumPy limitation* (component access on a user dtype), not bugs to be waited out. + +.. _gotcha-float64-boundary: + +Mixing float64 into multiprecision arrays +========================================== + +A Python ``float`` (or NumPy ``float64``) does **not** silently promote into a +multiprecision array. This is deliberate: ``0.1`` as a float64 is not the number you +typed -- it carries binary noise past the 16th digit -- so anywhere a double could +sneak into a high-precision computation, Bertini makes you say what you mean +(construct from a *string*: ``real_mp('0.1')``). + +The visible consequences: + +* ``np.isclose(a, b)`` / ``np.allclose(a, b)`` raise ``DTypePromotionError`` on mp + arrays, because their float tolerances (``rtol=1e-05``, ``atol=1e-08``) cannot + promote. +* ``np.round(a, decimals)`` with nonzero ``decimals`` fails for the same reason (it + scales by a float power of ten internally). Plain ``np.round(a)`` / + :func:`numpy.rint` work. +* arithmetic between an mp array and a float scalar/array fails with + ``ufunc ... not supported``. + +Use all-multiprecision operands instead. The closeness idiom: -NumPy reductions over multiprecision arrays -============================================ +.. doctest:: -Bertini 2 exposes :class:`~bertini.real_mp` (variable-precision real) and -:class:`~bertini.complex_mp` (variable-precision complex) as **custom NumPy -dtypes** (via eigenpy). Element-wise math, indexing, ``@`` / :func:`numpy.dot`, -:func:`numpy.linalg.norm`, :func:`numpy.cumsum` and friends all work on arrays of these -dtypes. + >>> import numpy as np + >>> from bertini.multiprec import real_mp, complex_mp + >>> a = np.array([complex_mp(3), complex_mp(4)]) + >>> b = np.array([complex_mp(3), complex_mp(4)]) + + >>> # ✓ the all-mp replacement for np.allclose(a, b) + >>> bool(np.all(np.abs(a - b) <= real_mp('1e-10'))) + True + + >>> # ✗ np.allclose itself cannot work: its float64 tolerances cannot promote + >>> np.allclose(a, b) + Traceback (most recent call last): + ... + numpy.exceptions.DTypePromotionError: The DType could not be promoted by . This means that no common DType exists for the given inputs. For example they cannot be stored in a single array unless the dtype is `object`. The full list of DTypes is: (, ) -The sharp edge is the **identity-seeded reductions** -- :func:`numpy.sum`, -:func:`numpy.prod`, :func:`numpy.mean`, and a bare ``ufunc.reduce`` with no ``initial=``: +Integers are fine in both directions of intent -- they are exact, so they convert and +promote freely (``a * 3``, ``np.power(a, 2)``, ``real_mp(7)``). And converting mp +*down* to double is available when you ask for it explicitly (``float(x)``, +``complex(z)``, ``arr.astype(complex)``) -- you are consciously truncating. -.. code-block:: python +.. _gotcha-complex-components: + +Component access on complex arrays: never ``.real`` / ``.imag`` +================================================================ + +NumPy does not know a user-defined dtype is complex-like, so the ndarray attributes +``.real`` and ``.imag`` (and :func:`numpy.real`, :func:`numpy.imag`, +:func:`numpy.angle`, which route through them) return **silently wrong values** on +``complex_mp`` arrays: ``.real`` returns the complex values themselves and ``.imag`` +returns zeros. This is a NumPy limitation for legacy user dtypes; Bertini cannot hook +those attributes. + +Use the :mod:`bertini.multiprec` component functions, which accept arrays: + +.. doctest:: >>> import numpy as np + >>> import bertini.multiprec as mp >>> from bertini.multiprec import complex_mp - >>> v = np.array([complex_mp(3), complex_mp(4)], dtype=complex_mp) - >>> np.sum(v) # ✗ DON'T -- may crash - Traceback (most recent call last): - ... - SystemError: returned NULL without setting an exception + >>> w = np.array([complex_mp(1, 2), complex_mp(3, 4)]) + + >>> # ✓ arrays of the real/imaginary parts, as real_mp + >>> [float(x) for x in mp.real(w)] + [1.0, 3.0] + >>> [float(x) for x in mp.imag(w)] + [2.0, 4.0] + + >>> # ✓ the np.angle replacement + >>> [round(float(x), 4) for x in mp.arg(w)] + [1.1071, 0.9273] + + >>> # ✗ the ndarray attribute is wrong (numpy returns zeros -- not Bertini's doing) + >>> complex(w.imag[0]) + 0j -Why this happens ----------------- +On complex *scalars* the properties are correct (``w[0].imag`` is exact); it is only +the *array* attributes that lie. -To reduce an array, NumPy needs the reduction's **identity element** -- ``0`` for a sum, -``1`` for a product -- materialized *in the array's dtype*. For these custom dtypes, -NumPy (2.x) tries to build that identity from a Python ``int`` and fails *inside its own -machinery*, before ever calling Bertini's dtype code, returning ``NULL`` without setting a -Python exception (hence the bare ``SystemError``). It is a limitation of the NumPy ↔ -eigenpy boundary for legacy user-defined dtypes, **not** something Bertini can patch in its -bindings, and whether it triggers depends on the exact NumPy / eigenpy / Eigen build -- so -it may "work" on one machine and crash on another. Treat it as always-unsupported. +.. _gotcha-complex-ordering: -What to do instead ------------------- +Complex is unordered +===================== -Give the reduction a value to start from, or use an operation that seeds itself from the -data (``dot`` / ``norm`` / a plain Python loop). All of these are stable across builds: +There is no ``<`` on complex numbers, so ``complex_mp`` arrays do not support +:func:`numpy.sort`, :func:`numpy.argmax`, :func:`numpy.minimum` / ``maximum``, or the +ordering comparisons. (NumPy's built-in ``complex128`` sorts lexicographically for +historical reasons; Bertini deliberately does not reproduce that.) Sort a derived real +quantity instead -- e.g. ``np.argsort(np.abs(w))``. + +Casts that do not exist +======================== + +There is deliberately **no** cast from the multiprecision types to any integer type. +``arr.astype(int)`` on an mp array will fail; go through double first if you truly +want it, accepting the truncation. + +.. _gotcha-numpy-reductions: + +A historical note on reductions +================================ + +Earlier versions of this page declared :func:`numpy.sum` / :func:`numpy.prod` / +:func:`numpy.mean` unsupported on mp arrays: on some older NumPy builds the +identity-seeded reduce failed inside NumPy with ``SystemError: ... returned NULL +without setting an exception``. With current NumPy (verified on 2.3 and 2.4 series; +regression-tested in CI on all platforms) they work: .. doctest:: >>> import numpy as np >>> from bertini.multiprec import real_mp, complex_mp - >>> v = np.array([complex_mp(3), complex_mp(4)], dtype=complex_mp) - >>> w = np.array([real_mp(1), real_mp(2), real_mp(3)], dtype=real_mp) + >>> v = np.array([complex_mp(3), complex_mp(4)]) + >>> w = np.array([real_mp(1), real_mp(2), real_mp(3)]) - >>> # ✓ sum: hand ufunc.reduce an explicit, correctly-typed identity - >>> complex(np.add.reduce(v, initial=complex_mp(0))) + >>> complex(np.sum(v)) (7+0j) - >>> float(np.add.reduce(w, initial=real_mp(0))) - 6.0 + >>> float(np.mean(w)) + 2.0 + +If you are pinned to an old NumPy and see that ``SystemError``, the old idioms all +still work and remain the portable fallback: - >>> # ✓ sum: or just use Python's built-in sum() +.. doctest:: + + >>> # explicit, correctly-typed identity + >>> complex(np.add.reduce(v, initial=complex_mp(0))) + (7+0j) + >>> # python's own sum >>> float(sum(w)) 6.0 - - >>> # ✓ Euclidean norm works directly (it routes through the dtype's dot slot) + >>> # norm / sum of squares route through the dtype's dot slot >>> float(np.linalg.norm(v)) 5.0 - - >>> # ✓ sum of squares: np.dot does not conjugate, so dot(v, v) == sum(v**2) >>> complex(np.dot(v, v)) (25+0j) - - >>> # ✓ mean: reduce with an identity, then divide by the count - >>> float(np.add.reduce(w, initial=real_mp(0)) / w.size) - 2.0 - -In short: anywhere you would reach for ``np.sum(a)`` / ``np.prod(a)`` / ``np.mean(a)`` on a -``real_mp`` or ``complex_mp`` array, reach for ``np.add.reduce(a, initial=...)`` (or -``np.multiply.reduce(a, initial=...)``), ``np.dot`` / ``np.linalg.norm``, or a plain Python -``sum`` instead. diff --git a/python/docs/source/tutorials/settings_and_precision/precision_models/index.rst b/python/docs/source/tutorials/settings_and_precision/precision_models/index.rst index 9bb6bff2c..5feb89a33 100644 --- a/python/docs/source/tutorials/settings_and_precision/precision_models/index.rst +++ b/python/docs/source/tutorials/settings_and_precision/precision_models/index.rst @@ -41,8 +41,8 @@ Reading solutions: convert with ``complex()`` The model changes the **type** of the numbers you get back. A double solve returns NumPy ``complex128``; a multiprecision solve returns arrays of :class:`bertini.complex_mp` (NumPy -arrays with ``dtype`` ``Complex``). The portable habit -- works in every model -- is to convert each -coordinate with :func:`complex`: +arrays with ``dtype`` ``complex_mp``). The portable habit -- works in every model -- is to convert +each coordinate with :func:`complex`: .. testcode:: diff --git a/python/test/tracking/endgame_test.py b/python/test/tracking/endgame_test.py index 2b2d9c917..7f423ecda 100644 --- a/python/test/tracking/endgame_test.py +++ b/python/test/tracking/endgame_test.py @@ -137,9 +137,8 @@ def test_using_total_degree_ss(): for soln in dehomogenized_solns: diff = exact_soln - soln - # NB: np.sum / np.prod / np.mean over multiprecision (mpfr/mpc) arrays can raise - # SystemError on some numpy + eigenpy builds -- numpy cannot construct the reduction - # identity element for these custom dtypes. See the "Known gotchas" page in the docs. - # np.dot(diff, diff) == sum(diff_i**2) (numpy's dot does not conjugate) and goes through - # the dtype's dot slot, which works everywhere; it preserves this assertion exactly. + # np.dot(diff, diff) == sum(diff_i**2) (numpy's dot does not conjugate); it goes + # through the dtype's dot slot. np.sum would work too on current numpy (see the + # reductions note on the "Known gotchas" docs page), but dot stays portable to + # older numpy builds where the identity-seeded reduce raised SystemError. assert mp.abs(np.sqrt(np.dot(diff, diff))) < 1e-10 From 0505a7d04e9a5c40580576749752de4a9a778972 Mon Sep 17 00:00:00 2001 From: silviana amethyst Date: Wed, 8 Jul 2026 18:38:24 +0000 Subject: [PATCH 4/8] feat(bindings): tolerance orderings vs float64; astype(complex) down-conversion Comparisons against a double tolerance -- np.abs(a - b) < 1e-10 -- are safe: an ordering yields a bool, so no float value ever flows into a multiprecision computation, and the comparison is exact (boost compares the number against the double directly; 1e-22 is not lost against a 1e-30 tolerance). This matches the C++ solvers' double ToleranceT and the line the scalar bindings already drew (GreatLessVisitor bound, double equality deliberately not). Register mixed mp-vs-float64 ordering loops (< <= > >=, both operand orders, real type only). Mixed EQUALITY with a float stays unregistered (exact equality against a float literal is the 0.1-intent trap the unsafe double->mp cast exists to block), as does mixed arithmetic; np.isclose/np.allclose still raise (they COMPUTE with float64 tolerances internally). The float boundary protects polynomial-system construction, not tolerance checks. Also register mp -> complex128 casts as unsafe, alongside the pre-existing mp -> double: arr.astype(complex) / astype(float) are the explicit, conscious truncations (astype(complex) had never actually been registered). Registration-order fix: casts now register 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 "registered/modified ... after the cast had been used"). Docs: the gotchas page float64-boundary section now states the values-vs-comparisons line and the plain-float closeness idiom; ADR-0051 decision 5 amended. Tests cover both operand orders, exactness of the mixed compare, the still-blocked equality/arithmetic/ isclose, and astype. Co-Authored-By: Claude Fable 5 --- ...-numpy-ufunc-coverage-and-owned-getitem.md | 28 +++-- python/docs/source/known_gotchas.rst | 37 ++++--- python/test/classes/numpy_ufuncs_test.py | 62 +++++++++-- .../include/eigenpy_interaction.hpp | 104 +++++++++++++++++- python_bindings/src/mpfr_export.cpp | 22 +++- 5 files changed, 214 insertions(+), 39 deletions(-) diff --git a/docs/adr/0051-numpy-ufunc-coverage-and-owned-getitem.md b/docs/adr/0051-numpy-ufunc-coverage-and-owned-getitem.md index ca95383b5..f86369029 100644 --- a/docs/adr/0051-numpy-ufunc-coverage-and-owned-getitem.md +++ b/docs/adr/0051-numpy-ufunc-coverage-and-owned-getitem.md @@ -73,14 +73,26 @@ that behavior faithfully. Consequences: `SystemError` and keep the `initial=` idiom as the fallback for older numpy. `pyproject.toml` keeps `numpy` unpinned. -5. **The float64 boundary stays closed** (reaffirmed 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. - Visible consequences, documented as gotchas rather than fixed: - `np.isclose`/`np.allclose` raise `DTypePromotionError` (sanctioned idiom: - `np.all(np.abs(a - b) <= real_mp('1e-10'))`), and `np.round(a, decimals≠0)` - fails (it scales by a float internally). +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` + 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 diff --git a/python/docs/source/known_gotchas.rst b/python/docs/source/known_gotchas.rst index 20c870ffe..52b78d18f 100644 --- a/python/docs/source/known_gotchas.rst +++ b/python/docs/source/known_gotchas.rst @@ -33,23 +33,29 @@ Mixing float64 into multiprecision arrays ========================================== A Python ``float`` (or NumPy ``float64``) does **not** silently promote into a -multiprecision array. This is deliberate: ``0.1`` as a float64 is not the number you +multiprecision *value*. This is deliberate: ``0.1`` as a float64 is not the number you typed -- it carries binary noise past the 16th digit -- so anywhere a double could sneak into a high-precision computation, Bertini makes you say what you mean (construct from a *string*: ``real_mp('0.1')``). -The visible consequences: - -* ``np.isclose(a, b)`` / ``np.allclose(a, b)`` raise ``DTypePromotionError`` on mp - arrays, because their float tolerances (``rtol=1e-05``, ``atol=1e-08``) cannot - promote. +The line is drawn at whether the float's *value* flows into the computation: + +* **Tolerance comparisons against a float are allowed.** ``np.abs(a - b) < 1e-10`` + works: an ordering comparison yields a ``bool``, so no float ever enters a + multiprecision value -- and the comparison is *exact* (the mp value is compared + against the double, not rounded down to double first). This matches the C++ + solvers, whose tolerances are doubles. +* **Mixed arithmetic is blocked** (``a + 0.1`` raises ``ufunc ... not supported``), + and so is **mixed equality** (``a == 0.1`` -- exact equality against a float + literal is precisely the trap the boundary exists for). +* ``np.isclose(a, b)`` / ``np.allclose(a, b)`` raise ``DTypePromotionError``: they + *compute* ``atol + rtol*np.abs(b)`` with float64 values internally, which is + arithmetic, not comparison. * ``np.round(a, decimals)`` with nonzero ``decimals`` fails for the same reason (it scales by a float power of ten internally). Plain ``np.round(a)`` / :func:`numpy.rint` work. -* arithmetic between an mp array and a float scalar/array fails with - ``ufunc ... not supported``. -Use all-multiprecision operands instead. The closeness idiom: +The closeness idiom is therefore just: .. doctest:: @@ -58,11 +64,15 @@ Use all-multiprecision operands instead. The closeness idiom: >>> a = np.array([complex_mp(3), complex_mp(4)]) >>> b = np.array([complex_mp(3), complex_mp(4)]) - >>> # ✓ the all-mp replacement for np.allclose(a, b) - >>> bool(np.all(np.abs(a - b) <= real_mp('1e-10'))) + >>> # ✓ the replacement for np.allclose(a, b): compare against a double tolerance + >>> bool(np.all(np.abs(a - b) < 1e-10)) True - >>> # ✗ np.allclose itself cannot work: its float64 tolerances cannot promote + >>> # ✓ the comparison is exact -- 1e-22 is not lost against a 1e-30 tolerance + >>> bool(np.all(np.array([real_mp('1e-22')]) < 1e-30)) + False + + >>> # ✗ np.allclose itself cannot work: it computes with float64 tolerances >>> np.allclose(a, b) Traceback (most recent call last): ... @@ -71,7 +81,8 @@ Use all-multiprecision operands instead. The closeness idiom: Integers are fine in both directions of intent -- they are exact, so they convert and promote freely (``a * 3``, ``np.power(a, 2)``, ``real_mp(7)``). And converting mp *down* to double is available when you ask for it explicitly (``float(x)``, -``complex(z)``, ``arr.astype(complex)``) -- you are consciously truncating. +``complex(z)``, ``arr.astype(float)``, ``arr.astype(complex)``) -- you are consciously +truncating. .. _gotcha-complex-components: diff --git a/python/test/classes/numpy_ufuncs_test.py b/python/test/classes/numpy_ufuncs_test.py index 18dad87fd..96edaf675 100644 --- a/python/test/classes/numpy_ufuncs_test.py +++ b/python/test/classes/numpy_ufuncs_test.py @@ -288,24 +288,70 @@ def test_reduce_with_explicit_initial_still_works(self, dtype): class TestCloseness: - """np.isclose/np.allclose need float64 promotion, which stays disabled by - design (double -> mp casts are deliberately unsafe: the user has to think - about float-literal intent). This is the sanctioned all-mp idiom.""" - - def test_allclose_idiom(self, dtype): + """The float64 boundary: tolerance ORDERINGS against a double are allowed + (a comparison yields a bool -- no float flows into a multiprecision value; + this mirrors the C++ solvers' double ToleranceT and the scalar + GreatLessVisitor). Everything that would let a float VALUE into + an mp computation stays closed: mixed equality, mixed arithmetic, + np.isclose's internal float tolerances.""" + + def test_tolerance_comparison_with_float(self, dtype): + v = _sample(dtype) + w = v + dtype('1e-20') + assert np.all(np.abs(v - w) <= 1e-10) + assert not np.all(np.abs(v - (w + dtype(1))) <= 1e-10) + # both operand orders, and float64 arrays as well as scalars + assert np.all(1e-10 >= np.abs(v - w)) + assert np.all(np.abs(v - w) < np.full(len(v), 1e-10)) + + def test_tolerance_comparison_is_exact_not_sloppy(self): + # the double is compared exactly (boost mixed compare), not by rounding + # the mp value down to double first + tiny = real_mp('1e-22') + assert np.all(np.array([tiny]) < 1e-10) + assert not np.any(np.array([tiny]) < 1e-30) + + def test_allclose_idiom_all_mp_still_works(self, dtype): v = _sample(dtype) w = v + dtype('1e-20') assert np.all(np.abs(v - w) <= real_mp('1e-10')) - assert not np.all(np.abs(v - (w + dtype(1))) <= real_mp('1e-10')) + + def test_float_equality_stays_blocked(self, dtype): + # exact equality against a float literal is the 0.1-intent trap; it is + # not bound at the scalar level either + v = _sample(dtype) + with pytest.raises(TypeError): + v == 0.1 + + def test_float_arithmetic_stays_blocked(self, dtype): + v = _sample(dtype) + with pytest.raises(TypeError): + v + 0.1 def test_isclose_itself_still_raises(self, dtype): - # if this ever starts passing, numpy grew user-dtype promotion -- - # revisit the known-gotchas page + # its internal float64 rtol/atol cannot promote; if this ever starts + # passing, numpy grew user-dtype promotion -- revisit the gotchas page v = _sample(dtype) with pytest.raises(TypeError): np.isclose(v, v) +class TestExplicitDownConversion: + """astype is the explicit, conscious truncation to double precision.""" + + def test_astype_float_and_complex(self): + v = np.array([real_mp('1.5'), real_mp(2)]) + w = np.array([complex_mp(1, 2), complex_mp(3, 4)]) + assert list(v.astype(float)) == [1.5, 2.0] + assert list(v.astype(complex)) == [1.5 + 0j, 2.0 + 0j] + assert list(w.astype(complex)) == [1 + 2j, 3 + 4j] + + def test_astype_int_stays_forbidden(self, dtype): + v = np.array([dtype(1)]) + with pytest.raises(TypeError): + v.astype(np.int64) + + class TestPrecisionPreservation: """Outputs carry the operands' precision, not the ambient default -- including through the mixed real/complex division inside sign and diff --git a/python_bindings/include/eigenpy_interaction.hpp b/python_bindings/include/eigenpy_interaction.hpp index b50b8eba6..b8ff5c022 100644 --- a/python_bindings/include/eigenpy_interaction.hpp +++ b/python_bindings/include/eigenpy_interaction.hpp @@ -7,6 +7,7 @@ #include "python_common.hpp" +#include #include #include @@ -173,12 +174,15 @@ namespace eigenpy struct op_subtract { template static T apply(T const& x, T const& y) { return T(x - y); } }; struct op_multiply { template static T apply(T const& x, T const& y) { return T(x * y); } }; struct op_divide { template static T apply(T const& x, T const& y) { return T(x / y); } }; - struct op_equal { template static bool apply(T const& x, T const& y) { return x == y; } }; - struct op_not_equal { template static bool apply(T const& x, T const& y) { return x != y; } }; - struct op_greater { template static bool apply(T const& x, T const& y) { return x > y; } }; - struct op_less { template static bool apply(T const& x, T const& y) { return x < y; } }; - struct op_greater_equal { template static bool apply(T const& x, T const& y) { return x >= y; } }; - struct op_less_equal { template static bool apply(T const& x, T const& y) { return x <= y; } }; + // comparison functors are heterogeneous (two type parameters) so the same + // functor serves the mp-vs-mp loops AND the mp-vs-double tolerance loops + // (boost::multiprecision compares a number against a double exactly). + struct op_equal { template static bool apply(T const& x, U const& y) { return x == y; } }; + struct op_not_equal { template static bool apply(T const& x, U const& y) { return x != y; } }; + struct op_greater { template static bool apply(T const& x, U const& y) { return x > y; } }; + struct op_less { template static bool apply(T const& x, U const& y) { return x < y; } }; + struct op_greater_equal { template static bool apply(T const& x, U const& y) { return x >= y; } }; + struct op_less_equal { template static bool apply(T const& x, U const& y) { return x <= y; } }; struct op_negative { template static T apply(T const& x) { return T(-x); } }; struct op_square { template static T apply(T const& x) { return T(x * x); } }; @@ -448,6 +452,44 @@ namespace eigenpy } } + // mixed mp-vs-float64 ORDERING comparison (Reversed swaps operand order: + // false = (mp, double), true = (double, mp)). Comparisons against a + // double tolerance -- np.abs(a - b) < 1e-10 -- are safe: boost compares a + // number against a double exactly, and the result is a bool, so no float + // ever flows INTO a multiprecision value. This mirrors the scalar + // bindings (GreatLessVisitor) and the C++ solvers' double + // ToleranceT. Deliberately orderings-only: mixed EQUALITY with a float + // literal is the 0.1-intent trap the unsafe double->mp cast exists to + // block, and it is not bound at the scalar level either. + template + void guarded_mixed_compare_op( + char **args, EIGENPY_NPY_CONST_UFUNC_ARG npy_intp *dimensions, + EIGENPY_NPY_CONST_UFUNC_ARG npy_intp *steps, void * /*data*/) + { + npy_intp is0 = steps[0], is1 = steps[1], os = steps[2], n = *dimensions; + char *i0 = args[0], *i1 = args[1], *o = args[2]; + const T zero(0); + for (npy_intp k = 0; k < n; ++k) + { + bool& res = *reinterpret_cast(o); + if constexpr (Reversed) + { + double const& x = *reinterpret_cast(i0); + T const& y = value_or_zero(*reinterpret_cast(i1), zero); + res = Op::apply(x, y); + } + else + { + T const& x = value_or_zero(*reinterpret_cast(i0), zero); + double const& y = *reinterpret_cast(i1); + res = Op::apply(x, y); + } + i0 += is0; + i1 += is1; + o += os; + } + } + template void guarded_unary_op( char **args, EIGENPY_NPY_CONST_UFUNC_ARG npy_intp *dimensions, @@ -649,6 +691,32 @@ namespace eigenpy } }; + // mp -> complex128, for the explicit down-conversion arr.astype(complex) + // (registered unsafe, like mp -> double: you consciously truncate). + // boost mp numbers have no conversion operator to std::complex, so go + // through the components. + template <> + struct cast> + { + static std::complex run(bertini::real_mp const& from) + { + if (internal::mpfr_slot::uninitialized(from)) + return {0.0, 0.0}; + return {from.convert_to(), 0.0}; + } + }; + + template <> + struct cast> + { + static std::complex run(bertini::complex_mp const& from) + { + if (internal::mpfr_slot::uninitialized(from)) + return {0.0, 0.0}; + return {from.real().convert_to(), from.imag().convert_to()}; + } + }; + // Install the zero-initialization setitem guard for an MPFR-backed dtype. // Call immediately after eigenpy::registerNewType(), before any arrays @@ -830,6 +898,30 @@ namespace eigenpy binary("greater_equal", &internal::guarded_compare_op, bool_code); binary("less_equal", &internal::guarded_compare_op, bool_code); + // mixed mp-vs-float64 orderings, both operand orders: the tolerance + // idiom `np.abs(a - b) < 1e-10`. Orderings ONLY -- see + // guarded_mixed_compare_op for why equality stays mp-vs-mp. + auto mixed_ordering = [&](char const* name, PyUFuncGenericFunction fwd, + PyUFuncGenericFunction rev) + { + int types_td[3] = {type_code, NPY_DOUBLE, bool_code}; + int types_dt[3] = {NPY_DOUBLE, type_code, bool_code}; + registerGuardedLoop(numpy, name, type_code, fwd, types_td, 3); + registerGuardedLoop(numpy, name, type_code, rev, types_dt, 3); + }; + mixed_ordering("greater", + &internal::guarded_mixed_compare_op, + &internal::guarded_mixed_compare_op); + mixed_ordering("less", + &internal::guarded_mixed_compare_op, + &internal::guarded_mixed_compare_op); + mixed_ordering("greater_equal", + &internal::guarded_mixed_compare_op, + &internal::guarded_mixed_compare_op); + mixed_ordering("less_equal", + &internal::guarded_mixed_compare_op, + &internal::guarded_mixed_compare_op); + // real-only unary: rounding family, fabs, real-only transcendentals unary("floor", &internal::guarded_unary_op, type_code); unary("ceil", &internal::guarded_unary_op, type_code); diff --git a/python_bindings/src/mpfr_export.cpp b/python_bindings/src/mpfr_export.cpp index ef3447ce5..484448c5e 100644 --- a/python_bindings/src/mpfr_export.cpp +++ b/python_bindings/src/mpfr_export.cpp @@ -478,9 +478,12 @@ namespace bertini{ eigenpy::HardenDotfunc(); // np.dot/np.inner guard — see eigenpy_interaction.hpp eigenpy::HardenCompare(); // element compare slot: np.sort/argsort/searchsorted eigenpy::HardenArgMinMax(); // argmax/argmin slots - // guarded loops (real type — orderings included); eigenpy's registerCommonUfunc - // loops read input slots unguarded and crash on never-written np.zeros/np.empty slots. - eigenpy::registerGuardedUfunct(); + + // casts must be registered BEFORE the ufunc loops: registering the + // mixed mp-vs-double ordering loops makes numpy query the mp<->double + // casts, and a cast first registered after it has been queried is + // ignored (numpy RuntimeWarning "registered/modified ... after the + // cast had been used"). // you can convert from integer types with no fear eigenpy::registerCast(true); @@ -497,6 +500,13 @@ namespace bertini{ // both directions are scary. eigenpy::registerCast(false); eigenpy::registerCast(false); + // explicit down-conversion arr.astype(complex) — unsafe, you consciously truncate + eigenpy::registerCast>(false); + + // guarded loops (real type — orderings included, plus the mp-vs-double + // tolerance orderings); eigenpy's registerCommonUfunc loops read input + // slots unguarded and crash on never-written np.zeros/np.empty slots. + eigenpy::registerGuardedUfunct(); IMPLICITLY_CONVERTIBLE(int,T); @@ -573,8 +583,8 @@ namespace bertini{ eigenpy::registerNewType(); eigenpy::HardenSetitem(); // zero slots before assignment — see eigenpy_interaction.hpp & ADR-0003 eigenpy::HardenDotfunc(); // np.dot/np.inner guard — see eigenpy_interaction.hpp - eigenpy::registerUfunct_without_comparitors(); + // casts before ufunc loops — see the note in ExposeFloat. // you can safely convert from integer types to Complex's, there's no loss possible eigenpy::registerCast(true); @@ -590,10 +600,14 @@ namespace bertini{ // these conversions are unsafe. you can, but you probably shouldn't eigenpy::registerCast(false); eigenpy::registerCast(false); + // explicit down-conversion arr.astype(complex) — unsafe, you consciously truncate + eigenpy::registerCast>(false); // it's ok to convert from variable precision Float to Complex, that's ok! eigenpy::registerCast(true); + eigenpy::registerUfunct_without_comparitors(); + IMPLICITLY_CONVERTIBLE(int,T); IMPLICITLY_CONVERTIBLE(long,T); IMPLICITLY_CONVERTIBLE(int64_t,T); From 7be279d2de6cc73c66f21d33f58fb6f41790f05c Mon Sep 17 00:00:00 2001 From: silviana amethyst Date: Wed, 8 Jul 2026 18:42:28 +0000 Subject: [PATCH 5/8] feat(bindings): np.rint / np.round on complex_mp (component-wise, half-to-even) numpy defines rint for complex dtypes (rounds real and imaginary parts independently) and np.round dispatches through it; complex_mp had no loop, so np.round(complex_mp array) raised. Register the component-wise rint for complex_mp, matching complex128 semantics exactly (half-to-even via mpfr_rint in MPFR_RNDN). Rounds out the np.abs/np.round ask of bertiniteam/b2#301. Co-Authored-By: Claude Fable 5 --- python/test/classes/numpy_ufuncs_test.py | 9 ++++++++ .../include/eigenpy_interaction.hpp | 22 +++++++++++++++---- 2 files changed, 27 insertions(+), 4 deletions(-) diff --git a/python/test/classes/numpy_ufuncs_test.py b/python/test/classes/numpy_ufuncs_test.py index 96edaf675..0406333e6 100644 --- a/python/test/classes/numpy_ufuncs_test.py +++ b/python/test/classes/numpy_ufuncs_test.py @@ -175,6 +175,15 @@ def test_rint_rounds_half_to_even(self): real_mp('-0.5'), real_mp('-2.5')]) assert [str(x) for x in np.rint(v)] == ['0', '2', '2', '-0', '-2'] + def test_rint_and_round_on_complex(self): + # numpy defines rint (and hence np.round) for complex, component-wise + w = np.array([complex_mp('1.5', '2.5'), complex_mp('-0.5', '3.4')]) + want = np.rint(np.array([1.5 + 2.5j, -0.5 + 3.4j])) + for got, ref in zip(np.rint(w), want): + assert complex(got) == ref + for got, ref in zip(np.round(w), want): + assert complex(got) == ref + def test_floor_ceil_trunc(self): v = np.array([real_mp('1.7'), real_mp('-1.7')]) assert [str(x) for x in np.floor(v)] == ['1', '-2'] diff --git a/python_bindings/include/eigenpy_interaction.hpp b/python_bindings/include/eigenpy_interaction.hpp index b8ff5c022..a83e017b8 100644 --- a/python_bindings/include/eigenpy_interaction.hpp +++ b/python_bindings/include/eigenpy_interaction.hpp @@ -285,15 +285,25 @@ namespace eigenpy struct op_trunc { template static T apply(T const& x) { return T(trunc(x)); } }; // numpy's rint is round-half-to-EVEN; boost's rint rounds half away from // zero, so call mpfr directly in MPFR_RNDN (nearest, ties to even). + // Unlike the rest of the rounding family, numpy defines rint for complex + // (component-wise) — np.round on a complex array goes through it. struct op_rint { - template static T apply(T const& x) + static bertini::real_mp rint_one(bertini::real_mp const& x) { - T out(0); + bertini::real_mp out(0); out.precision(x.precision()); mpfr_rint(out.backend().data(), x.backend().data(), MPFR_RNDN); return out; } + + template static T apply(T const& x) + { + if constexpr (is_complex_mp::value) + return at_precision_of(T(rint_one(x.real()), rint_one(x.imag())), x); + else + return rint_one(x); + } }; // ----- unary ops, cross-type output ----------------------------------- @@ -891,6 +901,10 @@ namespace eigenpy unary("isinf", &internal::guarded_unary_op_out, bool_code); unary("isfinite", &internal::guarded_unary_op_out, bool_code); + // rint is the one rounding ufunc numpy defines for complex too + // (component-wise) — np.round dispatches through it + unary("rint", &internal::guarded_unary_op, type_code); + if constexpr (WithOrderingComparitors) // the ordering-dependent set; NOT defined for complex types { binary("greater", &internal::guarded_compare_op, bool_code); @@ -922,11 +936,11 @@ namespace eigenpy &internal::guarded_mixed_compare_op, &internal::guarded_mixed_compare_op); - // real-only unary: rounding family, fabs, real-only transcendentals + // real-only unary: rounding family (sans rint, registered for both + // above), fabs, real-only transcendentals unary("floor", &internal::guarded_unary_op, type_code); unary("ceil", &internal::guarded_unary_op, type_code); unary("trunc", &internal::guarded_unary_op, type_code); - unary("rint", &internal::guarded_unary_op, type_code); unary("fabs", &internal::guarded_unary_op, type_code); unary("exp2", &internal::guarded_unary_op, type_code); unary("log2", &internal::guarded_unary_op, type_code); From 86c5688a6c08bfe2e301badecb053a6acf177afe Mon Sep 17 00:00:00 2001 From: silviana amethyst Date: Wed, 8 Jul 2026 19:23:55 +0000 Subject: [PATCH 6/8] feat(python): defend complex component access everywhere numpy lets us numpy's ndarray.real/.imag are C getsets gated on PyArray_ISCOMPLEX -- a hardwired builtin-type-number check (verified against numpy 2.4.x getset.c). No user dtype, legacy or NEP-42-style, can make the attributes correct: on complex_mp arrays .real returns the complex values and .imag returns zeros, silently. Decision: a crash is better than incorrect values. Defend every spelling within reach: - Solution overrides .real/.imag at the subclass level (a python property shadows the C getset), so solve results are simply correct -- for the mp dtype and for double solves alike. - bertini._numpy_guard, installed at import: np.real/np.imag raise a TypeError naming the right tool on plain mp-complex arrays (and on lists that would convert to them inside numpy); np.angle raises for EVERY mp-complex input, since it branches on dtype and dies in arctan2 where not even a subclass property can reach; all other inputs pass through to numpy untouched. - The raw .real/.imag attributes on a plain self-built ndarray remain the one spelling nothing can reach -- documented with a warning and a pinning test. Docs page and ADR-0051 updated accordingly. Co-Authored-By: Claude Fable 5 --- ...-numpy-ufunc-coverage-and-owned-getitem.md | 26 +++- python/bertini/__init__.py | 7 ++ python/bertini/_numpy_guard.py | 111 ++++++++++++++++++ python/bertini/records.py | 26 ++++ python/docs/source/numpy.rst | 26 ++-- python/test/classes/numpy_ufuncs_test.py | 68 +++++++++++ 6 files changed, 250 insertions(+), 14 deletions(-) create mode 100644 python/bertini/_numpy_guard.py diff --git a/docs/adr/0051-numpy-ufunc-coverage-and-owned-getitem.md b/docs/adr/0051-numpy-ufunc-coverage-and-owned-getitem.md index f86369029..1b22bfb01 100644 --- a/docs/adr/0051-numpy-ufunc-coverage-and-owned-getitem.md +++ b/docs/adr/0051-numpy-ufunc-coverage-and-owned-getitem.md @@ -96,12 +96,26 @@ that behavior faithfully. Consequences: 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, so the ndarray `.real`/`.imag` attributes - (and `np.real`/`np.imag`/`np.angle`) return silently wrong values — `.real` - gives the complex values, `.imag` gives zeros. Not hookable from the bindings; - documented, plus 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.) + 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.) ## Consequences diff --git a/python/bertini/__init__.py b/python/bertini/__init__.py index 509820f46..2e9760988 100644 --- a/python/bertini/__init__.py +++ b/python/bertini/__init__.py @@ -153,6 +153,13 @@ def __getattr__(name): # 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 diff --git a/python/bertini/_numpy_guard.py b/python/bertini/_numpy_guard.py new file mode 100644 index 000000000..e5b47b24e --- /dev/null +++ b/python/bertini/_numpy_guard.py @@ -0,0 +1,111 @@ +# This file is part of Bertini 2. +# +# python/bertini/_numpy_guard.py is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# python/bertini/_numpy_guard.py is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this file. If not, see . +# +# Copyright(C) Bertini2 Development Team + +"""Make numpy's component accessors fail LOUDLY on multiprecision complex arrays. + +numpy's ``ndarray.real`` / ``.imag`` are hardwired to its three built-in complex types +(``PyArray_ISCOMPLEX`` in ``getset.c``): on any other dtype, ``.real`` returns the array +itself and ``.imag`` returns zeros -- **silently wrong** for ``complex_mp``, and there is +no user-dtype hook to fix or even detect it at the numpy level. ``np.real`` / ``np.imag`` +/ ``np.angle`` are thin wrappers over those attributes and inherit the lie. + +The attributes themselves are C-level and untouchable, but the module *functions* are +plain Python -- so importing bertini wraps them: called on a plain ndarray of the +multiprecision complex dtype they raise ``TypeError`` (wrong-by-construction is converted +to loud), and every other input passes straight through to the original numpy functions. +A crash is better than incorrect values. + +Correct spellings, always: + +* ``bertini.real(x)`` / ``bertini.imag(x)`` / ``bertini.multiprec.arg(x)`` +* solution points (:class:`bertini.records.Solution`) override ``.real`` / ``.imag`` + at the subclass level and are simply correct -- the guard lets them through. + +The remaining untouchable spelling is the raw ``.real`` / ``.imag`` attribute on a plain +ndarray you built yourself; see the "Multiprecision numbers and NumPy" docs page. +""" + +import functools as _functools + +import numpy as _np + +from bertini.multiprec import complex_mp as _complex_mp + +_CPLX_MP = _np.dtype(_complex_mp) + +_MSG = ( + "numpy.{name}() cannot work on an array of the multiprecision complex dtype: " + "numpy's component access is hardwired to its built-in complex types and would " + "silently return WRONG values for user dtypes (there is no hook for bertini to fix " + "it). Use bertini.real(x) / bertini.imag(x) / bertini.multiprec.arg(x) instead -- " + "or hold a bertini Solution, whose .real/.imag are correct." +) + + +def _is_plain_complex_mp_array(x): + # exact-type check: subclasses (e.g. bertini's Solution) override .real/.imag + # correctly and must pass through + return type(x) is _np.ndarray and x.dtype == _CPLX_MP + + +def _would_lie(val): + if _is_plain_complex_mp_array(val): + return True + # a list/tuple of complex_mp converts to a plain mp-dtype array INSIDE numpy's + # real()/imag(), landing on the same wrong attribute path + if isinstance(val, (list, tuple)): + try: + return _is_plain_complex_mp_array(_np.asanyarray(val)) + except Exception: + return False + return False + + +def _is_complex_mp_anything(val): + # np.angle never consults .real/.imag attributes: it branches on the DTYPE and + # would die in arctan2 with a cryptic error for every mp-complex input -- + # Solutions included, since the subclass property cannot help it. Catch them + # all and say what to use instead. + if isinstance(val, _complex_mp): + return True + try: + return _np.asanyarray(val).dtype == _CPLX_MP + except Exception: + return False + + +def _guarded(orig, name, applies): + @_functools.wraps(orig) + def wrapper(val, *args, **kwargs): + if applies(val): + raise TypeError(_MSG.format(name=name)) + return orig(val, *args, **kwargs) + + wrapper._bertini_guarded = True + wrapper._bertini_original = orig + return wrapper + + +def install(): + """Wrap ``np.real`` / ``np.imag`` / ``np.angle``. Idempotent.""" + for name, applies in (("real", _would_lie), + ("imag", _would_lie), + ("angle", _is_complex_mp_anything)): + orig = getattr(_np, name) + if getattr(orig, "_bertini_guarded", False): + continue + setattr(_np, name, _guarded(orig, name, applies)) diff --git a/python/bertini/records.py b/python/bertini/records.py index 492c7340e..bd24e4037 100644 --- a/python/bertini/records.py +++ b/python/bertini/records.py @@ -144,6 +144,32 @@ def __array_wrap__(self, out_arr, context=None, return_scalar=False): result.annotations = {} return result + # numpy's ndarray .real/.imag are hardwired to its built-in complex types: on a + # multiprecision-complex array the base attributes return silently WRONG values + # (.real gives the complex values, .imag gives zeros), with no hook for user + # dtypes. A python subclass property CAN shadow the C-level attribute, so + # solution points -- the mp-complex arrays users actually hold -- are correct. + # Plain ndarrays are covered by the guarded np.real/np.imag/np.angle (see + # bertini._numpy_guard) and by bertini.real/imag. + + @property + def real(self): + """The real parts -- correct also for the multiprecision complex dtype.""" + from bertini.multiprec import complex_mp + if self.dtype == _np.dtype(complex_mp): + from bertini import _numpy_helpers as _nh + return _nh.real(_np.asarray(self)) + return _np.ndarray.real.__get__(self) + + @property + def imag(self): + """The imaginary parts -- correct also for the multiprecision complex dtype.""" + from bertini.multiprec import complex_mp + if self.dtype == _np.dtype(complex_mp): + from bertini import _numpy_helpers as _nh + return _nh.imag(_np.asarray(self)) + return _np.ndarray.imag.__get__(self) + class SolveResult: """What ``solve`` returns: the solutions plus a claim ticket on the recorded run. diff --git a/python/docs/source/numpy.rst b/python/docs/source/numpy.rst index 53f9ace82..49b2bdb9c 100644 --- a/python/docs/source/numpy.rst +++ b/python/docs/source/numpy.rst @@ -74,18 +74,28 @@ Component access on complex arrays =================================== The one place numpy itself cannot be taught about a user-defined dtype: the ndarray -attributes ``.real`` and ``.imag`` (and :func:`numpy.real` / :func:`numpy.imag` / -:func:`numpy.angle`, which route through them). +attributes ``.real`` and ``.imag``. numpy hardwires them to its own built-in complex +types -- on any other dtype ``.real`` returns the array itself and ``.imag`` returns +zeros, with no hook for the bindings to fix or even detect it. + +Bertini defends every spelling it can reach: + +* **Solution points are simply correct.** The arrays returned by ``bertini.solve`` + override ``.real`` / ``.imag`` at the subclass level, so ``sol.real`` / ``sol.imag`` + give the true components at full precision. +* **The numpy functions fail loudly.** Importing bertini wraps :func:`numpy.real` / + :func:`numpy.imag` / :func:`numpy.angle`: on a plain mp-complex array they raise a + ``TypeError`` naming the right tool, instead of silently returning wrong values. + All other inputs pass straight through to numpy. .. warning:: - On ``complex_mp`` **arrays**, ``.real`` returns the complex values themselves and - ``.imag`` returns zeros -- silently. numpy only knows its own built-in complex - types are complex-like; this is not hookable from the bindings. (Complex - *scalars* are fine: ``w[0].imag`` is exact.) + The raw ``.real`` / ``.imag`` **attributes on a plain ndarray** are the one + spelling nothing can reach: on a ``complex_mp`` array you built yourself (not a + Solution), they return wrong values silently. Complex *scalars* are fine -- + ``w[0].imag`` is exact. -Bertini provides the component accessors instead -- same results, full precision, -array-capable: +The component accessors -- same results, full precision, array-capable: .. doctest:: diff --git a/python/test/classes/numpy_ufuncs_test.py b/python/test/classes/numpy_ufuncs_test.py index 2f6998003..a01164b28 100644 --- a/python/test/classes/numpy_ufuncs_test.py +++ b/python/test/classes/numpy_ufuncs_test.py @@ -489,3 +489,71 @@ def test_ndarray_real_imag_attributes_are_untrustworthy(self): w = np.array([complex_mp(1, 2)]) assert w.real.dtype == np.dtype(complex_mp) # not real_mp! assert w.imag[0] == complex_mp(0) # wrong value, by numpy + + +class TestGuardedNumpyComponentFunctions: + """np.real/np.imag/np.angle raise on plain mp-complex arrays instead of + silently returning wrong values (bertini._numpy_guard) -- a crash is better + than incorrect values. Solution points override .real/.imag at the subclass + level and pass through correct.""" + + def test_np_real_imag_raise_on_plain_complex_mp_array(self): + w = np.array([complex_mp(1, 2), complex_mp(3, 4)]) + with pytest.raises(TypeError, match="bertini.real"): + np.real(w) + with pytest.raises(TypeError, match="bertini.real"): + np.imag(w) + with pytest.raises(TypeError, match="bertini.real"): + np.angle(w) + + def test_np_real_imag_raise_on_lists_of_complex_mp(self): + # a list converts to a plain mp array inside numpy, same wrong path + with pytest.raises(TypeError): + np.imag([complex_mp(1, 2)]) + + def test_guard_passes_everything_else_through(self): + # ordinary numpy is untouched + z = np.array([1 + 2j, 3 + 4j]) + assert list(np.real(z)) == [1.0, 3.0] + assert list(np.imag(z)) == [2.0, 4.0] + assert np.angle(np.array([1j]))[0] == pytest.approx(np.pi / 2) + # real_mp arrays are not complex: base semantics are already correct + v = np.array([real_mp(1), real_mp(2)]) + assert list(np.real(v)) == [real_mp(1), real_mp(2)] + assert list(np.imag(v)) == [real_mp(0), real_mp(0)] + # mp-complex SCALARS go through the (correct) scalar properties + assert np.real(complex_mp(1, 2)) == real_mp(1) + assert np.imag(complex_mp(1, 2)) == real_mp(2) + + def test_guard_is_idempotent(self): + import bertini._numpy_guard as guard + before = np.real + guard.install() + assert np.real is before + + def test_solution_real_imag_are_correct(self): + from bertini.records import Solution + s = Solution(np.array([complex_mp(1, 2), complex_mp(3, 4)])) + assert [str(x) for x in s.real] == ['1', '3'] + assert [str(x) for x in s.imag] == ['2', '4'] + assert s.real.dtype == np.dtype(real_mp) + # np.real/np.imag on a Solution route through the subclass property + assert [str(x) for x in np.real(s)] == ['1', '3'] + assert [str(x) for x in np.imag(s)] == ['2', '4'] + + def test_solution_real_imag_correct_for_double_solves_too(self): + from bertini.records import Solution + s = Solution(np.array([1 + 2j, 3 + 4j])) + assert list(s.real) == [1.0, 3.0] + assert list(s.imag) == [2.0, 4.0] + + def test_np_angle_raises_helpfully_for_all_mp_complex(self): + # np.angle branches on the DTYPE (never the .real/.imag attributes), so + # not even the Solution subclass can make it work -- the guard turns the + # cryptic arctan2 failure into a pointer at mp.arg, for every spelling + from bertini.records import Solution + for val in (np.array([complex_mp(1, 1)]), + Solution(np.array([complex_mp(1, 1)])), + complex_mp(1, 1)): + with pytest.raises(TypeError, match="arg"): + np.angle(val) From fc765e314e3687cb07341a0b343b12d2525c7033 Mon Sep 17 00:00:00 2001 From: silviana amethyst Date: Wed, 8 Jul 2026 19:36:40 +0000 Subject: [PATCH 7/8] feat(python): bertini.operators is the one-stop polymorphic math vocabulary `from bertini.operators import *` now gives functions that work on EVERYTHING, dispatched per argument: a symbolic Variable/expression builds a function-tree node; multiprecision scalars, mp-dtype numpy arrays, lists, and plain python numbers take the numeric path through the native precision-preserving loops. No more remembering that arg lives in bertini.multiprec while imag is at the top level. - sin/cos/tan/asin/acos/atan/exp/log/sqrt: fully polymorphic (symbolic twin exists). - sinh/cosh/tanh/asinh/acosh/atanh and abs/arg/real/imag/conj/round/ sum/norm/is_real: numeric, with a clear TypeError on symbolic input. - E/Pi/I ride along. abs/round/sum shadow the builtins only inside this opt-in star-import (they fall back to builtin behavior on plain python input); `from bertini import *` still never shadows. - the top-level elementary functions are rebound to the polymorphic versions (a strict superset: bertini.sin(x) now also accepts numbers and arrays), and bertini.arg + the hyperbolics join the top level. - new _numpy_helpers.arg rides mp.arg on complex arrays (real arrays go through the exact real->complex cast). Also fixes a latent #305 bug the new tests caught: _numpy_helpers' scalar fallbacks called bare abs()/round(), which resolve to the module's own shadowing functions at module scope -> infinite recursion on plain python input. Now explicitly builtins. Co-Authored-By: Claude Fable 5 --- python/bertini/__init__.py | 16 ++- python/bertini/_numpy_helpers.py | 27 ++++- python/bertini/operators.py | 113 +++++++++++++++++-- python/docs/source/numpy.rst | 15 +++ python/test/classes/operators_test.py | 153 ++++++++++++++++++++++++++ 5 files changed, 307 insertions(+), 17 deletions(-) create mode 100644 python/test/classes/operators_test.py diff --git a/python/bertini/__init__.py b/python/bertini/__init__.py index 2e9760988..840e73809 100644 --- a/python/bertini/__init__.py +++ b/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 @@ -163,12 +160,22 @@ def __getattr__(name): 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 @@ -178,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 @@ -189,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'] diff --git a/python/bertini/_numpy_helpers.py b/python/bertini/_numpy_helpers.py index 18f1be049..70d8aad81 100644 --- a/python/bertini/_numpy_helpers.py +++ b/python/bertini/_numpy_helpers.py @@ -34,6 +34,7 @@ these helpers or ``bertini.multiprec.real/imag/arg``. """ +import builtins as _builtins import numpy as _np from decimal import Decimal as _Decimal, ROUND_HALF_EVEN as _ROUND_HALF_EVEN @@ -88,7 +89,9 @@ def _imag_scalar(v): def _abs_scalar(v): if isinstance(v, (_complex_mp, _real_mp)): return _mp_abs(v) - return abs(v) + # explicitly the BUILTIN: this module's own `abs` shadows it at module scope, + # and the bare name recursed infinitely for plain python input + return _builtins.abs(v) def _conj_scalar(v): @@ -110,7 +113,7 @@ def _round_scalar(v, decimals): return _complex_mp(_round_real_mp(v.real, decimals), _round_real_mp(v.imag, decimals)) if isinstance(v, _real_mp): return _round_real_mp(v, decimals) - return round(v, decimals) + return _builtins.round(v, decimals) # explicitly the builtin (module `round` shadows it) # --- the public helpers ----------------------------------------------------------------------- @@ -153,6 +156,26 @@ def conj(x): return _elementwise(_conj_scalar, x) +def _arg_scalar(v): + if isinstance(v, _complex_mp): + return _mp.arg(v) + if isinstance(v, _real_mp): + return _mp.arg(_complex_mp(v)) + import cmath + return cmath.phase(complex(v)) + + +def arg(x): + """Argument(s) -- the angle from 0 -- as ``real_mp``, over a scalar / list / array + (the ``np.angle`` replacement; numpy's own cannot work on mp dtypes). Beware the + branch cut.""" + if _is_mp_array(x) and x.ndim == 1: + if x.dtype == _np.dtype(_real_mp): + x = x.astype(_np.dtype(_complex_mp)) # registered safe cast, exact + return _mp.arg(x) + return _elementwise(_arg_scalar, x) + + def round(x, decimals=0): """Round to ``decimals`` places, staying mp-native -- over a scalar / list / array (``np.round``).""" return _elementwise(lambda v: _round_scalar(v, decimals), x) diff --git a/python/bertini/operators.py b/python/bertini/operators.py index af44d7853..bbdb5e646 100644 --- a/python/bertini/operators.py +++ b/python/bertini/operators.py @@ -19,21 +19,112 @@ # as well as COPYING. Bertini2 is provided with permitted # additional terms in the b2/licenses/ directory. -"""The symbolic math vocabulary, gathered for ``from bertini.operators import *``. +"""The whole math vocabulary in one namespace -- symbols and numbers alike. -These are the elementary functions (``sin``, ``cos``, ...) and constants (``E``, ``Pi``, ``I``) used -to *build* symbolic systems on :class:`~bertini.Variable`\\ s. They also live at the top level -(``bertini.sin``, ``bertini.Pi``, ...); this module exists only so you can pull the math vocabulary -into your namespace without importing the rest of ``bertini``:: +One import gives you functions that work on *everything*: a symbolic +:class:`~bertini.Variable`/expression, a multiprecision number, a numpy array of them, +or a plain python number:: from bertini.operators import * - f = sin(x) + Pi*y - E -These are the *symbolic* operators; the numeric elementary functions (acting on multiprecision -numbers rather than expression nodes) live in :mod:`bertini.multiprec`. + f = sin(x) + Pi*y - E # symbolic (x, y Variables -> an expression) + v = sin(real_mp('0.5')) # numeric (full precision) + m = abs(solutions[0]) # numpy arrays (multiprecision dtypes included) + t = arg(complex_mp(1, 1)) # components (arg/real/imag/conj) + +Dispatch is by argument: a function-tree node builds a symbolic node; everything else +takes the numeric path (multiprecision scalars and mp-dtype numpy arrays go through the +native precision-preserving loops; python numbers and float arrays are plain numpy). +You never have to remember whether a name lives in ``bertini.multiprec`` or at the top +level -- it is here. + +The functions with no symbolic counterpart (``abs``, ``arg``, ``real``, ``imag``, +``conj``, ``round``, ``sum``, ``norm``, ``is_real``, and the hyperbolics) raise a clear +``TypeError`` when handed a symbolic expression. + +``abs``, ``round``, and ``sum`` shadow the python builtins **within your namespace** +when you star-import this module -- that is the point (they fall back to builtin +behavior on plain python input), but it is opt-in: ``from bertini import *`` never +shadows builtins. """ -# Re-exported from the top-level bertini package (defined there before this module is imported). -from . import sin, cos, tan, asin, acos, atan, exp, log, sqrt, E, Pi, I # noqa: F401 +import numpy as _np + +from bertini._pybertini.function_tree import AbstractNode as _AbstractNode + +from . import symbolics as _sym +from . import _numpy_helpers as _nh + +# the symbolic constants, ready to drop into expressions +from . import E, Pi, I # noqa: F401 + + +def _polymorphic(name, sym_fn, num_fn, doc): + def f(x): + if isinstance(x, _AbstractNode): + return sym_fn(x) + return num_fn(x) + f.__name__ = name + f.__qualname__ = name + f.__doc__ = doc + ("\n\nPolymorphic: builds a symbolic node for a function-tree " + "argument, computes numerically (precision-preserving, numpy " + "containers included) for everything else.") + return f + + +def _numeric_only(name, num_fn, doc): + def f(x, *args, **kwargs): + if isinstance(x, _AbstractNode): + raise TypeError( + f"{name}() is not defined for symbolic expressions -- it is a numeric " + "operation. Evaluate the expression first, or use the symbolic " + "functions (sin, cos, exp, ...) to build systems.") + return num_fn(x, *args, **kwargs) + f.__name__ = name + f.__qualname__ = name + f.__doc__ = doc + ("\n\nNumeric: multiprecision scalars, numpy arrays (mp dtypes " + "included), lists, and plain python numbers.") + return f + + +# --- the elementary functions with symbolic twins: full polymorphic dispatch ------------------- + +sin = _polymorphic('sin', _sym.sin, _np.sin, "Sine.") +cos = _polymorphic('cos', _sym.cos, _np.cos, "Cosine.") +tan = _polymorphic('tan', _sym.tan, _np.tan, "Tangent.") +asin = _polymorphic('asin', _sym.asin, _np.arcsin, "Arcsine.") +acos = _polymorphic('acos', _sym.acos, _np.arccos, "Arccosine.") +atan = _polymorphic('atan', _sym.atan, _np.arctan, "Arctangent.") +exp = _polymorphic('exp', _sym.exp, _np.exp, "Exponential, base e.") +log = _polymorphic('log', _sym.log, _np.log, "Natural logarithm.") +sqrt = _polymorphic('sqrt', _sym.sqrt, _np.sqrt, "Square root.") + +# --- numeric-only elementary functions (no symbolic node exists) ------------------------------- + +sinh = _numeric_only('sinh', _np.sinh, "Hyperbolic sine.") +cosh = _numeric_only('cosh', _np.cosh, "Hyperbolic cosine.") +tanh = _numeric_only('tanh', _np.tanh, "Hyperbolic tangent.") +asinh = _numeric_only('asinh', _np.arcsinh, "Hyperbolic arcsine.") +acosh = _numeric_only('acosh', _np.arccosh, "Hyperbolic arccosine.") +atanh = _numeric_only('atanh', _np.arctanh, "Hyperbolic arctangent.") + +# --- components, magnitudes, and friends (numeric-only) ---------------------------------------- + +abs = _numeric_only('abs', _nh.abs, "Magnitude(s), as real_mp for mp input.") # noqa: A001 +arg = _numeric_only('arg', _nh.arg, "Argument(s) (angle from 0), as real_mp. Beware the branch cut.") +real = _numeric_only('real', _nh.real, "Real part(s), as real_mp for mp input.") +imag = _numeric_only('imag', _nh.imag, "Imaginary part(s), as real_mp for mp input.") +conj = _numeric_only('conj', _nh.conj, "Complex conjugate(s).") +round = _numeric_only('round', _nh.round, "Round to N DECIMAL digits, staying mp-native.") # noqa: A001 +sum = _numeric_only('sum', _nh.sum, "Sum of a collection, staying mp-native.") # noqa: A001 +norm = _numeric_only('norm', _nh.norm, "Euclidean (2-)norm, as real_mp.") +is_real = _numeric_only('is_real', _nh.is_real, "Is every coordinate real (|imag| < tol)?") + -__all__ = ['sin', 'cos', 'tan', 'asin', 'acos', 'atan', 'exp', 'log', 'sqrt', 'E', 'Pi', 'I'] +__all__ = [ + 'sin', 'cos', 'tan', 'asin', 'acos', 'atan', + 'sinh', 'cosh', 'tanh', 'asinh', 'acosh', 'atanh', + 'exp', 'log', 'sqrt', + 'abs', 'arg', 'real', 'imag', 'conj', 'round', 'sum', 'norm', 'is_real', + 'E', 'Pi', 'I', +] diff --git a/python/docs/source/numpy.rst b/python/docs/source/numpy.rst index 49b2bdb9c..a587534b1 100644 --- a/python/docs/source/numpy.rst +++ b/python/docs/source/numpy.rst @@ -138,6 +138,21 @@ input they work element-wise. The builtin-shadowing names (``abs``, ``round``, never clobbers the Python builtins. For the tolerance point comparison behind de-duplication, see :func:`bertini.is_distinct_up_to`. +And for the whole vocabulary in one go -- the elementary functions *and* these +helpers, working on symbolic expressions, multiprecision numbers, and numpy +containers alike, dispatched per argument:: + + from bertini.operators import * + + f = sin(x) + Pi*y # symbolic (x, y Variables) + v = sin(real_mp('0.5')) # numeric, full precision + m = abs(solutions[0]) # numpy containers, mp dtypes included + t = arg(complex_mp(1, 1)) # components: arg/real/imag/conj + +This star-import *does* shadow ``abs``/``round``/``sum`` in your namespace -- that is +its point (they fall back to builtin behavior on plain python input), and it is +opt-in. + Boundaries by design ===================== diff --git a/python/test/classes/operators_test.py b/python/test/classes/operators_test.py new file mode 100644 index 000000000..5c54eff96 --- /dev/null +++ b/python/test/classes/operators_test.py @@ -0,0 +1,153 @@ +# This file is part of Bertini 2. +# +# python/test/classes/operators_test.py is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# python/test/classes/operators_test.py is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this file. If not, see . +# +# Copyright(C) Bertini2 Development Team + +"""bertini.operators: one star-import, functions that work on symbols AND numbers AND +numpy containers alike, dispatched per argument.""" + +import builtins + +import numpy as np +import pytest + +import bertini as pb +import bertini.multiprec as mp +import bertini.operators as ops +from bertini.multiprec import complex_mp, real_mp +from bertini._pybertini.function_tree import AbstractNode + + +class TestPolymorphicDispatch: + """sin & friends: symbolic on expressions, numeric on everything else.""" + + def test_symbolic_on_variables(self): + x = pb.Variable('x') + f = ops.sin(x) + ops.Pi * x - ops.E + assert isinstance(ops.sin(x), AbstractNode) + assert isinstance(f, AbstractNode) + + def test_numeric_on_mp_scalars(self): + v = real_mp('0.5') + assert ops.sin(v) == mp.sin(v) + assert ops.exp(v) == mp.exp(v) + z = complex_mp('0.5', '0.25') + assert ops.sqrt(z) == mp.sqrt(z) + + def test_numeric_on_mp_arrays(self): + v = np.array([real_mp('0.25'), real_mp('0.5')]) + out = ops.cos(v) + assert out.dtype == np.dtype(real_mp) + assert out[0] == mp.cos(v[0]) + + def test_numeric_on_lists_and_python_numbers(self): + assert ops.sin(0.0) == 0.0 + out = ops.tan([real_mp('0.25'), real_mp('0.5')]) + assert out[1] == mp.tan(real_mp('0.5')) + + def test_asin_maps_to_arcsin(self): + v = real_mp('0.5') + assert ops.asin(v) == mp.asin(v) + assert ops.acos(v) == mp.acos(v) + assert ops.atan(v) == mp.atan(v) + + def test_hyperbolics_numeric(self): + v = real_mp('0.5') + assert ops.sinh(v) == mp.sinh(v) + assert ops.atanh(v) == mp.atanh(v) + + def test_hyperbolics_reject_symbols(self): + x = pb.Variable('x') + with pytest.raises(TypeError, match="symbolic"): + ops.sinh(x) + + +class TestComponentsAndFriends: + """abs/arg/real/imag/conj/round/sum/norm/is_real, all in the same namespace.""" + + def test_components_on_arrays(self): + w = np.array([complex_mp(1, 2), complex_mp(3, 4)]) + assert [str(t) for t in ops.real(w)] == ['1', '3'] + assert [str(t) for t in ops.imag(w)] == ['2', '4'] + assert str(ops.abs(np.array([complex_mp(3, 4)]))[0]) == '5' + assert ops.arg(w)[0] == mp.arg(w[0]) + assert complex(ops.conj(w)[1]) == complex(3, -4) + + def test_arg_on_scalars_and_reals(self): + assert ops.arg(complex_mp(0, 1)) == mp.arg(complex_mp(0, 1)) + # arg of a negative real is pi + assert mp.abs(ops.arg(np.array([real_mp(-2)]))[0] - mp.arg(complex_mp(-2))) == 0 + # plain python numbers give floats + assert ops.arg(1j) == pytest.approx(np.pi / 2) + + def test_sum_norm_is_real(self): + v = np.array([real_mp(3), real_mp(4)]) + assert ops.sum(v) == real_mp(7) + assert ops.norm(v) == real_mp(5) + assert ops.is_real(np.array([complex_mp(1)])) is True + + def test_round_stays_decimal(self): + assert str(ops.round(real_mp('2.34567'), 2)) == '2.35' + + def test_numeric_only_reject_symbols(self): + x = pb.Variable('x') + for fn in (ops.abs, ops.arg, ops.real, ops.imag, ops.conj, + ops.round, ops.sum, ops.norm, ops.is_real): + with pytest.raises(TypeError, match="symbolic"): + fn(x) + + def test_builtin_fallback_on_plain_python(self): + # the shadowing names still behave sanely on plain python input + assert ops.abs(-3) == 3 + assert ops.sum([1, 2, 3]) == 6 + assert ops.round(2.345, 1) == builtins.round(2.345, 1) + + +class TestStarImportSurface: + def test_star_import_gives_the_whole_vocabulary(self): + ns = {} + exec("from bertini.operators import *", ns) + for name in ('sin', 'cos', 'tan', 'asin', 'acos', 'atan', + 'sinh', 'cosh', 'tanh', 'asinh', 'acosh', 'atanh', + 'exp', 'log', 'sqrt', + 'abs', 'arg', 'real', 'imag', 'conj', + 'round', 'sum', 'norm', 'is_real', + 'E', 'Pi', 'I'): + assert name in ns, name + + def test_one_import_covers_symbols_and_numbers(self): + # the point of the module, as a single flow + ns = {} + exec("from bertini.operators import *", ns) + x = pb.Variable('x') + f = ns['sin'](x) # symbolic + assert isinstance(f, AbstractNode) + val = ns['sin'](real_mp('0.5')) # mp scalar + assert val == mp.sin(real_mp('0.5')) + w = np.array([complex_mp(1, 2)]) # numpy container + assert ns['imag'](w)[0] == real_mp(2) + + def test_top_level_functions_are_polymorphic_too(self): + x = pb.Variable('x') + assert isinstance(pb.sin(x), AbstractNode) + assert pb.sin(real_mp('0.5')) == mp.sin(real_mp('0.5')) + assert pb.arg(complex_mp(1, 1)) == mp.arg(complex_mp(1, 1)) + + def test_from_bertini_star_still_never_shadows_builtins(self): + ns = {} + exec("from bertini import *", ns) + assert 'abs' not in ns + assert 'round' not in ns + assert 'sum' not in ns From eb238eef6a28a219e4595a597eee91799d6f7fb3 Mon Sep 17 00:00:00 2001 From: silviana amethyst Date: Wed, 8 Jul 2026 23:48:18 +0000 Subject: [PATCH 8/8] fix(bindings): never free or write through numpy-owned slots (numpy 2.5 reduce UAF) CI segfaulted on Linux py3.12-3.14 (numpy 2.5.1; py3.10/3.11 got numpy 2.2/2.4 and passed). Root cause, pinned locally with valgrind after reproducing in a numpy-2.5.1 venv: numpy 2.5 initializes the accumulator of an IDENTITYLESS reduce (np.min/np.max) by memcpy of element 0, so the accumulator slot and v[0] share one mpfr allocation. Our loops' `res = Op::apply(...)` move-assignment freed the slot's old limbs -- freeing v[0]'s storage out from under it (np.max(v) also silently rewrote v[0] through the shared limbs: v[0] went from 3 to 3.5). The freed block was re-read by the next reduce: use-after-free, double-free, corrupted mimalloc metadata, and a SIGSEGV three tests later inside np.median -- the classic action-at-a-distance crash. Fix: slot_write is now the only sanctioned store into an mp output slot. It computes the value first (the slot may also be an input), memsets the slot to BMP's uninitialized sentinel, and move-assigns the fresh value in -- no existing allocation is ever freed or written through. Applied to every loop that writes mp slots (binary, unary, cross-type unary, matmul, dotfunc); bool outputs stay plain stores. Same crash-into-bounded-leak doctrine as HardenSetitem (ADR-0006); recorded as ADR-0051 decision 7. Verified: numpy 2.5.1 -- min/max/median correct, input arrays intact, valgrind 0 errors (previously UAF pair between the maximum and minimum loops + fatal invalid read inside mi_malloc); numpy 2.4.6 -- full suite green. New regression test runs the reduces repeatedly and asserts the input survives. Also fixes the all-Windows wheel failure (separate, infrastructural): clang rejects a PCH whose input mtime changed even with identical content, which happens when the wheel build re-configures into a shared bld/. Add -Xclang -fno-pch-timestamp to the two PCH targets for clang builds. Co-Authored-By: Claude Fable 5 --- core/CMakeLists.txt | 9 +++++ ...-numpy-ufunc-coverage-and-owned-getitem.md | 18 +++++++++ python/test/classes/numpy_ufuncs_test.py | 15 +++++++ python_bindings/CMakeLists.txt | 7 ++++ .../include/eigenpy_interaction.hpp | 40 +++++++++++++++---- 5 files changed, 82 insertions(+), 7 deletions(-) diff --git a/core/CMakeLists.txt b/core/CMakeLists.txt index e5b171634..0f558f1d2 100644 --- a/core/CMakeLists.txt +++ b/core/CMakeLists.txt @@ -570,6 +570,15 @@ target_precompile_headers(bertini2_exe PRIVATE ) +# 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 "$<$:SHELL:-Xclang -fno-pch-timestamp>") +endif() + # todo: this should be made a devmode thing #target_compile_options(bertini2 PRIVATE -Wall -Wextra) diff --git a/docs/adr/0051-numpy-ufunc-coverage-and-owned-getitem.md b/docs/adr/0051-numpy-ufunc-coverage-and-owned-getitem.md index 1b22bfb01..c1dccf848 100644 --- a/docs/adr/0051-numpy-ufunc-coverage-and-owned-getitem.md +++ b/docs/adr/0051-numpy-ufunc-coverage-and-owned-getitem.md @@ -117,6 +117,24 @@ that behavior faithfully. Consequences: (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 diff --git a/python/test/classes/numpy_ufuncs_test.py b/python/test/classes/numpy_ufuncs_test.py index a01164b28..75f72bbe4 100644 --- a/python/test/classes/numpy_ufuncs_test.py +++ b/python/test/classes/numpy_ufuncs_test.py @@ -248,6 +248,21 @@ def test_argmax_argmin_max_min(self): assert str(np.max(v)) == '7' assert str(np.min(v)) == '1' + def test_identityless_reduce_does_not_corrupt_input(self): + # numpy >= 2.5 initializes the accumulator of an identityless reduce + # (np.min/np.max) by BITWISE copy of element 0, so the accumulator and + # v[0] share one mpfr allocation. The loops must never free or write + # through an output slot's existing allocation (slot_write): before + # that rule, np.max silently rewrote v[0] and freed its storage, and + # the next reduce crashed the interpreter (use-after-free -> corrupted + # allocator). Values AND the input array must survive, repeatedly. + v = np.array([real_mp(3), real_mp(1), real_mp(7), real_mp(2)]) + for _ in range(3): + assert str(np.max(v)) == '7' + assert [str(x) for x in v] == ['3', '1', '7', '2'] + assert str(np.min(v)) == '1' + assert [str(x) for x in v] == ['3', '1', '7', '2'] + def test_argmax_nan_wins(self): # numpy float semantics: the first nan is the arg-extremum v = np.array([real_mp(1), real_mp('nan'), real_mp(3)]) diff --git a/python_bindings/CMakeLists.txt b/python_bindings/CMakeLists.txt index 8252361e5..a20f30372 100644 --- a/python_bindings/CMakeLists.txt +++ b/python_bindings/CMakeLists.txt @@ -293,6 +293,13 @@ endif() # from scratch, causing peak-RAM spikes during parallel builds and inflating log output. target_precompile_headers(_pybertini PRIVATE "include/python_common.hpp") +# clang rejects a PCH whose input file's mtime changed, even with identical content -- +# which happens when the wheel build re-configures into a shared bld/. See the note +# beside bertini2_exe's PCH in core/CMakeLists.txt. +if(CMAKE_CXX_COMPILER_ID MATCHES "Clang") + target_compile_options(_pybertini PRIVATE "$<$:SHELL:-Xclang -fno-pch-timestamp>") +endif() + cmake_print_variables(SKBUILD) if(${SKBUILD}) # see https://stackoverflow.com/questions/1242904/finding-python-site-packages-directory-with-cmake diff --git a/python_bindings/include/eigenpy_interaction.hpp b/python_bindings/include/eigenpy_interaction.hpp index a83e017b8..00d8219fc 100644 --- a/python_bindings/include/eigenpy_interaction.hpp +++ b/python_bindings/include/eigenpy_interaction.hpp @@ -167,8 +167,31 @@ namespace eigenpy // These replace eigenpy's EIGENPY_REGISTER_{BINARY,UNARY}_UFUNC loop bodies // (and its gufunc_matrix_multiply), which read input slots unguarded and // segfault inside libmpfr/libmpc on never-written np.zeros/np.empty slots. - // Writes into the output slot go through BMP operator=, which initializes - // a zeroed destination itself. + // + // Writing an output slot: NEVER through BMP operator= on the slot's + // existing value. 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 assignment move-frees the slot's old limbs (freeing v[0]'s + // storage out from under it: use-after-free, double-free, corrupted + // allocator, SIGSEGV a few calls later) and writes through shared + // storage (silently mutating v[0]'s value). slot_write below is the + // only sanctioned store: compute the value FIRST (the slot may also be + // an input), memset the slot to BMP's uninitialized sentinel, then + // move the fresh value in -- no existing allocation is ever freed or + // written through. If the slot held a uniquely-owned value, its + // allocation leaks (numpy never destructs user-dtype elements anyway); + // same crash-into-bounded-leak trade as HardenSetitem / ADR-0006. + + // store `val` into a numpy-managed mp slot without freeing or writing + // through the slot's existing (possibly aliased) allocation. + template + inline void slot_write(T& slot, T val) + { + std::memset(static_cast(&slot), 0, sizeof(T)); + slot = std::move(val); // move into the sentinel: steals val's limbs, frees nothing + } struct op_add { template static T apply(T const& x, T const& y) { return T(x + y); } }; struct op_subtract { template static T apply(T const& x, T const& y) { return T(x - y); } }; @@ -435,7 +458,7 @@ namespace eigenpy T const& x = value_or_zero(*reinterpret_cast(i0), zero); T const& y = value_or_zero(*reinterpret_cast(i1), zero); T& res = *reinterpret_cast(o); - res = Op::apply(x, y); + slot_write(res, Op::apply(x, y)); // never plain operator= -- see slot_write i0 += is0; i1 += is1; o += os; @@ -512,7 +535,7 @@ namespace eigenpy { T const& x = value_or_zero(*reinterpret_cast(i), zero); T& res = *reinterpret_cast(o); - res = Op::apply(x); + slot_write(res, Op::apply(x)); // never plain operator= -- see slot_write i += is; o += os; } @@ -534,7 +557,10 @@ namespace eigenpy { T const& x = value_or_zero(*reinterpret_cast(i), zero); OutT& res = *reinterpret_cast(o); - res = Op::apply(x); + if constexpr (std::is_trivially_copyable_v) + res = Op::apply(x); + else + slot_write(res, Op::apply(x)); // never plain operator= -- see slot_write i += is; o += os; } @@ -567,7 +593,7 @@ namespace eigenpy b += is2_n; } T& res = *reinterpret_cast(op); - res = sum; + slot_write(res, std::move(sum)); // never plain operator= -- see slot_write ip2 += is2_p; op += os_p; } @@ -623,7 +649,7 @@ namespace eigenpy p0 += is0; p1 += is1; } - *reinterpret_cast(op) = acc; + slot_write(*reinterpret_cast(op), std::move(acc)); // never plain operator= -- see slot_write } // guarded element comparison for the PyArray_ArrFuncs `compare` slot