Add complex type handling

Author: antonwolfy

dpnp/backend/extensions/ufunc/elementwise_functions/sinc.cpp

@@ -59,11 +59,13 @@ namespace td_ns = dpctl::tensor::type_dispatch;
 template <typename T>
 struct OutputType
 {
-    using value_type =
-        typename std::disjunction<td_ns::TypeMapResultEntry<T, sycl::half>,
-                                  td_ns::TypeMapResultEntry<T, float>,
-                                  td_ns::TypeMapResultEntry<T, double>,
-                                  td_ns::DefaultResultEntry<void>>::result_type;
+    using value_type = typename std::disjunction<
+        td_ns::TypeMapResultEntry<T, sycl::half>,
+        td_ns::TypeMapResultEntry<T, float>,
+        td_ns::TypeMapResultEntry<T, double>,
+        td_ns::TypeMapResultEntry<T, std::complex<float>>,
+        td_ns::TypeMapResultEntry<T, std::complex<double>>,
+        td_ns::DefaultResultEntry<void>>::result_type;
 };
 
 using dpnp::kernels::sinc::SincFunctor;

dpnp/backend/kernels/elementwise_functions/sinc.hpp

@@ -25,32 +25,176 @@
 
 #pragma once
 
+#define SYCL_EXT_ONEAPI_COMPLEX
+#if __has_include(<sycl/ext/oneapi/experimental/sycl_complex.hpp>)
+#include <sycl/ext/oneapi/experimental/sycl_complex.hpp>
+#else
+#include <sycl/ext/oneapi/experimental/complex/complex.hpp>
+#endif
+
 #include <sycl/sycl.hpp>
 
