From 0f0832ecda8ff546b69892722c2322092a00d8a7 Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Wed, 14 Oct 2026 16:06:40 +0200 Subject: [PATCH 1/9] Add oneMKL power loops for float32 and float64 Arrays use v?Pow and a scalar exponent uses v?Powx, above the usual 8192-element threshold. The loop keeps NumPy's semantics: - the scalar exponent fast path of NumPy >= 2.3 (-1, 0, 0.5, 1, 2) is mirrored when built against it - v?Powx gets -0.0 and -inf wrong for exponents -0.5, 0.5 and 1.5, so those elements are recomputed with pow - MKL raises spurious invalid/overflow flags for some special inputs and drops flags raised in its worker threads, so when the VM error status or the flags report a problem they are recomputed with pow - in-place calls go through a buffer so the inputs survive for that Complex loops are left to NumPy: v?Pow is not exact for integer exponents (1j**2). --- CHANGELOG.md | 1 + mkl_umath/generate_umath.py | 11 +- mkl_umath/src/mkl_umath_loops.c.src | 158 ++++++++++++++++++++++++++++ mkl_umath/src/mkl_umath_loops.h.src | 2 +- mkl_umath/tests/test_basic.py | 72 ++++++++++++- 5 files changed, 241 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1e2fdd3e..088e4e00 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [0.6.0] - 2026-10-15 ### Added +* Added oneMKL-backed `power` loops for float32 and float64 arrays, using `v?Pow` and `v?Powx` for a scalar exponent * Enabled support of Python 3.15 [gh-279](https://github.com/IntelPython/mkl_umath/pull/279) * Added support for free-threaded (GIL-disabled) CPython builds: the Cython extension is compiled with `freethreading_compatible=True` and `_ufuncs` declares `Py_MOD_GIL_NOT_USED`, so importing `mkl_umath` no longer re-enables the GIL [gh-255](https://github.com/IntelPython/mkl_umath/pull/255) diff --git a/mkl_umath/generate_umath.py b/mkl_umath/generate_umath.py index ca6a578b..15a80886 100644 --- a/mkl_umath/generate_umath.py +++ b/mkl_umath/generate_umath.py @@ -501,7 +501,16 @@ def english_upper(s): TD(inexactvec + cmplxvec), ), # '_ones_like': - # 'power': + # 'power': complex loops are not added since v?Pow computes exp(y*log(x)), + # which is not exact for integer exponents (e.g. 1j**2) + "power": Ufunc( + 2, + 1, + None, + docstrings.get("numpy._core.umath.power"), + None, + TD(inexactvec), + ), # 'float_power': "absolute": Ufunc( 1, diff --git a/mkl_umath/src/mkl_umath_loops.c.src b/mkl_umath/src/mkl_umath_loops.c.src index c1a41660..fbbe432e 100644 --- a/mkl_umath/src/mkl_umath_loops.c.src +++ b/mkl_umath/src/mkl_umath_loops.c.src @@ -71,6 +71,11 @@ #define VML_ASM_THRESHOLD 100000 #define VML_D_THRESHOLD 8000 +/* floating-point exceptions reported by NumPy that v?Pow does not raise reliably */ +#define POW_FPE_FLAGS (FE_INVALID | FE_DIVBYZERO | FE_OVERFLOW) +/* block size, in elements, for in-place power, which goes through a buffer */ +#define POW_INPLACE_BLOCK ((npy_intp) 1 << 16) + #define MKL_INT_MAX ((npy_intp) ((~((MKL_UINT) 0)) >> 1)) #define CHUNKED_VML_CALL2(vml_func, n, type, mkl_type, mkl_ftype, in1, op1) \ @@ -107,6 +112,22 @@ } \ } while(0) +#define CHUNKED_VML_POWX_CALL(vml_func, n, type, in1, in2, op1) \ + do { \ + npy_intp _n_ = (n); \ + const npy_intp _chunk_size = MKL_INT_MAX; \ + type *in1p = (type *) (in1); \ + const type _in2 = *((type *) (in2)); \ + type *op1p = (type *) (op1); \ + while (_n_ > 0) { \ + npy_intp _current_chunk = (_n_ > _chunk_size) ? _chunk_size : _n_; \ + vml_func((MKL_INT) _current_chunk, in1p, _in2, op1p); \ + _n_ -= _current_chunk; \ + in1p += _current_chunk; \ + op1p += _current_chunk; \ + } \ + } while(0) + #define CHUNKED_VML_LINEARFRAC_CALL(vml_func, n, type, in1, op1, scaleA, shiftA, scaleB, shiftB) \ do { \ @@ -837,6 +858,143 @@ mkl_umath_@TYPE@_@kind@(char **args, const npy_intp *dimensions, const npy_intp } /**end repeat1**/ +/* + * v@s@Pow and v@s@Powx return the same values as pow@c@, except v@s@Powx for + * -0.0 and -inf raised to -0.5, 0.5 and 1.5, which are recomputed. They do + * not raise the same floating-point exception flags though: they raise + * spurious invalid and overflow flags for some special inputs (e.g. 0 ** -1), + * and flags raised in MKL worker threads never reach the calling thread. The + * VM error status does reach it, so whenever MKL reports a problem the flags + * are recomputed with pow@c@. r must not alias x or y, which are read again. + */ +static void +@TYPE@_power_vml(const @type@ *x, const @type@ *y, const int y_is_scalar, @type@ *r, const npy_intp n) +{ + /* flags raised by earlier calls of the loop are NumPy's to report */ + const int prev_flags = fetestexcept(POW_FPE_FLAGS); + + feclearexcept(POW_FPE_FLAGS); + vmlClearErrStatus(); + if (y_is_scalar) { + CHUNKED_VML_POWX_CALL(v@s@Powx, n, @type@, x, y, r); + } + else { + CHUNKED_VML_CALL3(v@s@Pow, n, @type@, @type@, x, y, r); + } + if (vmlGetErrStatus() != VML_STATUS_OK || fetestexcept(POW_FPE_FLAGS)) { + npy_intp i; + + feclearexcept(POW_FPE_FLAGS); + for (i = 0; i < n; i++) { + /* volatile keeps the calls, which are only made for their flags */ + volatile @type@ res = pow@c@(x[i], y_is_scalar ? *y : y[i]); + (void) res; + } + } + if (y_is_scalar && (*y == -0.5 || *y == 0.5 || *y == 1.5)) { + const @type@ ninf = -(@type@) NPY_INFINITY; + int candidates = 0; + npy_intp i; + + /* branch-free so that it vectorizes, the fix-up below rarely runs */ + for (i = 0; i < n; i++) { + candidates |= (x[i] == 0) | (x[i] == ninf); + } + if (candidates) { + for (i = 0; i < n; i++) { + if ((x[i] == 0 || x[i] == ninf) && signbit(x[i])) { + r[i] = pow@c@(x[i], *y); + } + } + } + } + if (prev_flags) { + feraiseexcept(prev_flags); + } +} + +void +mkl_umath_@TYPE@_power(char **args, const npy_intp *dimensions, const npy_intp *steps, void *NPY_UNUSED(func)) +{ + const npy_intp n = dimensions[0]; + +#if defined(NPY_2_3_API_VERSION) && NPY_API_VERSION >= NPY_2_3_API_VERSION + /* same fast path for a scalar exponent as in NumPy >= 2.3 */ + if (steps[1] == 0) { + const @type@ in2 = *(@type@ *)args[1]; + + if (in2 == -1.0 || in2 == 0.0 || in2 == 0.5 || in2 == 1.0 || in2 == 2.0) { + BINARY_LOOP { + const @type@ in1 = *(@type@ *)ip1; + + if (in2 == -1.0) { + *(@type@ *)op1 = 1.0 / in1; + } + else if (in2 == 0.0) { + *(@type@ *)op1 = 1.0; + } + else if (in2 == 0.5) { + *(@type@ *)op1 = sqrt@c@(in1); + } + else if (in2 == 1.0) { + *(@type@ *)op1 = in1; + } + else { + *(@type@ *)op1 = in1 * in1; + } + } + return; + } + } +#endif + + { + const int contig = IS_BINARY_CONT(@type@, @type@); + const int contig_scalar = IS_BINARY_CONT_S2(@type@, @type@); + const int disjoint_or_same1 = DISJOINT_OR_SAME(args[0], args[2], n, sizeof(@type@)); + const int disjoint_or_same2 = DISJOINT_OR_SAME(args[1], args[2], n, sizeof(@type@)); + + if (n > VML_TRANSCEDENTAL_THRESHOLD && disjoint_or_same1 && + ((contig && disjoint_or_same2) || contig_scalar)) { + const @type@ *x = (const @type@ *)args[0]; + const @type@ *y = (const @type@ *)args[1]; + @type@ *r = (@type@ *)args[2]; + const int y_is_scalar = !contig; + + if (r != x && (y_is_scalar || r != y)) { + @TYPE@_power_vml(x, y, y_is_scalar, r, n); + return; + } + else { + /* in-place: keep the inputs intact until their block is done */ + const npy_intp block = (n < POW_INPLACE_BLOCK) ? n : POW_INPLACE_BLOCK; + @type@ *buf = (@type@ *) malloc(block * sizeof(@type@)); + + if (buf != NULL) { + npy_intp off; + + for (off = 0; off < n; off += block) { + const npy_intp len = (n - off < block) ? n - off : block; + + @TYPE@_power_vml(x + off, y_is_scalar ? y : y + off, y_is_scalar, buf, len); + memcpy(r + off, buf, len * sizeof(@type@)); + } + free(buf); + return; + } + } + } + } + + { + BINARY_LOOP { + const @type@ in1 = *(@type@ *)ip1; + const @type@ in2 = *(@type@ *)ip2; + *((@type@ *)op1) = pow@c@(in1, in2); + } + } +} + /**begin repeat1 * # kind = cos, sin, tan, arccos, arcsin, arctan, cosh, sinh, tanh, arccosh, arcsinh, arctanh, fabs, floor, ceil, rint, trunc, cbrt, sqrt, expm1, log, log2, log1p, log10# * # func = cos, sin, tan, acos, asin, atan, cosh, sinh, tanh, acosh, asinh, atanh, fabs, floor, ceil, rint, trunc, cbrt, sqrt, expm1, log, log2, log1p, log10# diff --git a/mkl_umath/src/mkl_umath_loops.h.src b/mkl_umath/src/mkl_umath_loops.h.src index 087a97cb..4c97ea1e 100644 --- a/mkl_umath/src/mkl_umath_loops.h.src +++ b/mkl_umath/src/mkl_umath_loops.h.src @@ -58,7 +58,7 @@ copysign, cos, cosh, divide, divmod, isfinite, isinf, isnan, equal, exp, exp2, expm1, fabs, floor, fmax, fmin, frexp, greater, greater_equal, ldexp, less, less_equal, log, log2, log10, log1p, logical_and, logical_not, logical_or, logical_xor, multiply, modf, - negative, nextafter, not_equal, positive, reciprocal, rint, sign, signbit, sin, + negative, nextafter, not_equal, positive, power, reciprocal, rint, sign, signbit, sin, sinh, spacing, sqrt, square, subtract, tan, tanh, trunc# */ MKL_UMATH_API diff --git a/mkl_umath/tests/test_basic.py b/mkl_umath/tests/test_basic.py index 8433428c..477dce67 100644 --- a/mkl_umath/tests/test_basic.py +++ b/mkl_umath/tests/test_basic.py @@ -23,6 +23,8 @@ # OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE # OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +import warnings + import numpy as np import pytest @@ -70,7 +72,7 @@ def get_args(args_str, size, low, high): if umath in ["arccos", "arcsin", "arctanh"]: low = -1 high = 1 - elif umath in ["log", "log10", "log1p", "log2", "sqrt"]: + elif umath in ["log", "log10", "log1p", "log2", "power", "sqrt"]: low = 0 high = 10 elif umath == "arccosh": @@ -202,3 +204,71 @@ def test_absolute_nan_signbit(size, dtype): f"absolute({sign}nan) must be +nan for " f"{dtype.__name__} at size {size}" ) + + +_power_specials = [0.0, -0.0, 0.5, -0.5, 1.0, -1.0, 1.5, -1.5, 2.0, -2.0] +_power_specials += [3.0, -3.0, 1e-30, 1e30, np.inf, -np.inf, np.nan] + + +def _assert_power_equal(mkl_res, np_res): + # MKL may differ from libm by 1 ULP, special values must match exactly + eps = np.spacing(np_res.dtype.type(1)) + np.testing.assert_allclose(mkl_res, np_res, rtol=2 * eps, equal_nan=True) + assert (np.signbit(mkl_res) == np.signbit(np_res))[~np.isnan(np_res)].all() + + +def _power_warnings(func, *args): + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter("always") + func(*args) + return sorted(str(x.message) for x in w) + + +@pytest.mark.parametrize("size", [100, 8192 + 1]) +@pytest.mark.parametrize("dtype", [np.float32, np.float64]) +def test_power_special_values(size, dtype): + # size 100 takes the fallback loop, 8193 takes MKL's v?Pow and v?Powx + specials = np.array(_power_specials, dtype=dtype) + x, y = (a.ravel() for a in np.meshgrid(specials, specials)) + with np.errstate(all="ignore"): + x = np.resize(x, max(size, x.size)) + y = np.resize(y, max(size, y.size)) + _assert_power_equal(mu.power(x, y), np.power(x, y)) + for e in specials: + _assert_power_equal(mu.power(x, e), np.power(x, e)) + + +@pytest.mark.parametrize("dtype", [np.float32, np.float64]) +@pytest.mark.parametrize("exponent", [None, 2.5, 1.5, -0.5]) +def test_power_inplace(dtype, exponent): + x = np.random.uniform(0, 10, 10**5 + 1).astype(dtype) + y = np.random.uniform(-3, 3, 10**5 + 1).astype(dtype) + if exponent is not None: + y = dtype(exponent) + np_res = np.power(x, y) + + mkl_res = x.copy() + mu.power(mkl_res, y, out=mkl_res) + _assert_power_equal(mkl_res, np_res) + if exponent is None: + mkl_res = y.copy() + mu.power(x, mkl_res, out=mkl_res) + _assert_power_equal(mkl_res, np_res) + + +@pytest.mark.parametrize("dtype", [np.float32, np.float64]) +@pytest.mark.parametrize( + "base, exponent", + [(-1.0, 1.5), (0.0, -1.0), (-0.0, -0.5), (1e30, 3.0), (0.0, -np.inf)], +) +@pytest.mark.parametrize("pos", [0, -1]) +def test_power_warnings(dtype, base, exponent, pos): + # MKL raises spurious flags for some of these, and flags raised in its + # worker threads (for the element at the end) never reach the caller + x = np.full(10**6, 2.0, dtype=dtype) + y = np.full(10**6, 1.5, dtype=dtype) + x[pos], y[pos] = base, exponent + assert _power_warnings(mu.power, x, y) == _power_warnings(np.power, x, y) + assert _power_warnings(mu.power, x, y[pos]) == _power_warnings( + np.power, x, y[pos] + ) From c4373dfc74dedaed8757f27b8869c7a544eda4fb Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Wed, 14 Oct 2026 16:32:10 +0200 Subject: [PATCH 2/9] Add PR link to the changelog entry --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 088e4e00..fd10468b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,7 +7,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [0.6.0] - 2026-10-15 ### Added -* Added oneMKL-backed `power` loops for float32 and float64 arrays, using `v?Pow` and `v?Powx` for a scalar exponent +* Added oneMKL-backed `power` loops for float32 and float64 arrays, using `v?Pow` and `v?Powx` for a scalar exponent [gh-285](https://github.com/IntelPython/mkl_umath/pull/285) * Enabled support of Python 3.15 [gh-279](https://github.com/IntelPython/mkl_umath/pull/279) * Added support for free-threaded (GIL-disabled) CPython builds: the Cython extension is compiled with `freethreading_compatible=True` and `_ufuncs` declares `Py_MOD_GIL_NOT_USED`, so importing `mkl_umath` no longer re-enables the GIL [gh-255](https://github.com/IntelPython/mkl_umath/pull/255) From 5143651b6746d0da7bbae274811db76c6f422409 Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Wed, 14 Oct 2026 16:54:06 +0200 Subject: [PATCH 3/9] Make power's FP warnings independent of the CPU test_power_warnings failed on some CI runners: it compared against np.power, whose warnings for edge cases such as 0 ** -inf depend on NumPy's CPU dispatch, and the flag recheck re-ran pow, whose flags depend on the math library and CPU as well. - derive the flags in the recheck from inputs and results as C99 (and glibc) raise them, instead of calling pow again - broadcast a scalar exponent with v?PowI and a zero increment instead of v?Powx, which never sets the VM error status (so flags from worker threads were lost) and needed a fix-up for -0.0 and -inf - spell out the expected warnings in the test, and report them on failure --- mkl_umath/src/mkl_umath_loops.c.src | 56 ++++++++++++++--------------- mkl_umath/tests/test_basic.py | 29 +++++++++------ 2 files changed, 46 insertions(+), 39 deletions(-) diff --git a/mkl_umath/src/mkl_umath_loops.c.src b/mkl_umath/src/mkl_umath_loops.c.src index fbbe432e..3fd34b34 100644 --- a/mkl_umath/src/mkl_umath_loops.c.src +++ b/mkl_umath/src/mkl_umath_loops.c.src @@ -112,16 +112,17 @@ } \ } while(0) -#define CHUNKED_VML_POWX_CALL(vml_func, n, type, in1, in2, op1) \ +/* strided VM call with a zero increment for in2, which broadcasts its first element */ +#define CHUNKED_VML_CALL3_S2(vml_func, n, type, in1, in2, op1) \ do { \ npy_intp _n_ = (n); \ const npy_intp _chunk_size = MKL_INT_MAX; \ type *in1p = (type *) (in1); \ - const type _in2 = *((type *) (in2)); \ + type *in2p = (type *) (in2); \ type *op1p = (type *) (op1); \ while (_n_ > 0) { \ npy_intp _current_chunk = (_n_ > _chunk_size) ? _chunk_size : _n_; \ - vml_func((MKL_INT) _current_chunk, in1p, _in2, op1p); \ + vml_func((MKL_INT) _current_chunk, in1p, 1, in2p, 0, op1p, 1); \ _n_ -= _current_chunk; \ in1p += _current_chunk; \ op1p += _current_chunk; \ @@ -859,13 +860,16 @@ mkl_umath_@TYPE@_@kind@(char **args, const npy_intp *dimensions, const npy_intp /**end repeat1**/ /* - * v@s@Pow and v@s@Powx return the same values as pow@c@, except v@s@Powx for - * -0.0 and -inf raised to -0.5, 0.5 and 1.5, which are recomputed. They do - * not raise the same floating-point exception flags though: they raise - * spurious invalid and overflow flags for some special inputs (e.g. 0 ** -1), - * and flags raised in MKL worker threads never reach the calling thread. The - * VM error status does reach it, so whenever MKL reports a problem the flags - * are recomputed with pow@c@. r must not alias x or y, which are read again. + * v@s@Pow returns the same values as pow@c@, a scalar exponent is broadcast + * with v@s@PowI rather than passed to v@s@Powx, which gets -0.0 and -inf + * wrong for some exponents and does not set the VM error status. v@s@Pow does + * not raise the same floating-point exception flags as pow@c@ though: it + * raises spurious invalid and overflow flags for some special inputs (e.g. + * 0 ** -1), and flags raised in MKL worker threads never reach the calling + * thread. The VM error status does reach it, so whenever MKL reports a problem + * the flags are derived from the inputs and results the way C99 (and glibc) + * raise them, which, unlike calling pow@c@ again, does not depend on the math + * library or the CPU. r must not alias x or y, which are read again. */ static void @TYPE@_power_vml(const @type@ *x, const @type@ *y, const int y_is_scalar, @type@ *r, const npy_intp n) @@ -876,36 +880,30 @@ static void feclearexcept(POW_FPE_FLAGS); vmlClearErrStatus(); if (y_is_scalar) { - CHUNKED_VML_POWX_CALL(v@s@Powx, n, @type@, x, y, r); + CHUNKED_VML_CALL3_S2(v@s@PowI, n, @type@, x, y, r); } else { CHUNKED_VML_CALL3(v@s@Pow, n, @type@, @type@, x, y, r); } if (vmlGetErrStatus() != VML_STATUS_OK || fetestexcept(POW_FPE_FLAGS)) { + int flags = 0; npy_intp i; feclearexcept(POW_FPE_FLAGS); for (i = 0; i < n; i++) { - /* volatile keeps the calls, which are only made for their flags */ - volatile @type@ res = pow@c@(x[i], y_is_scalar ? *y : y[i]); - (void) res; - } - } - if (y_is_scalar && (*y == -0.5 || *y == 0.5 || *y == 1.5)) { - const @type@ ninf = -(@type@) NPY_INFINITY; - int candidates = 0; - npy_intp i; + const @type@ in1 = x[i]; + const @type@ in2 = y_is_scalar ? *y : y[i]; - /* branch-free so that it vectorizes, the fix-up below rarely runs */ - for (i = 0; i < n; i++) { - candidates |= (x[i] == 0) | (x[i] == ninf); - } - if (candidates) { - for (i = 0; i < n; i++) { - if ((x[i] == 0 || x[i] == ninf) && signbit(x[i])) { - r[i] = pow@c@(x[i], *y); - } + if (isnan(r[i]) && !isnan(in1) && !isnan(in2)) { + flags |= FE_INVALID; } + else if (isinf(r[i]) && isfinite(in1) && isfinite(in2)) { + /* 0 ** -inf is +inf too, but C99 does not require a flag there */ + flags |= (in1 == 0) ? FE_DIVBYZERO : FE_OVERFLOW; + } + } + if (flags) { + feraiseexcept(flags); } } if (prev_flags) { diff --git a/mkl_umath/tests/test_basic.py b/mkl_umath/tests/test_basic.py index 477dce67..49788c76 100644 --- a/mkl_umath/tests/test_basic.py +++ b/mkl_umath/tests/test_basic.py @@ -221,7 +221,7 @@ def _power_warnings(func, *args): with warnings.catch_warnings(record=True) as w: warnings.simplefilter("always") func(*args) - return sorted(str(x.message) for x in w) + return sorted(str(x.message).split(" encountered")[0] for x in w) @pytest.mark.parametrize("size", [100, 8192 + 1]) @@ -258,17 +258,26 @@ def test_power_inplace(dtype, exponent): @pytest.mark.parametrize("dtype", [np.float32, np.float64]) @pytest.mark.parametrize( - "base, exponent", - [(-1.0, 1.5), (0.0, -1.0), (-0.0, -0.5), (1e30, 3.0), (0.0, -np.inf)], + "base, exponent, expected", + [ + (2.0, 1.5, []), + (-1.0, 1.5, ["invalid value"]), + (0.0, -1.0, ["divide by zero"]), + (-0.0, -0.5, ["divide by zero"]), + (1e30, 20.0, ["overflow"]), + # C99 allows, but does not require, divide by zero here + (0.0, -np.inf, []), + ], ) @pytest.mark.parametrize("pos", [0, -1]) -def test_power_warnings(dtype, base, exponent, pos): - # MKL raises spurious flags for some of these, and flags raised in its - # worker threads (for the element at the end) never reach the caller +def test_power_warnings(dtype, base, exponent, expected, pos): + # The flags C99 requires, as glibc raises them. They are spelled out since + # NumPy's own loops do not raise them the same way on every CPU. MKL raises + # spurious flags for some of these, and flags raised in its worker threads + # (for the element at the end) never reach the caller. x = np.full(10**6, 2.0, dtype=dtype) y = np.full(10**6, 1.5, dtype=dtype) x[pos], y[pos] = base, exponent - assert _power_warnings(mu.power, x, y) == _power_warnings(np.power, x, y) - assert _power_warnings(mu.power, x, y[pos]) == _power_warnings( - np.power, x, y[pos] - ) + for args in ((x, y), (x, y[pos])): + got = _power_warnings(mu.power, *args) + assert got == expected, f"power{(base, exponent)} warned {got}" From 121d3a0138886260f51db02a3139fc0e3b6c008c Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Wed, 14 Oct 2026 17:01:34 +0200 Subject: [PATCH 4/9] Derive power's FP flags from its results on every call On some CI runners (most likely AVX-512) MKL's float32 v?Pow raises neither a flag nor the VM error status for a domain error such as (-1.0) ** 1.5, so the warning was lost. Instead of trusting MKL to report problems, clear its flags after every call, scan the results for nan and inf, and only then derive the C99 flags from inputs and results. --- mkl_umath/src/mkl_umath_loops.c.src | 30 ++++++++++++++++------------- 1 file changed, 17 insertions(+), 13 deletions(-) diff --git a/mkl_umath/src/mkl_umath_loops.c.src b/mkl_umath/src/mkl_umath_loops.c.src index 3fd34b34..1bfc3e61 100644 --- a/mkl_umath/src/mkl_umath_loops.c.src +++ b/mkl_umath/src/mkl_umath_loops.c.src @@ -862,34 +862,38 @@ mkl_umath_@TYPE@_@kind@(char **args, const npy_intp *dimensions, const npy_intp /* * v@s@Pow returns the same values as pow@c@, a scalar exponent is broadcast * with v@s@PowI rather than passed to v@s@Powx, which gets -0.0 and -inf - * wrong for some exponents and does not set the VM error status. v@s@Pow does - * not raise the same floating-point exception flags as pow@c@ though: it - * raises spurious invalid and overflow flags for some special inputs (e.g. - * 0 ** -1), and flags raised in MKL worker threads never reach the calling - * thread. The VM error status does reach it, so whenever MKL reports a problem - * the flags are derived from the inputs and results the way C99 (and glibc) - * raise them, which, unlike calling pow@c@ again, does not depend on the math - * library or the CPU. r must not alias x or y, which are read again. + * wrong for some exponents. The floating-point exception flags MKL raises + * cannot be relied on though: they are spurious for some special inputs + * (e.g. 0 ** -1 also raises invalid), those raised in MKL worker threads never + * reach the calling thread, and on some CPUs float32 domain errors raise + * nothing, not even in the VM error status. So MKL's flags are replaced with + * the ones C99 (and glibc) raise, derived from the inputs and results: each + * such case gives a nan or an inf, which a cheap scan of the results finds. + * r must not alias x or y, which are read again. */ static void @TYPE@_power_vml(const @type@ *x, const @type@ *y, const int y_is_scalar, @type@ *r, const npy_intp n) { /* flags raised by earlier calls of the loop are NumPy's to report */ const int prev_flags = fetestexcept(POW_FPE_FLAGS); + int nonfinite = 0; + npy_intp i; - feclearexcept(POW_FPE_FLAGS); - vmlClearErrStatus(); if (y_is_scalar) { CHUNKED_VML_CALL3_S2(v@s@PowI, n, @type@, x, y, r); } else { CHUNKED_VML_CALL3(v@s@Pow, n, @type@, @type@, x, y, r); } - if (vmlGetErrStatus() != VML_STATUS_OK || fetestexcept(POW_FPE_FLAGS)) { + feclearexcept(POW_FPE_FLAGS); + + /* branch-free so that it vectorizes */ + for (i = 0; i < n; i++) { + nonfinite |= !isfinite(r[i]); + } + if (nonfinite) { int flags = 0; - npy_intp i; - feclearexcept(POW_FPE_FLAGS); for (i = 0; i < n; i++) { const @type@ in1 = x[i]; const @type@ in2 = y_is_scalar ? *y : y[i]; From c494c0a81d103608ab2b26c5de791c8a79012725 Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Wed, 14 Oct 2026 17:15:39 +0200 Subject: [PATCH 5/9] Update the changelog entry for v?PowI --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index fd10468b..182a5ccd 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,7 +7,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [0.6.0] - 2026-10-15 ### Added -* Added oneMKL-backed `power` loops for float32 and float64 arrays, using `v?Pow` and `v?Powx` for a scalar exponent [gh-285](https://github.com/IntelPython/mkl_umath/pull/285) +* Added oneMKL-backed `power` loops for float32 and float64 arrays, using `v?Pow`, and `v?PowI` to broadcast a scalar exponent [gh-285](https://github.com/IntelPython/mkl_umath/pull/285) * Enabled support of Python 3.15 [gh-279](https://github.com/IntelPython/mkl_umath/pull/279) * Added support for free-threaded (GIL-disabled) CPython builds: the Cython extension is compiled with `freethreading_compatible=True` and `_ufuncs` declares `Py_MOD_GIL_NOT_USED`, so importing `mkl_umath` no longer re-enables the GIL [gh-255](https://github.com/IntelPython/mkl_umath/pull/255) From dd828be0cc2e42b5be9023e3e08095dd4e72d0ea Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Wed, 14 Oct 2026 17:54:01 +0200 Subject: [PATCH 6/9] Allow a 1-ULP difference for subnormal power results in the tests On the Windows Jenkins jobs NumPy returns 0 for float32 1e-30 ** 1.5 and 1e30 ** -1.5, where MKL returns the correctly rounded smallest subnormal. A relative tolerance cannot absorb that, so add the dtype's smallest subnormal as absolute tolerance. --- mkl_umath/tests/test_basic.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/mkl_umath/tests/test_basic.py b/mkl_umath/tests/test_basic.py index 49788c76..960ccc10 100644 --- a/mkl_umath/tests/test_basic.py +++ b/mkl_umath/tests/test_basic.py @@ -211,9 +211,14 @@ def test_absolute_nan_signbit(size, dtype): def _assert_power_equal(mkl_res, np_res): - # MKL may differ from libm by 1 ULP, special values must match exactly - eps = np.spacing(np_res.dtype.type(1)) - np.testing.assert_allclose(mkl_res, np_res, rtol=2 * eps, equal_nan=True) + # MKL may differ from libm by 1 ULP, also for subnormal results (NumPy on + # Windows returns 0 for float32 1e-30 ** 1.5), special values must match + # exactly + one = np_res.dtype.type(1) + eps, tiny = np.spacing(one), np.nextafter(0 * one, one) + np.testing.assert_allclose( + mkl_res, np_res, rtol=2 * eps, atol=tiny, equal_nan=True + ) assert (np.signbit(mkl_res) == np.signbit(np_res))[~np.isnan(np_res)].all() From 670f19fca880fbc3dca94165f28485251b34a1b0 Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Thu, 15 Oct 2026 20:04:33 +0200 Subject: [PATCH 7/9] Check the loaded NumPy for power's scalar-exponent fast path The fast path mirrors NumPy >= 2.3, but was gated on the NumPy headers mkl_umath was built against. A build against 2.3+ headers running on NumPy <= 2.2 then gave sqrt semantics for an exponent of 0.5, e.g. (-inf)**0.5 = nan, and a build against older headers running on 2.3+ did not. _ufuncs now passes the runtime C-API feature version to the loops library on import. --- mkl_umath/src/mkl_umath_loops.c.src | 23 +++++++++++++++++++---- mkl_umath/src/mkl_umath_loops.h.src | 8 ++++++++ mkl_umath/src/ufuncsmodule.c | 4 ++++ 3 files changed, 31 insertions(+), 4 deletions(-) diff --git a/mkl_umath/src/mkl_umath_loops.c.src b/mkl_umath/src/mkl_umath_loops.c.src index 1bfc3e61..398825fa 100644 --- a/mkl_umath/src/mkl_umath_loops.c.src +++ b/mkl_umath/src/mkl_umath_loops.c.src @@ -76,6 +76,20 @@ /* block size, in elements, for in-place power, which goes through a buffer */ #define POW_INPLACE_BLOCK ((npy_intp) 1 << 16) +/* not defined by the NumPy < 2.3 headers */ +#ifndef NPY_2_3_API_VERSION +#define NPY_2_3_API_VERSION 0x00000014 +#endif + +/* C-API feature version of the NumPy loaded at runtime, set by _ufuncs on import */ +static int numpy_runtime_version = 0; + +void +mkl_umath_set_numpy_runtime_version(int version) +{ + numpy_runtime_version = version; +} + #define MKL_INT_MAX ((npy_intp) ((~((MKL_UINT) 0)) >> 1)) #define CHUNKED_VML_CALL2(vml_func, n, type, mkl_type, mkl_ftype, in1, op1) \ @@ -920,9 +934,11 @@ mkl_umath_@TYPE@_power(char **args, const npy_intp *dimensions, const npy_intp * { const npy_intp n = dimensions[0]; -#if defined(NPY_2_3_API_VERSION) && NPY_API_VERSION >= NPY_2_3_API_VERSION - /* same fast path for a scalar exponent as in NumPy >= 2.3 */ - if (steps[1] == 0) { + /* + * same fast path for a scalar exponent as in NumPy >= 2.3, checked against + * the NumPy loaded rather than the headers mkl_umath was built against + */ + if (steps[1] == 0 && numpy_runtime_version >= NPY_2_3_API_VERSION) { const @type@ in2 = *(@type@ *)args[1]; if (in2 == -1.0 || in2 == 0.0 || in2 == 0.5 || in2 == 1.0 || in2 == 2.0) { @@ -948,7 +964,6 @@ mkl_umath_@TYPE@_power(char **args, const npy_intp *dimensions, const npy_intp * return; } } -#endif { const int contig = IS_BINARY_CONT(@type@, @type@); diff --git a/mkl_umath/src/mkl_umath_loops.h.src b/mkl_umath/src/mkl_umath_loops.h.src index 4c97ea1e..91d9342e 100644 --- a/mkl_umath/src/mkl_umath_loops.h.src +++ b/mkl_umath/src/mkl_umath_loops.h.src @@ -48,6 +48,14 @@ #define USE_NUMPY_2 #endif +/* + * Records the C-API feature version of the NumPy loaded at runtime, for the + * loops whose results depend on it (power in NumPy >= 2.3). + */ +MKL_UMATH_API +void +mkl_umath_set_numpy_runtime_version(int version); + /**begin repeat * Float types * #TYPE = FLOAT, DOUBLE# diff --git a/mkl_umath/src/ufuncsmodule.c b/mkl_umath/src/ufuncsmodule.c index 12564ab9..7625f355 100644 --- a/mkl_umath/src/ufuncsmodule.c +++ b/mkl_umath/src/ufuncsmodule.c @@ -45,6 +45,10 @@ PyMODINIT_FUNC PyInit__ufuncs(void) import_array(); import_umath(); + /* PyArray_RUNTIME_VERSION, which the NumPy 1.x headers do not provide */ + mkl_umath_set_numpy_runtime_version( + (int)PyArray_GetNDArrayCFeatureVersion()); + m = PyModule_Create(&_ufuncs_module); if (m == NULL) return NULL; From 553f678a4bd808327b510cd46283bd7e38fc347c Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Thu, 15 Oct 2026 20:04:35 +0200 Subject: [PATCH 8/9] Mention v?PowI in the power special-values test comment --- mkl_umath/tests/test_basic.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mkl_umath/tests/test_basic.py b/mkl_umath/tests/test_basic.py index 960ccc10..e13a7e38 100644 --- a/mkl_umath/tests/test_basic.py +++ b/mkl_umath/tests/test_basic.py @@ -232,7 +232,7 @@ def _power_warnings(func, *args): @pytest.mark.parametrize("size", [100, 8192 + 1]) @pytest.mark.parametrize("dtype", [np.float32, np.float64]) def test_power_special_values(size, dtype): - # size 100 takes the fallback loop, 8193 takes MKL's v?Pow and v?Powx + # size 100 takes the fallback loop, 8193 takes MKL's v?Pow and v?PowI specials = np.array(_power_specials, dtype=dtype) x, y = (a.ravel() for a in np.meshgrid(specials, specials)) with np.errstate(all="ignore"): From 2582d50decdb7fe91606cde878f888e1c305af4c Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Thu, 15 Oct 2026 20:04:40 +0200 Subject: [PATCH 9/9] Move the changelog entry under a new [dev] section --- CHANGELOG.md | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 182a5ccd..bd96a346 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,10 +4,14 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). -## [0.6.0] - 2026-10-15 +## [dev] - YYYY-MM-DD ### Added * Added oneMKL-backed `power` loops for float32 and float64 arrays, using `v?Pow`, and `v?PowI` to broadcast a scalar exponent [gh-285](https://github.com/IntelPython/mkl_umath/pull/285) + +## [0.6.0] - 2026-10-15 + +### Added * Enabled support of Python 3.15 [gh-279](https://github.com/IntelPython/mkl_umath/pull/279) * Added support for free-threaded (GIL-disabled) CPython builds: the Cython extension is compiled with `freethreading_compatible=True` and `_ufuncs` declares `Py_MOD_GIL_NOT_USED`, so importing `mkl_umath` no longer re-enables the GIL [gh-255](https://github.com/IntelPython/mkl_umath/pull/255)