From 34edb6e44892910e0b42a38daf3a0f7fc6716666 Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Tue, 19 May 2026 17:10:28 +0200 Subject: [PATCH 01/11] Add dedicated header which defines complex math function to be used in ufunc kernels --- .../elementwise_functions/complex_math.hpp | 106 ++++++++++++++++++ 1 file changed, 106 insertions(+) create mode 100644 dpnp/tensor/libtensor/include/kernels/elementwise_functions/complex_math.hpp diff --git a/dpnp/tensor/libtensor/include/kernels/elementwise_functions/complex_math.hpp b/dpnp/tensor/libtensor/include/kernels/elementwise_functions/complex_math.hpp new file mode 100644 index 00000000000..82dbed83fc3 --- /dev/null +++ b/dpnp/tensor/libtensor/include/kernels/elementwise_functions/complex_math.hpp @@ -0,0 +1,106 @@ +//***************************************************************************** +// Copyright (c) 2026, Intel Corporation +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// - Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// - Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// - Neither the name of the copyright holder nor the names of its contributors +// may be used to endorse or promote products derived from this software +// without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, 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. +//***************************************************************************** +/// +/// \file +/// This file defines complex math functions. +//===----------------------------------------------------------------------===// + +#pragma once + +#include +#include +#include + +#include + +#include "sycl_complex.hpp" // for exprm_ns + +namespace dpnp::tensor::kernels::complex_math +{ +static constexpr double ln2 = 0.6931471805599453094172321214581765L; +static constexpr double pi = 3.1415926535897932384626433832795029L; + +template +T cacos(const T &in) +{ + using realT = typename T::value_type; + using sycl_complexT = exprm_ns::complex; + + static constexpr realT q_nan = std::numeric_limits::quiet_NaN(); + + const realT x = std::real(in); + const realT y = std::imag(in); + + if (std::isnan(x)) { + // acos(NaN + I*+-Inf) = NaN + I*-+Inf + if (std::isinf(y)) { + return T(q_nan, -y); + } + + // all other cases involving NaN return NaN + I*NaN + return T(q_nan, q_nan); + } + + if (std::isnan(y)) { + // acos(+-Inf + I*NaN) = NaN + I*opt(-)Inf + if (std::isinf(x)) { + return T(q_nan, -std::numeric_limits::infinity()); + } + + // acos(0 + I*NaN) = PI/2 + I*NaN with inexact + if (x == realT(0)) { + static constexpr realT pio2 = realT(pi) / realT(2); // PI/2 + return T(pio2, q_nan); + } + + // all other cases involving NaN return NaN + I*NaN + return T(q_nan, q_nan); + } + + /* + * For large x or y including acos(+-Inf + I*+-Inf). + * exprm_ns::acos(x) is based on calculating log(x + sqrt(x^2 - 1)), + * so r_eps = sqrt(1/eps) is appropriate precision loss point. + */ + const realT r_eps = + sycl::sqrt(realT(1) / std::numeric_limits::epsilon()); + if (sycl::fabs(x) > r_eps || sycl::fabs(y) > r_eps) { + sycl_complexT log_in = exprm_ns::log(sycl_complexT(in)); + + const realT wx = log_in.real(); + const realT wy = log_in.imag(); + const realT rx = sycl::fabs(wy); + + realT ry = wx + realT(ln2); + return T(rx, (sycl::signbit(y)) ? ry : -ry); + } + + // ordinary cases + return exprm_ns::acos(sycl_complexT(in)); // sycl::acos(in); +} +} // namespace dpnp::tensor::kernels::complex_math From 70e05c053c721248e86b143dddc9b6ee5ce95859 Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Tue, 19 May 2026 17:16:16 +0200 Subject: [PATCH 02/11] Update acos and acosh to use hew helper dpnp::tensor::kernels::complex_math::cacos() function --- .../kernels/elementwise_functions/acos.hpp | 57 +---------------- .../kernels/elementwise_functions/acosh.hpp | 62 +------------------ 2 files changed, 5 insertions(+), 114 deletions(-) diff --git a/dpnp/tensor/libtensor/include/kernels/elementwise_functions/acos.hpp b/dpnp/tensor/libtensor/include/kernels/elementwise_functions/acos.hpp index 8176d6de62d..4bb165b2940 100644 --- a/dpnp/tensor/libtensor/include/kernels/elementwise_functions/acos.hpp +++ b/dpnp/tensor/libtensor/include/kernels/elementwise_functions/acos.hpp @@ -33,15 +33,13 @@ //===---------------------------------------------------------------------===// #pragma once -#include #include #include #include -#include #include #include -#include "sycl_complex.hpp" +#include "complex_math.hpp" #include "vec_size_util.hpp" #include "kernels/dpnp_tensor_types.hpp" @@ -75,58 +73,7 @@ struct AcosFunctor resT operator()(const argT &in) const { if constexpr (is_complex::value) { - using realT = typename argT::value_type; - - static constexpr realT q_nan = - std::numeric_limits::quiet_NaN(); - - const realT x = std::real(in); - const realT y = std::imag(in); - - if (std::isnan(x)) { - /* acos(NaN + I*+-Inf) = NaN + I*-+Inf */ - if (std::isinf(y)) { - return resT{q_nan, -y}; - } - - /* all other cases involving NaN return NaN + I*NaN. */ - return resT{q_nan, q_nan}; - } - if (std::isnan(y)) { - /* acos(+-Inf + I*NaN) = NaN + I*opt(-)Inf */ - if (std::isinf(x)) { - return resT{q_nan, -std::numeric_limits::infinity()}; - } - /* acos(0 + I*NaN) = PI/2 + I*NaN with inexact */ - if (x == realT(0)) { - const realT res_re = sycl::atan(realT(1)) * 2; // PI/2 - return resT{res_re, q_nan}; - } - - /* all other cases involving NaN return NaN + I*NaN. */ - return resT{q_nan, q_nan}; - } - - /* - * For large x or y including acos(+-Inf + I*+-Inf) - */ - static constexpr realT r_eps = - realT(1) / std::numeric_limits::epsilon(); - if (sycl::fabs(x) > r_eps || sycl::fabs(y) > r_eps) { - using sycl_complexT = exprm_ns::complex; - sycl_complexT log_in = - exprm_ns::log(exprm_ns::complex(in)); - - const realT wx = log_in.real(); - const realT wy = log_in.imag(); - const realT rx = sycl::fabs(wy); - - realT ry = wx + sycl::log(realT(2)); - return resT{rx, (sycl::signbit(y)) ? ry : -ry}; - } - - /* ordinary cases */ - return exprm_ns::acos(exprm_ns::complex(in)); // acos(in); + return dpnp::tensor::kernels::complex_math::cacos(in); } else { static_assert(std::is_floating_point_v || diff --git a/dpnp/tensor/libtensor/include/kernels/elementwise_functions/acosh.hpp b/dpnp/tensor/libtensor/include/kernels/elementwise_functions/acosh.hpp index 559796984b9..5b352077d22 100644 --- a/dpnp/tensor/libtensor/include/kernels/elementwise_functions/acosh.hpp +++ b/dpnp/tensor/libtensor/include/kernels/elementwise_functions/acosh.hpp @@ -37,13 +37,12 @@ #include #include #include -#include #include #include #include -#include "sycl_complex.hpp" +#include "complex_math.hpp" #include "vec_size_util.hpp" #include "kernels/dpnp_tensor_types.hpp" @@ -79,63 +78,8 @@ struct AcoshFunctor if constexpr (is_complex::value) { using realT = typename argT::value_type; - static constexpr realT q_nan = - std::numeric_limits::quiet_NaN(); - /* - * acosh(in) = I*acos(in) or -I*acos(in) - * where the sign is chosen so Re(acosh(in)) >= 0. - * So, we first calculate acos(in) and then acosh(in). - */ - const realT x = std::real(in); - const realT y = std::imag(in); - - resT acos_in; - if (std::isnan(x)) { - /* acos(NaN + I*+-Inf) = NaN + I*-+Inf */ - if (std::isinf(y)) { - acos_in = resT{q_nan, -y}; - } - else { - acos_in = resT{q_nan, q_nan}; - } - } - else if (std::isnan(y)) { - /* acos(+-Inf + I*NaN) = NaN + I*opt(-)Inf */ - static constexpr realT inf = - std::numeric_limits::infinity(); - - if (std::isinf(x)) { - acos_in = resT{q_nan, -inf}; - } - /* acos(0 + I*NaN) = Pi/2 + I*NaN with inexact */ - else if (x == realT(0)) { - const realT pi_half = sycl::atan(realT(1)) * 2; - acos_in = resT{pi_half, q_nan}; - } - else { - acos_in = resT{q_nan, q_nan}; - } - } - - static constexpr realT r_eps = - realT(1) / std::numeric_limits::epsilon(); - /* - * For large x or y including acos(+-Inf + I*+-Inf) - */ - if (sycl::fabs(x) > r_eps || sycl::fabs(y) > r_eps) { - using sycl_complexT = typename exprm_ns::complex; - const sycl_complexT log_in = exprm_ns::log(sycl_complexT(in)); - const realT wx = log_in.real(); - const realT wy = log_in.imag(); - const realT rx = sycl::fabs(wy); - realT ry = wx + sycl::log(realT(2)); - acos_in = resT{rx, (sycl::signbit(y)) ? ry : -ry}; - } - else { - /* ordinary cases */ - acos_in = - exprm_ns::acos(exprm_ns::complex(in)); // acos(in); - } + /* calculate acos value */ + resT acos_in = dpnp::tensor::kernels::complex_math::cacos(in); /* Now we calculate acosh(z) */ const realT rx = std::real(acos_in); From 13d7a9620e2882fbccda6dba1f61375f1f4181c4 Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Wed, 20 May 2026 14:36:42 +0200 Subject: [PATCH 03/11] Add tests for acos and acosh functions --- .../tensor/elementwise/test_hyperbolic.py | 71 +++++++++++++++++++ .../tensor/elementwise/test_trigonometric.py | 49 +++++++++++-- 2 files changed, 115 insertions(+), 5 deletions(-) diff --git a/dpnp/tests/tensor/elementwise/test_hyperbolic.py b/dpnp/tests/tensor/elementwise/test_hyperbolic.py index d25a9845358..62b7da9e2c6 100644 --- a/dpnp/tests/tensor/elementwise/test_hyperbolic.py +++ b/dpnp/tests/tensor/elementwise/test_hyperbolic.py @@ -219,3 +219,74 @@ def test_acosh_zero_nan(dtype): assert np.isnan(Y_dpt.real).all() assert_allclose(np.abs(Y_dpt.imag), np.pi / 2, atol=1e-6, strict=False) + + +@pytest.mark.parametrize("np_call, dpt_call", [(np.arccosh, dpt.acosh)]) +# @pytest.mark.parametrize("np_call, dpt_call", _inv_hyper_funcs) +@pytest.mark.parametrize("dtype", _complex_fp_dtypes) +def test_inv_hyper_large_negative_real(np_call, dpt_call, dtype): + """ + Test inverse hyperbolic functions for large negative real parts. + + Regression acosh test for threshold bug where catastrophic cancellation in + z + sqrt(z² - 1) caused log(0) = -infinity for |Re(z)| in [4e7, 9e7+] + with negative real parts (gh-2924). + """ + + q = get_queue_or_skip() + skip_if_dtype_not_supported(dtype, q) + + # input values that previously returned infinity + thr = np.sqrt(1 / dpt.finfo(dtype).eps) / 2 + x = [ + complex(-4e7, 1.0), # Boundary of bug zone + complex(-9e7, 1.0), # Middle of bug zone + complex(-1e8, 1.0), # Upper range + complex(-4.45712982e8, 1.0), # Original reported value + complex(thr, 1.0), # Exact threshold value + complex(np.nextafter(thr, np.inf), 1.0), # Next after threshold + ] + + xf = np.asarray(x, dtype=dtype) + yf = dpt.asarray(x, dtype=dtype, sycl_queue=q) + + result = dpt_call(yf) + expected = np_call(xf) + + tol = 8 * dpt.finfo(dtype).resolution + assert not dpt.any(dpt.isinf(result)) + assert_allclose(dpt.asnumpy(result), expected, atol=tol, rtol=tol) + + +@pytest.mark.parametrize("np_call, dpt_call", [(np.arccosh, dpt.acosh)]) +# @pytest.mark.parametrize("np_call, dpt_call", _inv_hyper_funcs) +@pytest.mark.parametrize("dtype", _complex_fp_dtypes) +@pytest.mark.parametrize("magnitude", [4e7, 9e7, 1e8, 4.45e8]) +def test_inv_hyper_all_quadrants_large(np_call, dpt_call, dtype, magnitude): + """ + Test inverse hyperbolic functions for large complex values in all quadrants. + + Ensures the threshold fix works correctly for all sign combinations + of real and imaginary parts. + """ + + q = get_queue_or_skip() + skip_if_dtype_not_supported(dtype, q) + + # input values with four quadrants with large magnitude + x = [ + complex(magnitude, 1.0), # +Re, +Im + complex(-magnitude, 1.0), # -Re, +Im (bug zone) + complex(magnitude, -1.0), # +Re, -Im + complex(-magnitude, -1.0), # -Re, -Im (bug zone) + ] + + xf = np.asarray(x, dtype=dtype) + yf = dpt.asarray(x, dtype=dtype, sycl_queue=q) + + result = dpt_call(yf) + expected = np_call(xf) + + tol = 8 * dpt.finfo(dtype).resolution + assert not dpt.any(dpt.isinf(result)) + assert_allclose(dpt.asnumpy(result), expected, atol=tol, rtol=tol) diff --git a/dpnp/tests/tensor/elementwise/test_trigonometric.py b/dpnp/tests/tensor/elementwise/test_trigonometric.py index 49743236030..07715162e2d 100644 --- a/dpnp/tests/tensor/elementwise/test_trigonometric.py +++ b/dpnp/tests/tensor/elementwise/test_trigonometric.py @@ -38,7 +38,9 @@ ) from .utils import ( _all_dtypes, + _complex_fp_dtypes, _map_to_device_dtype, + _real_fp_dtypes, ) _trig_funcs = [(np.sin, dpt.sin), (np.cos, dpt.cos), (np.tan, dpt.tan)] @@ -63,7 +65,7 @@ def test_trig_out_type(np_call, dpt_call, dtype): @pytest.mark.parametrize("np_call, dpt_call", _all_funcs) -@pytest.mark.parametrize("dtype", ["f2", "f4", "f8"]) +@pytest.mark.parametrize("dtype", _real_fp_dtypes) def test_trig_real_contig(np_call, dpt_call, dtype): q = get_queue_or_skip() skip_if_dtype_not_supported(dtype, q) @@ -96,7 +98,7 @@ def test_trig_real_contig(np_call, dpt_call, dtype): @pytest.mark.parametrize("np_call, dpt_call", _all_funcs) -@pytest.mark.parametrize("dtype", ["c8", "c16"]) +@pytest.mark.parametrize("dtype", _complex_fp_dtypes) def test_trig_complex_contig(np_call, dpt_call, dtype): q = get_queue_or_skip() skip_if_dtype_not_supported(dtype, q) @@ -134,7 +136,7 @@ def test_trig_complex_contig(np_call, dpt_call, dtype): @pytest.mark.parametrize("np_call, dpt_call", _all_funcs) -@pytest.mark.parametrize("dtype", ["f2", "f4", "f8"]) +@pytest.mark.parametrize("dtype", _real_fp_dtypes) def test_trig_real_strided(np_call, dpt_call, dtype): q = get_queue_or_skip() skip_if_dtype_not_supported(dtype, q) @@ -168,7 +170,7 @@ def test_trig_real_strided(np_call, dpt_call, dtype): @pytest.mark.parametrize("np_call, dpt_call", _all_funcs) -@pytest.mark.parametrize("dtype", ["c8", "c16"]) +@pytest.mark.parametrize("dtype", _complex_fp_dtypes) def test_trig_complex_strided(np_call, dpt_call, dtype): q = get_queue_or_skip() skip_if_dtype_not_supported(dtype, q) @@ -216,7 +218,7 @@ def test_trig_complex_strided(np_call, dpt_call, dtype): @pytest.mark.parametrize("np_call, dpt_call", _all_funcs) -@pytest.mark.parametrize("dtype", ["f2", "f4", "f8"]) +@pytest.mark.parametrize("dtype", _real_fp_dtypes) def test_trig_real_special_cases(np_call, dpt_call, dtype): q = get_queue_or_skip() skip_if_dtype_not_supported(dtype, q) @@ -232,3 +234,40 @@ def test_trig_real_special_cases(np_call, dpt_call, dtype): tol = 8 * dpt.finfo(dtype).resolution Y = dpt_call(yf) assert_allclose(dpt.asnumpy(Y), Y_np, atol=tol, rtol=tol) + + +@pytest.mark.parametrize("np_call, dpt_call", [(np.arccos, dpt.acos)]) +# @pytest.mark.parametrize("np_call, dpt_call", _inv_trig_funcs) +@pytest.mark.parametrize("dtype", _complex_fp_dtypes) +def test_inv_trig_large_negative_real(np_call, dpt_call, dtype): + """ + Test inverse trigonometric functions for large negative real parts. + + Regression acos test for threshold bug where catastrophic cancellation in + z + sqrt(z² - 1) caused log(0) = -infinity for |Re(z)| in [4e7, 9e7+] + with negative real parts. + """ + + q = get_queue_or_skip() + skip_if_dtype_not_supported(dtype, q) + + # input values that previously returned infinity + thr = np.sqrt(1 / dpt.finfo(dtype).eps) / 2 + x = [ + complex(-4e7, 1.0), # Boundary of bug zone + complex(-9e7, 1.0), # Middle of bug zone + complex(-1e8, 1.0), # Upper range + complex(-4.45712982e8, 1.0), # Original reported value + complex(thr, 1.0), # Exact threshold value + complex(np.nextafter(thr, np.inf), 1.0), # Next after threshold + ] + + xf = np.asarray(x, dtype=dtype) + yf = dpt.asarray(x, dtype=dtype, sycl_queue=q) + + result = dpt_call(yf) + expected = np_call(xf) + + tol = 8 * dpt.finfo(dtype).resolution + assert not dpt.any(dpt.isinf(result)) + assert_allclose(dpt.asnumpy(result), expected, atol=tol, rtol=tol) From 2e06457244b27530ff38289a8f0e605672fe9520 Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Wed, 20 May 2026 14:37:42 +0200 Subject: [PATCH 04/11] Change threashold to sqrt(1/eps)/2 --- .../include/kernels/elementwise_functions/complex_math.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/dpnp/tensor/libtensor/include/kernels/elementwise_functions/complex_math.hpp b/dpnp/tensor/libtensor/include/kernels/elementwise_functions/complex_math.hpp index 82dbed83fc3..ce03b656d04 100644 --- a/dpnp/tensor/libtensor/include/kernels/elementwise_functions/complex_math.hpp +++ b/dpnp/tensor/libtensor/include/kernels/elementwise_functions/complex_math.hpp @@ -85,10 +85,10 @@ T cacos(const T &in) /* * For large x or y including acos(+-Inf + I*+-Inf). * exprm_ns::acos(x) is based on calculating log(x + sqrt(x^2 - 1)), - * so r_eps = sqrt(1/eps) is appropriate precision loss point. + * so r_eps = sqrt(1/eps)/2 is appropriate precision loss point. */ const realT r_eps = - sycl::sqrt(realT(1) / std::numeric_limits::epsilon()); + sycl::sqrt(realT(1) / std::numeric_limits::epsilon()) / 2; if (sycl::fabs(x) > r_eps || sycl::fabs(y) > r_eps) { sycl_complexT log_in = exprm_ns::log(sycl_complexT(in)); From b31c452cc14fe67171df239e4d100fe3e6b30ad9 Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Wed, 20 May 2026 15:34:57 +0200 Subject: [PATCH 05/11] Add casin and casinh functions --- .../elementwise_functions/complex_math.hpp | 94 +++++++++++++++++-- 1 file changed, 86 insertions(+), 8 deletions(-) diff --git a/dpnp/tensor/libtensor/include/kernels/elementwise_functions/complex_math.hpp b/dpnp/tensor/libtensor/include/kernels/elementwise_functions/complex_math.hpp index ce03b656d04..0c4e9797862 100644 --- a/dpnp/tensor/libtensor/include/kernels/elementwise_functions/complex_math.hpp +++ b/dpnp/tensor/libtensor/include/kernels/elementwise_functions/complex_math.hpp @@ -45,41 +45,42 @@ namespace dpnp::tensor::kernels::complex_math static constexpr double ln2 = 0.6931471805599453094172321214581765L; static constexpr double pi = 3.1415926535897932384626433832795029L; +template +static constexpr realT q_nan = std::numeric_limits::quiet_NaN(); + template T cacos(const T &in) { using realT = typename T::value_type; using sycl_complexT = exprm_ns::complex; - static constexpr realT q_nan = std::numeric_limits::quiet_NaN(); - const realT x = std::real(in); const realT y = std::imag(in); if (std::isnan(x)) { // acos(NaN + I*+-Inf) = NaN + I*-+Inf if (std::isinf(y)) { - return T(q_nan, -y); + return T{q_nan, -y}; } // all other cases involving NaN return NaN + I*NaN - return T(q_nan, q_nan); + return T{q_nan, q_nan}; } if (std::isnan(y)) { // acos(+-Inf + I*NaN) = NaN + I*opt(-)Inf if (std::isinf(x)) { - return T(q_nan, -std::numeric_limits::infinity()); + return T{q_nan, -std::numeric_limits::infinity()}; } // acos(0 + I*NaN) = PI/2 + I*NaN with inexact if (x == realT(0)) { static constexpr realT pio2 = realT(pi) / realT(2); // PI/2 - return T(pio2, q_nan); + return T{pio2, q_nan}; } // all other cases involving NaN return NaN + I*NaN - return T(q_nan, q_nan); + return T{q_nan, q_nan}; } /* @@ -97,10 +98,87 @@ T cacos(const T &in) const realT rx = sycl::fabs(wy); realT ry = wx + realT(ln2); - return T(rx, (sycl::signbit(y)) ? ry : -ry); + return T{rx, (sycl::signbit(y)) ? ry : -ry}; } // ordinary cases return exprm_ns::acos(sycl_complexT(in)); // sycl::acos(in); } + +template +T casinh(const T &in) +{ + using realT = typename T::value_type; + using sycl_complexT = exprm_ns::complex; + + const realT x = std::real(in); + const realT y = std::imag(in); + + if (std::isnan(x)) { + // asinh(NaN + I*+-Inf) = opt(+-)Inf + I*NaN + if (std::isinf(y)) { + return T{y, q_nan}; + } + // asinh(NaN + I*0) = NaN + I*0 + if (y == realT(0)) { + return T{q_nan, y}; + } + // all other cases involving NaN return NaN + I*NaN + return T{q_nan, q_nan}; + } + + if (std::isnan(y)) { + // asinh(+-Inf + I*NaN) = +-Inf + I*NaN + if (std::isinf(x)) { + return T{x, q_nan}; + } + // all other cases involving NaN return NaN + I*NaN + return T{q_nan, q_nan}; + } + + /* + * For large x or y including asinh(+-Inf + I*+-Inf) + * asinh(in) = sign(x)*log(sign(x)*in) + O(1/in^2) as in -> infinity + * The above formula works for the imaginary part as well, because + * Im(asinh(in)) = sign(x)*atan2(sign(x)*y, fabs(x)) + O(y/in^3) + * as in -> infinity, uniformly in y. + * + * exprm_ns::asinh(x) is based on calculating log(x + sqrt(x^2 + 1)), + * so r_eps = sqrt(1/eps)/2 is appropriate precision loss point. + */ + const realT r_eps = + sycl::sqrt(realT(1) / std::numeric_limits::epsilon()) / 2; + if (sycl::fabs(x) > r_eps || sycl::fabs(y) > r_eps) { + sycl_complexT log_in = (sycl::signbit(x)) + ? exprm_ns::log(sycl_complexT(-in)) + : exprm_ns::log(sycl_complexT(in)); + realT wx = log_in.real() + realT(ln2); + realT wy = log_in.imag(); + + const realT res_re = sycl::copysign(wx, x); + const realT res_im = sycl::copysign(wy, y); + return T{res_re, res_im}; + } + + // ordinary cases + return exprm_ns::asinh(sycl_complexT(in)); +} + +template +T casin(const T &in) +{ + /* + * casin(z) = reverse(casinh(reverse(z))), + * where reverse(x + I*y) = y + I*x = I*conj(z) + */ + + // reverse(z): swap real and imaginary parts + T reversed{std::imag(in), std::real(in)}; + + // compute asinh of reversed input + T w = casinh(reversed); + + // reverse result back: swap real and imaginary parts + return T{std::imag(w), std::real(w)}; +} } // namespace dpnp::tensor::kernels::complex_math From 2b6717b4739a28951ae55db0362ed72ec8656899 Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Wed, 20 May 2026 15:35:49 +0200 Subject: [PATCH 06/11] Use casin and casinh functions in appropriate kernels --- .../kernels/elementwise_functions/acos.hpp | 2 +- .../kernels/elementwise_functions/acosh.hpp | 2 +- .../kernels/elementwise_functions/asin.hpp | 80 +------------------ .../kernels/elementwise_functions/asinh.hpp | 63 +-------------- 4 files changed, 8 insertions(+), 139 deletions(-) diff --git a/dpnp/tensor/libtensor/include/kernels/elementwise_functions/acos.hpp b/dpnp/tensor/libtensor/include/kernels/elementwise_functions/acos.hpp index 4bb165b2940..7d1de3469a7 100644 --- a/dpnp/tensor/libtensor/include/kernels/elementwise_functions/acos.hpp +++ b/dpnp/tensor/libtensor/include/kernels/elementwise_functions/acos.hpp @@ -39,11 +39,11 @@ #include #include +#include "common.hpp" #include "complex_math.hpp" #include "vec_size_util.hpp" #include "kernels/dpnp_tensor_types.hpp" -#include "kernels/elementwise_functions/common.hpp" #include "utils/type_dispatch_building.hpp" #include "utils/type_utils.hpp" diff --git a/dpnp/tensor/libtensor/include/kernels/elementwise_functions/acosh.hpp b/dpnp/tensor/libtensor/include/kernels/elementwise_functions/acosh.hpp index 5b352077d22..47b50da03d3 100644 --- a/dpnp/tensor/libtensor/include/kernels/elementwise_functions/acosh.hpp +++ b/dpnp/tensor/libtensor/include/kernels/elementwise_functions/acosh.hpp @@ -42,11 +42,11 @@ #include +#include "common.hpp" #include "complex_math.hpp" #include "vec_size_util.hpp" #include "kernels/dpnp_tensor_types.hpp" -#include "kernels/elementwise_functions/common.hpp" #include "utils/type_dispatch_building.hpp" #include "utils/type_utils.hpp" diff --git a/dpnp/tensor/libtensor/include/kernels/elementwise_functions/asin.hpp b/dpnp/tensor/libtensor/include/kernels/elementwise_functions/asin.hpp index 7ab7083498b..ac5b4f52283 100644 --- a/dpnp/tensor/libtensor/include/kernels/elementwise_functions/asin.hpp +++ b/dpnp/tensor/libtensor/include/kernels/elementwise_functions/asin.hpp @@ -33,21 +33,18 @@ //===---------------------------------------------------------------------===// #pragma once -#include -#include #include #include -#include #include #include #include -#include "sycl_complex.hpp" +#include "common.hpp" +#include "complex_math.hpp" #include "vec_size_util.hpp" #include "kernels/dpnp_tensor_types.hpp" -#include "kernels/elementwise_functions/common.hpp" #include "utils/type_dispatch_building.hpp" #include "utils/type_utils.hpp" @@ -77,78 +74,7 @@ struct AsinFunctor resT operator()(const argT &in) const { if constexpr (is_complex::value) { - using realT = typename argT::value_type; - - static constexpr realT q_nan = - std::numeric_limits::quiet_NaN(); - - /* - * asin(in) = I * conj( asinh(I * conj(in)) ) - * so we first calculate w = asinh(I * conj(in)) with - * x = real(I * conj(in)) = imag(in) - * y = imag(I * conj(in)) = real(in) - * and then return {imag(w), real(w)} which is asin(in) - */ - const realT x = std::imag(in); - const realT y = std::real(in); - - if (std::isnan(x)) { - /* asinh(NaN + I*+-Inf) = opt(+-)Inf + I*NaN */ - if (std::isinf(y)) { - const realT asinh_re = y; - const realT asinh_im = q_nan; - return resT{asinh_im, asinh_re}; - } - /* asinh(NaN + I*0) = NaN + I*0 */ - if (y == realT(0)) { - const realT asinh_re = q_nan; - const realT asinh_im = y; - return resT{asinh_im, asinh_re}; - } - /* All other cases involving NaN return NaN + I*NaN. */ - return resT{q_nan, q_nan}; - } - else if (std::isnan(y)) { - /* asinh(+-Inf + I*NaN) = +-Inf + I*NaN */ - if (std::isinf(x)) { - const realT asinh_re = x; - const realT asinh_im = q_nan; - return resT{asinh_im, asinh_re}; - } - /* All other cases involving NaN return NaN + I*NaN. */ - return resT{q_nan, q_nan}; - } - - /* - * For large x or y including asinh(+-Inf + I*+-Inf) - * asinh(in) = sign(x)*log(sign(x)*in) + O(1/in^2) as in -> - * infinity The above formula works for the imaginary part as well, - * because Im(asinh(in)) = sign(x)*atan2(sign(x)*y, fabs(x)) + - * O(y/in^3) as in -> infinity, uniformly in y - */ - static constexpr realT r_eps = - realT(1) / std::numeric_limits::epsilon(); - if (sycl::fabs(x) > r_eps || sycl::fabs(y) > r_eps) { - using sycl_complexT = exprm_ns::complex; - const sycl_complexT z{x, y}; - realT wx, wy; - if (!sycl::signbit(x)) { - const auto log_z = exprm_ns::log(z); - wx = log_z.real() + sycl::log(realT(2)); - wy = log_z.imag(); - } - else { - const auto log_mz = exprm_ns::log(-z); - wx = log_mz.real() + sycl::log(realT(2)); - wy = log_mz.imag(); - } - const realT asinh_re = sycl::copysign(wx, x); - const realT asinh_im = sycl::copysign(wy, y); - return resT{asinh_im, asinh_re}; - } - /* ordinary cases */ - return exprm_ns::asin( - exprm_ns::complex(in)); // sycl::asin(in); + return complex_math::casin(in); } else { static_assert(std::is_floating_point_v || diff --git a/dpnp/tensor/libtensor/include/kernels/elementwise_functions/asinh.hpp b/dpnp/tensor/libtensor/include/kernels/elementwise_functions/asinh.hpp index 1c78ed4f8af..bbb66a6bbf9 100644 --- a/dpnp/tensor/libtensor/include/kernels/elementwise_functions/asinh.hpp +++ b/dpnp/tensor/libtensor/include/kernels/elementwise_functions/asinh.hpp @@ -33,21 +33,18 @@ //===---------------------------------------------------------------------===// #pragma once -#include -#include #include #include -#include #include #include #include -#include "sycl_complex.hpp" +#include "common.hpp" +#include "complex_math.hpp" #include "vec_size_util.hpp" #include "kernels/dpnp_tensor_types.hpp" -#include "kernels/elementwise_functions/common.hpp" #include "utils/type_dispatch_building.hpp" #include "utils/type_utils.hpp" @@ -77,61 +74,7 @@ struct AsinhFunctor resT operator()(const argT &in) const { if constexpr (is_complex::value) { - using realT = typename argT::value_type; - - static constexpr realT q_nan = - std::numeric_limits::quiet_NaN(); - - const realT x = std::real(in); - const realT y = std::imag(in); - - if (std::isnan(x)) { - /* asinh(NaN + I*+-Inf) = opt(+-)Inf + I*NaN */ - if (std::isinf(y)) { - return resT{y, q_nan}; - } - /* asinh(NaN + I*0) = NaN + I*0 */ - if (y == realT(0)) { - return resT{q_nan, y}; - } - /* All other cases involving NaN return NaN + I*NaN. */ - return resT{q_nan, q_nan}; - } - - if (std::isnan(y)) { - /* asinh(+-Inf + I*NaN) = +-Inf + I*NaN */ - if (std::isinf(x)) { - return resT{x, q_nan}; - } - /* All other cases involving NaN return NaN + I*NaN. */ - return resT{q_nan, q_nan}; - } - - /* - * For large x or y including asinh(+-Inf + I*+-Inf) - * asinh(in) = sign(x)*log(sign(x)*in) + O(1/in^2) as in -> - * infinity The above formula works for the imaginary part as well, - * because Im(asinh(in)) = sign(x)*atan2(sign(x)*y, fabs(x)) + - * O(y/in^3) as in -> infinity, uniformly in y - */ - static constexpr realT r_eps = - realT(1) / std::numeric_limits::epsilon(); - - if (sycl::fabs(x) > r_eps || sycl::fabs(y) > r_eps) { - using sycl_complexT = exprm_ns::complex; - sycl_complexT log_in = (sycl::signbit(x)) - ? exprm_ns::log(sycl_complexT(-in)) - : exprm_ns::log(sycl_complexT(in)); - realT wx = log_in.real() + sycl::log(realT(2)); - realT wy = log_in.imag(); - - const realT res_re = sycl::copysign(wx, x); - const realT res_im = sycl::copysign(wy, y); - return resT{res_re, res_im}; - } - - /* ordinary cases */ - return exprm_ns::asinh(exprm_ns::complex(in)); // asinh(in); + return complex_math::casinh(in); } else { static_assert(std::is_floating_point_v || From 66b4d2cc4edcc89e91e016ed760391784ab01921 Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Wed, 20 May 2026 15:36:57 +0200 Subject: [PATCH 07/11] Extend testing to cover asin and asinh --- dpnp/tests/tensor/elementwise/test_hyperbolic.py | 8 ++++++-- dpnp/tests/tensor/elementwise/test_trigonometric.py | 4 +++- 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/dpnp/tests/tensor/elementwise/test_hyperbolic.py b/dpnp/tests/tensor/elementwise/test_hyperbolic.py index 62b7da9e2c6..2b28ec3347a 100644 --- a/dpnp/tests/tensor/elementwise/test_hyperbolic.py +++ b/dpnp/tests/tensor/elementwise/test_hyperbolic.py @@ -221,7 +221,9 @@ def test_acosh_zero_nan(dtype): assert_allclose(np.abs(Y_dpt.imag), np.pi / 2, atol=1e-6, strict=False) -@pytest.mark.parametrize("np_call, dpt_call", [(np.arccosh, dpt.acosh)]) +@pytest.mark.parametrize( + "np_call, dpt_call", [(np.arcsinh, dpt.asinh), (np.arccosh, dpt.acosh)] +) # @pytest.mark.parametrize("np_call, dpt_call", _inv_hyper_funcs) @pytest.mark.parametrize("dtype", _complex_fp_dtypes) def test_inv_hyper_large_negative_real(np_call, dpt_call, dtype): @@ -258,7 +260,9 @@ def test_inv_hyper_large_negative_real(np_call, dpt_call, dtype): assert_allclose(dpt.asnumpy(result), expected, atol=tol, rtol=tol) -@pytest.mark.parametrize("np_call, dpt_call", [(np.arccosh, dpt.acosh)]) +@pytest.mark.parametrize( + "np_call, dpt_call", [(np.arcsinh, dpt.asinh), (np.arccosh, dpt.acosh)] +) # @pytest.mark.parametrize("np_call, dpt_call", _inv_hyper_funcs) @pytest.mark.parametrize("dtype", _complex_fp_dtypes) @pytest.mark.parametrize("magnitude", [4e7, 9e7, 1e8, 4.45e8]) diff --git a/dpnp/tests/tensor/elementwise/test_trigonometric.py b/dpnp/tests/tensor/elementwise/test_trigonometric.py index 07715162e2d..1225557102a 100644 --- a/dpnp/tests/tensor/elementwise/test_trigonometric.py +++ b/dpnp/tests/tensor/elementwise/test_trigonometric.py @@ -236,7 +236,9 @@ def test_trig_real_special_cases(np_call, dpt_call, dtype): assert_allclose(dpt.asnumpy(Y), Y_np, atol=tol, rtol=tol) -@pytest.mark.parametrize("np_call, dpt_call", [(np.arccos, dpt.acos)]) +@pytest.mark.parametrize( + "np_call, dpt_call", [(np.arcsin, dpt.asin), (np.arccos, dpt.acos)] +) # @pytest.mark.parametrize("np_call, dpt_call", _inv_trig_funcs) @pytest.mark.parametrize("dtype", _complex_fp_dtypes) def test_inv_trig_large_negative_real(np_call, dpt_call, dtype): From 997e3ec63f7c74f52d4facbf28bb0d904c1768bb Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Wed, 20 May 2026 16:08:09 +0200 Subject: [PATCH 08/11] Add catan and catanh functions --- .../elementwise_functions/complex_math.hpp | 74 +++++++++++++++++++ 1 file changed, 74 insertions(+) diff --git a/dpnp/tensor/libtensor/include/kernels/elementwise_functions/complex_math.hpp b/dpnp/tensor/libtensor/include/kernels/elementwise_functions/complex_math.hpp index 0c4e9797862..a820e32e328 100644 --- a/dpnp/tensor/libtensor/include/kernels/elementwise_functions/complex_math.hpp +++ b/dpnp/tensor/libtensor/include/kernels/elementwise_functions/complex_math.hpp @@ -181,4 +181,78 @@ T casin(const T &in) // reverse result back: swap real and imaginary parts return T{std::imag(w), std::real(w)}; } + +template +T catanh(const T &in) +{ + using realT = typename T::value_type; + + static constexpr realT pio2 = realT(pi) / realT(2); // PI/2 + + const realT x = std::real(in); + const realT y = std::imag(in); + + if (std::isnan(x)) { + // atanh(NaN + I*+-Inf) = sign(NaN)*0 + I*+-PI/2 + if (std::isinf(y)) { + const realT res_re = sycl::copysign(realT(0), x); + const realT res_im = sycl::copysign(pio2, y); + return T{res_re, res_im}; + } + + // all other cases involving NaN return NaN + I*NaN + return T{q_nan, q_nan}; + } + + if (std::isnan(y)) { + // atanh(+-Inf + I*NaN) = +-0 + I*NaN + if (std::isinf(x)) { + const realT res_re = sycl::copysign(realT(0), x); + return T{res_re, q_nan}; + } + + // atanh(+-0 + I*NaN) = +-0 + I*NaN + if (x == realT(0)) { + return T{x, q_nan}; + } + + // all other cases involving NaN return NaN + I*NaN + return T{q_nan, q_nan}; + } + + /* + * For large x or y including atanh(+-Inf + I*+-Inf) = 0 + I*+-PI/2 + * The sign of PI/2 depends on the sign of imaginary part of the input. + * + * exprm_ns::atanh(x) is based on calculating log((1 + x) / (1 - x)) / 2, + * so r_eps = 1/eps is appropriate precision loss point. + */ + const realT r_eps = realT(1) / std::numeric_limits::epsilon(); + if (sycl::fabs(x) > r_eps || sycl::fabs(y) > r_eps) { + const realT res_re = realT(0); + const realT res_im = sycl::copysign(pio2, y); + return T{res_re, res_im}; + } + + // ordinary cases + return exprm_ns::atanh(exprm_ns::complex(in)); +} + +template +T catan(const T &in) +{ + /* + * catan(z) = reverse(catanh(reverse(z))), + * where reverse(x + I*y) = y + I*x = I*conj(z) + */ + + // reverse(z): swap real and imaginary parts + T reversed{std::imag(in), std::real(in)}; + + // compute atanh of reversed input + T w = catanh(reversed); + + // reverse result back: swap real and imaginary parts + return T{std::imag(w), std::real(w)}; +} } // namespace dpnp::tensor::kernels::complex_math From c8173769c480dee89d8102dc03db63bafe32b7f6 Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Wed, 20 May 2026 16:08:48 +0200 Subject: [PATCH 09/11] Use catan and catanh functions in appropriate kernels --- .../kernels/elementwise_functions/atan.hpp | 69 +------------------ .../kernels/elementwise_functions/atanh.hpp | 63 +---------------- 2 files changed, 6 insertions(+), 126 deletions(-) diff --git a/dpnp/tensor/libtensor/include/kernels/elementwise_functions/atan.hpp b/dpnp/tensor/libtensor/include/kernels/elementwise_functions/atan.hpp index 2cef9ba11e1..5ff9e342c77 100644 --- a/dpnp/tensor/libtensor/include/kernels/elementwise_functions/atan.hpp +++ b/dpnp/tensor/libtensor/include/kernels/elementwise_functions/atan.hpp @@ -33,21 +33,18 @@ //===---------------------------------------------------------------------===// #pragma once -#include -#include #include #include -#include #include #include #include -#include "sycl_complex.hpp" +#include "common.hpp" +#include "complex_math.hpp" #include "vec_size_util.hpp" #include "kernels/dpnp_tensor_types.hpp" -#include "kernels/elementwise_functions/common.hpp" #include "utils/type_dispatch_building.hpp" #include "utils/type_utils.hpp" @@ -80,67 +77,7 @@ struct AtanFunctor resT operator()(const argT &in) const { if constexpr (is_complex::value) { - using realT = typename argT::value_type; - - static constexpr realT q_nan = - std::numeric_limits::quiet_NaN(); - /* - * atan(in) = I * conj( atanh(I * conj(in)) ) - * so we first calculate w = atanh(I * conj(in)) with - * x = real(I * conj(in)) = imag(in) - * y = imag(I * conj(in)) = real(in) - * and then return {imag(w), real(w)} which is atan(in) - */ - const realT x = std::imag(in); - const realT y = std::real(in); - if (std::isnan(x)) { - /* atanh(NaN + I*+-Inf) = sign(NaN)*0 + I*+-Pi/2 */ - if (std::isinf(y)) { - const realT pi_half = sycl::atan(realT(1)) * 2; - - const realT atanh_re = sycl::copysign(realT(0), x); - const realT atanh_im = sycl::copysign(pi_half, y); - return resT{atanh_im, atanh_re}; - } - /* - * All other cases involving NaN return NaN + I*NaN. - */ - return resT{q_nan, q_nan}; - } - else if (std::isnan(y)) { - /* atanh(+-Inf + I*NaN) = +-0 + I*NaN */ - if (std::isinf(x)) { - const realT atanh_re = sycl::copysign(realT(0), x); - const realT atanh_im = q_nan; - return resT{atanh_im, atanh_re}; - } - /* atanh(+-0 + I*NaN) = +-0 + I*NaN */ - if (x == realT(0)) { - return resT{q_nan, x}; - } - /* - * All other cases involving NaN return NaN + I*NaN. - */ - return resT{q_nan, q_nan}; - } - - /* - * For large x or y including - * atanh(+-Inf + I*+-Inf) = 0 + I*+-PI/2 - * The sign of pi/2 depends on the sign of imaginary part of the - * input. - */ - static constexpr realT r_eps = - realT(1) / std::numeric_limits::epsilon(); - if (sycl::fabs(x) > r_eps || sycl::fabs(y) > r_eps) { - const realT pi_half = sycl::atan(realT(1)) * 2; - - const realT atanh_re = realT(0); - const realT atanh_im = sycl::copysign(pi_half, y); - return resT{atanh_im, atanh_re}; - } - /* ordinary cases */ - return exprm_ns::atan(exprm_ns::complex(in)); // atan(in); + return complex_math::catan(in); } else { static_assert(std::is_floating_point_v || diff --git a/dpnp/tensor/libtensor/include/kernels/elementwise_functions/atanh.hpp b/dpnp/tensor/libtensor/include/kernels/elementwise_functions/atanh.hpp index 7c74dc95cae..e349053a3f7 100644 --- a/dpnp/tensor/libtensor/include/kernels/elementwise_functions/atanh.hpp +++ b/dpnp/tensor/libtensor/include/kernels/elementwise_functions/atanh.hpp @@ -34,21 +34,18 @@ #pragma once -#include -#include #include #include -#include #include #include #include -#include "sycl_complex.hpp" +#include "common.hpp" +#include "complex_math.hpp" #include "vec_size_util.hpp" #include "kernels/dpnp_tensor_types.hpp" -#include "kernels/elementwise_functions/common.hpp" #include "utils/type_dispatch_building.hpp" #include "utils/type_utils.hpp" @@ -78,61 +75,7 @@ struct AtanhFunctor resT operator()(const argT &in) const { if constexpr (is_complex::value) { - using realT = typename argT::value_type; - static constexpr realT q_nan = - std::numeric_limits::quiet_NaN(); - - const realT x = std::real(in); - const realT y = std::imag(in); - - if (std::isnan(x)) { - /* atanh(NaN + I*+-Inf) = sign(NaN)0 + I*+-PI/2 */ - if (std::isinf(y)) { - const realT pi_half = sycl::atan(realT(1)) * 2; - - const realT res_re = sycl::copysign(realT(0), x); - const realT res_im = sycl::copysign(pi_half, y); - return resT{res_re, res_im}; - } - /* - * All other cases involving NaN return NaN + I*NaN. - */ - return resT{q_nan, q_nan}; - } - else if (std::isnan(y)) { - /* atanh(+-Inf + I*NaN) = +-0 + I*NaN */ - if (std::isinf(x)) { - const realT res_re = sycl::copysign(realT(0), x); - return resT{res_re, q_nan}; - } - /* atanh(+-0 + I*NaN) = +-0 + I*NaN */ - if (x == realT(0)) { - return resT{x, q_nan}; - } - /* - * All other cases involving NaN return NaN + I*NaN. - */ - return resT{q_nan, q_nan}; - } - - /* - * For large x or y including - * atanh(+-Inf + I*+-Inf) = 0 + I*+-PI/2 - * The sign of PI/2 depends on the sign of imaginary part of the - * input. - */ - const realT RECIP_EPSILON = - realT(1) / std::numeric_limits::epsilon(); - if (sycl::fabs(x) > RECIP_EPSILON || - sycl::fabs(y) > RECIP_EPSILON) { - const realT pi_half = sycl::atan(realT(1)) * 2; - - const realT res_re = realT(0); - const realT res_im = sycl::copysign(pi_half, y); - return resT{res_re, res_im}; - } - /* ordinary cases */ - return exprm_ns::atanh(exprm_ns::complex(in)); // atanh(in); + return complex_math::catanh(in); } else { static_assert(std::is_floating_point_v || From 2b8414f0452f2387291db35116c09c69c4d1bb52 Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Wed, 20 May 2026 16:09:48 +0200 Subject: [PATCH 10/11] Extend testing to cover atan and atanh --- .../tensor/elementwise/test_hyperbolic.py | 37 +++++++++++-------- .../tensor/elementwise/test_trigonometric.py | 16 ++++---- 2 files changed, 31 insertions(+), 22 deletions(-) diff --git a/dpnp/tests/tensor/elementwise/test_hyperbolic.py b/dpnp/tests/tensor/elementwise/test_hyperbolic.py index 2b28ec3347a..16a4ca2439b 100644 --- a/dpnp/tests/tensor/elementwise/test_hyperbolic.py +++ b/dpnp/tests/tensor/elementwise/test_hyperbolic.py @@ -221,10 +221,7 @@ def test_acosh_zero_nan(dtype): assert_allclose(np.abs(Y_dpt.imag), np.pi / 2, atol=1e-6, strict=False) -@pytest.mark.parametrize( - "np_call, dpt_call", [(np.arcsinh, dpt.asinh), (np.arccosh, dpt.acosh)] -) -# @pytest.mark.parametrize("np_call, dpt_call", _inv_hyper_funcs) +@pytest.mark.parametrize("np_call, dpt_call", _inv_hyper_funcs) @pytest.mark.parametrize("dtype", _complex_fp_dtypes) def test_inv_hyper_large_negative_real(np_call, dpt_call, dtype): """ @@ -239,14 +236,19 @@ def test_inv_hyper_large_negative_real(np_call, dpt_call, dtype): skip_if_dtype_not_supported(dtype, q) # input values that previously returned infinity - thr = np.sqrt(1 / dpt.finfo(dtype).eps) / 2 + thr1 = 1 / dpt.finfo(dtype).eps # acosh, asinh + thr2 = np.sqrt(thr1) / 2 # atanh x = [ complex(-4e7, 1.0), # Boundary of bug zone complex(-9e7, 1.0), # Middle of bug zone complex(-1e8, 1.0), # Upper range complex(-4.45712982e8, 1.0), # Original reported value - complex(thr, 1.0), # Exact threshold value - complex(np.nextafter(thr, np.inf), 1.0), # Next after threshold + # Exact threshold values + complex(thr1, 1.0), + complex(thr2, 1.0), + # Next values after threshold + complex(np.nextafter(thr1, np.inf), 1.0), + complex(np.nextafter(thr2, np.inf), 1.0), ] xf = np.asarray(x, dtype=dtype) @@ -260,10 +262,7 @@ def test_inv_hyper_large_negative_real(np_call, dpt_call, dtype): assert_allclose(dpt.asnumpy(result), expected, atol=tol, rtol=tol) -@pytest.mark.parametrize( - "np_call, dpt_call", [(np.arcsinh, dpt.asinh), (np.arccosh, dpt.acosh)] -) -# @pytest.mark.parametrize("np_call, dpt_call", _inv_hyper_funcs) +@pytest.mark.parametrize("np_call, dpt_call", _inv_hyper_funcs) @pytest.mark.parametrize("dtype", _complex_fp_dtypes) @pytest.mark.parametrize("magnitude", [4e7, 9e7, 1e8, 4.45e8]) def test_inv_hyper_all_quadrants_large(np_call, dpt_call, dtype, magnitude): @@ -278,11 +277,19 @@ def test_inv_hyper_all_quadrants_large(np_call, dpt_call, dtype, magnitude): skip_if_dtype_not_supported(dtype, q) # input values with four quadrants with large magnitude + thr1 = 1 / dpt.finfo(dtype).eps # acosh, asinh + thr2 = np.sqrt(thr1) / 2 # atanh x = [ - complex(magnitude, 1.0), # +Re, +Im - complex(-magnitude, 1.0), # -Re, +Im (bug zone) - complex(magnitude, -1.0), # +Re, -Im - complex(-magnitude, -1.0), # -Re, -Im (bug zone) + complex(-4e7, 1.0), # Boundary of bug zone + complex(-9e7, 1.0), # Middle of bug zone + complex(-1e8, 1.0), # Upper range + complex(-4.45712982e8, 1.0), # Original reported value + # Exact threshold values + complex(thr1, 1.0), + complex(thr2, 1.0), + # Next values after threshold + complex(np.nextafter(thr1, np.inf), 1.0), + complex(np.nextafter(thr2, np.inf), 1.0), ] xf = np.asarray(x, dtype=dtype) diff --git a/dpnp/tests/tensor/elementwise/test_trigonometric.py b/dpnp/tests/tensor/elementwise/test_trigonometric.py index 1225557102a..0c262f32200 100644 --- a/dpnp/tests/tensor/elementwise/test_trigonometric.py +++ b/dpnp/tests/tensor/elementwise/test_trigonometric.py @@ -236,10 +236,7 @@ def test_trig_real_special_cases(np_call, dpt_call, dtype): assert_allclose(dpt.asnumpy(Y), Y_np, atol=tol, rtol=tol) -@pytest.mark.parametrize( - "np_call, dpt_call", [(np.arcsin, dpt.asin), (np.arccos, dpt.acos)] -) -# @pytest.mark.parametrize("np_call, dpt_call", _inv_trig_funcs) +@pytest.mark.parametrize("np_call, dpt_call", _inv_trig_funcs) @pytest.mark.parametrize("dtype", _complex_fp_dtypes) def test_inv_trig_large_negative_real(np_call, dpt_call, dtype): """ @@ -254,14 +251,19 @@ def test_inv_trig_large_negative_real(np_call, dpt_call, dtype): skip_if_dtype_not_supported(dtype, q) # input values that previously returned infinity - thr = np.sqrt(1 / dpt.finfo(dtype).eps) / 2 + thr1 = 1 / dpt.finfo(dtype).eps # acos, asin + thr2 = np.sqrt(thr1) / 2 # atan x = [ complex(-4e7, 1.0), # Boundary of bug zone complex(-9e7, 1.0), # Middle of bug zone complex(-1e8, 1.0), # Upper range complex(-4.45712982e8, 1.0), # Original reported value - complex(thr, 1.0), # Exact threshold value - complex(np.nextafter(thr, np.inf), 1.0), # Next after threshold + # Exact threshold values + complex(thr1, 1.0), + complex(thr2, 1.0), + # Next values after threshold + complex(np.nextafter(thr1, np.inf), 1.0), + complex(np.nextafter(thr2, np.inf), 1.0), ] xf = np.asarray(x, dtype=dtype) From 00851825d4364c4762062e1a96a16848f0f67e3e Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Wed, 20 May 2026 17:52:18 +0200 Subject: [PATCH 11/11] Update changelog --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2a738267346..bffdced21f6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -30,6 +30,7 @@ This release is compatible with NumPy 2.4.5. * Fixed incorrect `dpnp.tensor.acosh` result for `complex(±0, NaN)` special case to match the Python Array API specification [#2914](https://github.com/IntelPython/dpnp/pull/2914) * Fixed fork PR documentation workflow failures by implementing conditional publishing strategy: upstream PRs publish to GitHub Pages with comment, fork PRs upload artifacts [#2910](https://github.com/IntelPython/dpnp/pull/2910) * Fixed missing `libtensor` headers in the installed `dpnp` package [#2915](https://github.com/IntelPython/dpnp/pull/2915) +* Fixed `dpnp.tensor.acosh` and `dpnp.tensor.acos` returning infinity for complex numbers with large negative real parts [#2928](https://github.com/IntelPython/dpnp/pull/2928) ### Security