+// dpctl tensor headers
+#include "utils/type_utils.hpp"
+
 namespace dpnp::kernels::sinc
 {
-template <typename argT, typename resT>
+namespace tu_ns = dpctl::tensor::type_utils;
+
+namespace impl
+{
+namespace exprm_ns = sycl::ext::oneapi::experimental;
+
+template <typename Tp>
+inline Tp sin(const Tp &in)
+{
+    if constexpr (tu_ns::is_complex<Tp>::value) {
+        using realTp = typename Tp::value_type;
+
+        constexpr realTp q_nan = std::numeric_limits<realTp>::quiet_NaN();
+
+        realTp const &in_re = std::real(in);
+        realTp const &in_im = std::imag(in);
+
+        const bool in_re_finite = sycl::isfinite(in_re);
+        const bool in_im_finite = sycl::isfinite(in_im);
+        /*
+         * Handle the nearly-non-exceptional cases where
+         * real and imaginary parts of input are finite.
+         */
+        if (in_re_finite && in_im_finite) {
+            Tp res = exprm_ns::sin(exprm_ns::complex<realTp>(in)); // sin(in);
+            if (in_re == realTp(0)) {
+                res.real(sycl::copysign(realTp(0), in_re));
+            }
+            return res;
+        }
+
+        /*
+         * since sin(in) = -I * sinh(I * in), for special cases,
+         * we calculate real and imaginary parts of z = sinh(I * in) and
+         * then return { imag(z) , -real(z) } which is sin(in).
+         */
+        const realTp x = -in_im;
+        const realTp y = in_re;
+        const bool xfinite = in_im_finite;
+        const bool yfinite = in_re_finite;
+        /*
+         * sinh(+-0 +- I Inf) = sign(d(+-0, dNaN))0 + I dNaN.
+         * The sign of 0 in the result is unspecified.  Choice = normally
+         * the same as dNaN.
+         *
+         * sinh(+-0 +- I NaN) = sign(d(+-0, NaN))0 + I d(NaN).
+         * The sign of 0 in the result is unspecified.  Choice = normally
+         * the same as d(NaN).
+         */
+        if (x == realTp(0) && !yfinite) {
+            const realTp sinh_im = q_nan;
+            const realTp sinh_re = sycl::copysign(realTp(0), x * sinh_im);
+            return Tp{sinh_im, -sinh_re};
+        }
+
+        /*
+         * sinh(+-Inf +- I 0) = +-Inf + I +-0.
+         *
+         * sinh(NaN +- I 0)   = d(NaN) + I +-0.
+         */
+        if (y == realTp(0) && !xfinite) {
+            if (std::isnan(x)) {
+                const realTp sinh_re = x;
+                const realTp sinh_im = y;
+                return Tp{sinh_im, -sinh_re};
+            }
+            const realTp sinh_re = x;
+            const realTp sinh_im = sycl::copysign(realTp(0), y);
+            return Tp{sinh_im, -sinh_re};
+        }
+
+        /*
+         * sinh(x +- I Inf) = dNaN + I dNaN.
+         *
+         * sinh(x + I NaN) = d(NaN) + I d(NaN).
+         */
+        if (xfinite && !yfinite) {
+            const realTp sinh_re = q_nan;
+            const realTp sinh_im = x * sinh_re;
+            return Tp{sinh_im, -sinh_re};
+        }
+
+        /*
+         * sinh(+-Inf + I NaN)  = +-Inf + I d(NaN).
+         * The sign of Inf in the result is unspecified.  Choice = normally
+         * the same as d(NaN).
+         *
+         * sinh(+-Inf +- I Inf) = +Inf + I dNaN.
+         * The sign of Inf in the result is unspecified.
+         * Choice = always - here for sinh to have positive result for
+         * imaginary part of sin.
+         *
+         * sinh(+-Inf + I y)   = +-Inf cos(y) + I Inf sin(y)
+         */
+        if (std::isinf(x)) {
+            if (!yfinite) {
+                const realTp sinh_re = -x * x;
+                const realTp sinh_im = x * (y - y);
+                return Tp{sinh_im, -sinh_re};
+            }
+            const realTp sinh_re = x * sycl::cos(y);
+            const realTp sinh_im =
+                std::numeric_limits<realTp>::infinity() * sycl::sin(y);
+            return Tp{sinh_im, -sinh_re};
+        }
+
+        /*
+         * sinh(NaN + I NaN)  = d(NaN) + I d(NaN).
+         *
+         * sinh(NaN +- I Inf) = d(NaN) + I d(NaN).
+         *
+         * sinh(NaN + I y)    = d(NaN) + I d(NaN).
+         */
+        const realTp y_m_y = (y - y);
+        const realTp sinh_re = (x * x) * y_m_y;
+        const realTp sinh_im = (x + x) * y_m_y;
+        return Tp{sinh_im, -sinh_re};
+    }
+    else {
+        if (in == Tp(0)) {
+            return in;
+        }
+        return sycl::sin(in);
+    }
+}
+} // namespace impl
+
+template <typename argT, typename Tp>
 struct SincFunctor
 {
     // is function constant for given argT
     using is_constant = typename std::false_type;
     // constant value, if constant
-    // constexpr resT constant_value = resT{};
+    // constexpr Tp constant_value = Tp{};
     // is function defined for sycl::vec
     using supports_vec = typename std::false_type;
-    // do both argT and resT support subgroup store/load operation
-    using supports_sg_loadstore = typename std::true_type;
+    // do both argT and Tp support subgroup store/load operation
+    using supports_sg_loadstore = typename std::negation<
+        std::disjunction<tu_ns::is_complex<Tp>, tu_ns::is_complex<argT>>>;
 
-    resT operator()(const argT &x) const
+    Tp operator()(const argT &x) const
     {
         constexpr argT pi =
             static_cast<argT>(3.1415926535897932384626433832795029L);
         const argT y = pi * x;
 
         if (y == argT(0)) {
-            return argT(1);
+            return Tp(1);
+        }
+
+        if constexpr (tu_ns::is_complex<argT>::value) {
+            return impl::sin(y) / y;
+        }
+        else {
+            return sycl::sinpi(x) / y;
         }
-        return sycl::sinpi(x) / y;
     }
 };
 } // namespace dpnp::kernels::sinc