diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/abs.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/abs.hpp
index 9c42ba1bb4..4327a9f3ff 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/abs.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/abs.hpp
@@ -33,6 +33,7 @@
 
 #include "kernels/elementwise_functions/common.hpp"
 
+#include "utils/math_utils.hpp"
 #include "utils/offset_utils.hpp"
 #include "utils/type_dispatch.hpp"
 #include "utils/type_utils.hpp"
@@ -49,6 +50,7 @@ namespace abs
 
 namespace py = pybind11;
 namespace td_ns = dpctl::tensor::type_dispatch;
+namespace mu_ns = dpctl::tensor::math_utils;
 
 using dpctl::tensor::type_utils::is_complex;
 
@@ -106,16 +108,16 @@ template <typename argT, typename resT> struct AbsFunctor
         constexpr realT q_nan = std::numeric_limits<realT>::quiet_NaN();
         constexpr realT p_inf = std::numeric_limits<realT>::infinity();
 
-        if (std::isinf(x)) {
+        if (mu_ns::isinf(x)) {
             return p_inf;
         }
-        else if (std::isinf(y)) {
+        else if (mu_ns::isinf(y)) {
             return p_inf;
         }
-        else if (std::isnan(x)) {
+        else if (mu_ns::isnan(x)) {
             return q_nan;
         }
-        else if (std::isnan(y)) {
+        else if (mu_ns::isnan(y)) {
             return q_nan;
         }
         else {
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/acos.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/acos.hpp
index 5dcd1268b5..aaff3ee4c8 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/acos.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/acos.hpp
@@ -31,6 +31,7 @@
 
 #include "kernels/elementwise_functions/common.hpp"
 
+#include "utils/math_utils.hpp"
 #include "utils/offset_utils.hpp"
 #include "utils/type_dispatch.hpp"
 #include "utils/type_utils.hpp"
@@ -47,6 +48,7 @@ namespace acos
 
 namespace py = pybind11;
 namespace td_ns = dpctl::tensor::type_dispatch;
+namespace mu_ns = dpctl::tensor::math_utils;
 
 using dpctl::tensor::type_utils::is_complex;
 
@@ -73,18 +75,18 @@ template <typename argT, typename resT> struct AcosFunctor
             const realT x = std::real(in);
             const realT y = std::imag(in);
 
-            if (std::isnan(x)) {
+            if (mu_ns::isnan(x)) {
                 /* acos(NaN + I*+-Inf) = NaN + I*-+Inf */
-                if (std::isinf(y)) {
+                if (mu_ns::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)) {
+            if (mu_ns::isnan(y)) {
                 /* acos(+-Inf + I*NaN) = NaN + I*opt(-)Inf */
-                if (std::isinf(x)) {
+                if (mu_ns::isinf(x)) {
                     return resT{q_nan, -std::numeric_limits<realT>::infinity()};
                 }
                 /* acos(0 + I*NaN) = PI/2 + I*NaN with inexact */
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/acosh.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/acosh.hpp
index 1dca4d36d6..ed16372dba 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/acosh.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/acosh.hpp
@@ -31,6 +31,7 @@
 
 #include "kernels/elementwise_functions/common.hpp"
 
+#include "utils/math_utils.hpp"
 #include "utils/offset_utils.hpp"
 #include "utils/type_dispatch.hpp"
 #include "utils/type_utils.hpp"
@@ -47,6 +48,7 @@ namespace acosh
 
 namespace py = pybind11;
 namespace td_ns = dpctl::tensor::type_dispatch;
+namespace mu_ns = dpctl::tensor::math_utils;
 
 using dpctl::tensor::type_utils::is_complex;
 
@@ -78,20 +80,20 @@ template <typename argT, typename resT> struct AcoshFunctor
             const realT y = std::imag(in);
 
             resT acos_in;
-            if (std::isnan(x)) {
+            if (mu_ns::isnan(x)) {
                 /* acos(NaN + I*+-Inf) = NaN + I*-+Inf */
-                if (std::isinf(y)) {
+                if (mu_ns::isinf(y)) {
                     acos_in = resT{q_nan, -y};
                 }
                 else {
                     acos_in = resT{q_nan, q_nan};
                 }
             }
-            else if (std::isnan(y)) {
+            else if (mu_ns::isnan(y)) {
                 /* acos(+-Inf + I*NaN) = NaN + I*opt(-)Inf */
                 constexpr realT inf = std::numeric_limits<realT>::infinity();
 
-                if (std::isinf(x)) {
+                if (mu_ns::isinf(x)) {
                     acos_in = resT{q_nan, -inf};
                 }
                 /* acos(0 + I*NaN) = Pi/2 + I*NaN with inexact */
@@ -126,16 +128,16 @@ template <typename argT, typename resT> struct AcoshFunctor
             const realT ry = std::imag(acos_in);
 
             /* acosh(NaN + I*NaN) = NaN + I*NaN */
-            if (std::isnan(rx) && std::isnan(ry)) {
+            if (mu_ns::isnan(rx) && mu_ns::isnan(ry)) {
                 return resT{ry, rx};
             }
             /* acosh(NaN + I*+-Inf) = +Inf + I*NaN */
             /* acosh(+-Inf + I*NaN) = +Inf + I*NaN */
-            if (std::isnan(rx)) {
+            if (mu_ns::isnan(rx)) {
                 return resT{std::abs(ry), rx};
             }
             /* acosh(0 + I*NaN) = NaN + I*NaN */
-            if (std::isnan(ry)) {
+            if (mu_ns::isnan(ry)) {
                 return resT{ry, ry};
             }
             /* ordinary cases */
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/asin.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/asin.hpp
index a150b8e2e1..eb75daf946 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/asin.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/asin.hpp
@@ -31,6 +31,7 @@
 
 #include "kernels/elementwise_functions/common.hpp"
 
+#include "utils/math_utils.hpp"
 #include "utils/offset_utils.hpp"
 #include "utils/type_dispatch.hpp"
 #include "utils/type_utils.hpp"
@@ -47,6 +48,7 @@ namespace asin
 
 namespace py = pybind11;
 namespace td_ns = dpctl::tensor::type_dispatch;
+namespace mu_ns = dpctl::tensor::math_utils;
 
 using dpctl::tensor::type_utils::is_complex;
 
@@ -80,9 +82,9 @@ template <typename argT, typename resT> struct AsinFunctor
             const realT x = std::imag(in);
             const realT y = std::real(in);
 
-            if (std::isnan(x)) {
+            if (mu_ns::isnan(x)) {
                 /* asinh(NaN + I*+-Inf) = opt(+-)Inf + I*NaN */
-                if (std::isinf(y)) {
+                if (mu_ns::isinf(y)) {
                     const realT asinh_re = y;
                     const realT asinh_im = q_nan;
                     return resT{asinh_im, asinh_re};
@@ -96,9 +98,9 @@ template <typename argT, typename resT> struct AsinFunctor
                 /* All other cases involving NaN return NaN + I*NaN. */
                 return resT{q_nan, q_nan};
             }
-            else if (std::isnan(y)) {
+            else if (mu_ns::isnan(y)) {
                 /* asinh(+-Inf + I*NaN) = +-Inf + I*NaN */
-                if (std::isinf(x)) {
+                if (mu_ns::isinf(x)) {
                     const realT asinh_re = x;
                     const realT asinh_im = q_nan;
                     return resT{asinh_im, asinh_re};
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/asinh.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/asinh.hpp
index a42ea02598..c198d68c61 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/asinh.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/asinh.hpp
@@ -31,6 +31,7 @@
 
 #include "kernels/elementwise_functions/common.hpp"
 
+#include "utils/math_utils.hpp"
 #include "utils/offset_utils.hpp"
 #include "utils/type_dispatch.hpp"
 #include "utils/type_utils.hpp"
@@ -47,6 +48,7 @@ namespace asinh
 
 namespace py = pybind11;
 namespace td_ns = dpctl::tensor::type_dispatch;
+namespace mu_ns = dpctl::tensor::math_utils;
 
 using dpctl::tensor::type_utils::is_complex;
 
@@ -73,9 +75,9 @@ template <typename argT, typename resT> struct AsinhFunctor
             const realT x = std::real(in);
             const realT y = std::imag(in);
 
-            if (std::isnan(x)) {
+            if (mu_ns::isnan(x)) {
                 /* asinh(NaN + I*+-Inf) = opt(+-)Inf + I*NaN */
-                if (std::isinf(y)) {
+                if (mu_ns::isinf(y)) {
                     return resT{y, q_nan};
                 }
                 /* asinh(NaN + I*0) = NaN + I*0 */
@@ -86,9 +88,9 @@ template <typename argT, typename resT> struct AsinhFunctor
                 return resT{q_nan, q_nan};
             }
 
-            if (std::isnan(y)) {
+            if (mu_ns::isnan(y)) {
                 /* asinh(+-Inf + I*NaN) = +-Inf + I*NaN */
-                if (std::isinf(x)) {
+                if (mu_ns::isinf(x)) {
                     return resT{x, q_nan};
                 }
                 /* All other cases involving NaN return NaN + I*NaN. */
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/atan.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/atan.hpp
index cc988bd56d..f838aab890 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/atan.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/atan.hpp
@@ -32,6 +32,7 @@
 
 #include "kernels/elementwise_functions/common.hpp"
 
+#include "utils/math_utils.hpp"
 #include "utils/offset_utils.hpp"
 #include "utils/type_dispatch.hpp"
 #include "utils/type_utils.hpp"
@@ -48,6 +49,7 @@ namespace atan
 
 namespace py = pybind11;
 namespace td_ns = dpctl::tensor::type_dispatch;
+namespace mu_ns = dpctl::tensor::math_utils;
 
 using dpctl::tensor::type_utils::is_complex;
 
@@ -79,9 +81,9 @@ template <typename argT, typename resT> struct AtanFunctor
              */
             const realT x = std::imag(in);
             const realT y = std::real(in);
-            if (std::isnan(x)) {
+            if (mu_ns::isnan(x)) {
                 /* atanh(NaN + I*+-Inf) = sign(NaN)*0 + I*+-Pi/2 */
-                if (std::isinf(y)) {
+                if (mu_ns::isinf(y)) {
                     const realT pi_half = std::atan(realT(1)) * 2;
 
                     const realT atanh_re = std::copysign(realT(0), x);
@@ -93,9 +95,9 @@ template <typename argT, typename resT> struct AtanFunctor
                  */
                 return resT{q_nan, q_nan};
             }
-            else if (std::isnan(y)) {
+            else if (mu_ns::isnan(y)) {
                 /* atanh(+-Inf + I*NaN) = +-0 + I*NaN */
-                if (std::isinf(x)) {
+                if (mu_ns::isinf(x)) {
                     const realT atanh_re = std::copysign(realT(0), x);
                     const realT atanh_im = q_nan;
                     return resT{atanh_im, atanh_re};
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/atan2.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/atan2.hpp
index 613e2bff85..ebd710ec7e 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/atan2.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/atan2.hpp
@@ -29,6 +29,7 @@
 #include <cstdint>
 #include <type_traits>
 
+#include "utils/math_utils.hpp"
 #include "utils/offset_utils.hpp"
 #include "utils/type_dispatch.hpp"
 #include "utils/type_utils.hpp"
@@ -48,6 +49,7 @@ namespace atan2
 namespace py = pybind11;
 namespace td_ns = dpctl::tensor::type_dispatch;
 namespace tu_ns = dpctl::tensor::type_utils;
+namespace mu_ns = dpctl::tensor::math_utils;
 
 template <typename argT1, typename argT2, typename resT> struct Atan2Functor
 {
@@ -57,8 +59,8 @@ template <typename argT1, typename argT2, typename resT> struct Atan2Functor
 
     resT operator()(const argT1 &in1, const argT2 &in2)
     {
-        if (std::isinf(in2) && !std::signbit(in2)) {
-            if (std::isfinite(in1)) {
+        if (mu_ns::isinf(in2) && !std::signbit(in2)) {
+            if (mu_ns::isfinite(in1)) {
                 return std::copysign(resT(0), in1);
             }
         }
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/atanh.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/atanh.hpp
index 4842f785a2..30dea5f482 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/atanh.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/atanh.hpp
@@ -32,6 +32,7 @@
 
 #include "kernels/elementwise_functions/common.hpp"
 
+#include "utils/math_utils.hpp"
 #include "utils/offset_utils.hpp"
 #include "utils/type_dispatch.hpp"
 #include "utils/type_utils.hpp"
@@ -48,6 +49,7 @@ namespace atanh
 
 namespace py = pybind11;
 namespace td_ns = dpctl::tensor::type_dispatch;
+namespace mu_ns = dpctl::tensor::math_utils;
 
 using dpctl::tensor::type_utils::is_complex;
 
@@ -73,9 +75,9 @@ template <typename argT, typename resT> struct AtanhFunctor
             const realT x = std::real(in);
             const realT y = std::imag(in);
 
-            if (std::isnan(x)) {
+            if (mu_ns::isnan(x)) {
                 /* atanh(NaN + I*+-Inf) = sign(NaN)0 + I*+-PI/2 */
-                if (std::isinf(y)) {
+                if (mu_ns::isinf(y)) {
                     const realT pi_half = std::atan(realT(1)) * 2;
 
                     const realT res_re = std::copysign(realT(0), x);
@@ -87,9 +89,9 @@ template <typename argT, typename resT> struct AtanhFunctor
                  */
                 return resT{q_nan, q_nan};
             }
-            else if (std::isnan(y)) {
+            else if (mu_ns::isnan(y)) {
                 /* atanh(+-Inf + I*NaN) = +-0 + I*NaN */
-                if (std::isinf(x)) {
+                if (mu_ns::isinf(x)) {
                     const realT res_re = std::copysign(realT(0), x);
                     return resT{res_re, q_nan};
                 }
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/cos.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/cos.hpp
index 313fb565e8..98986c813f 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/cos.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/cos.hpp
@@ -31,6 +31,7 @@
 
 #include "kernels/elementwise_functions/common.hpp"
 
+#include "utils/math_utils.hpp"
 #include "utils/offset_utils.hpp"
 #include "utils/type_dispatch.hpp"
 #include "utils/type_utils.hpp"
@@ -47,6 +48,7 @@ namespace cos
 
 namespace py = pybind11;
 namespace td_ns = dpctl::tensor::type_dispatch;
+namespace mu_ns = dpctl::tensor::math_utils;
 
 using dpctl::tensor::type_utils::is_complex;
 
@@ -73,8 +75,8 @@ template <typename argT, typename resT> struct CosFunctor
             realT const &in_re = std::real(in);
             realT const &in_im = std::imag(in);
 
-            const bool in_re_finite = std::isfinite(in_re);
-            const bool in_im_finite = std::isfinite(in_im);
+            const bool in_re_finite = mu_ns::isfinite(in_re);
+            const bool in_im_finite = mu_ns::isfinite(in_im);
 
             /*
              * Handle the nearly-non-exceptional cases where
@@ -137,7 +139,7 @@ template <typename argT, typename resT> struct CosFunctor
              *
              * cosh(+-Inf + I y)   = +Inf cos(y) +- I Inf sin(y)
              */
-            if (std::isinf(x)) {
+            if (mu_ns::isinf(x)) {
                 if (!yfinite) {
                     return resT{x * x, std::copysign(q_nan, x)};
                 }
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/cosh.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/cosh.hpp
index 317cd101af..3ed5610fe9 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/cosh.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/cosh.hpp
@@ -31,6 +31,7 @@
 
 #include "kernels/elementwise_functions/common.hpp"
 
+#include "utils/math_utils.hpp"
 #include "utils/offset_utils.hpp"
 #include "utils/type_dispatch.hpp"
 #include "utils/type_utils.hpp"
@@ -47,6 +48,7 @@ namespace cosh
 
 namespace py = pybind11;
 namespace td_ns = dpctl::tensor::type_dispatch;
+namespace mu_ns = dpctl::tensor::math_utils;
 
 using dpctl::tensor::type_utils::is_complex;
 
@@ -73,8 +75,8 @@ template <typename argT, typename resT> struct CoshFunctor
             const realT x = std::real(in);
             const realT y = std::imag(in);
 
-            const bool xfinite = std::isfinite(x);
-            const bool yfinite = std::isfinite(y);
+            const bool xfinite = mu_ns::isfinite(x);
+            const bool yfinite = mu_ns::isfinite(y);
 
             /*
              * Handle the nearly-non-exceptional cases where
@@ -126,7 +128,7 @@ template <typename argT, typename resT> struct CoshFunctor
              *
              * cosh(+-Inf + I y)   = +Inf cos(y) +- I Inf sin(y)
              */
-            if (std::isinf(x)) {
+            if (mu_ns::isinf(x)) {
                 if (!yfinite) {
                     return resT{x * x, x * q_nan};
                 }
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/exp.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/exp.hpp
index 2f7657f5b4..2fc1843960 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/exp.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/exp.hpp
@@ -31,6 +31,7 @@
 
 #include "kernels/elementwise_functions/common.hpp"
 
+#include "utils/math_utils.hpp"
 #include "utils/offset_utils.hpp"
 #include "utils/type_dispatch.hpp"
 #include "utils/type_utils.hpp"
@@ -47,6 +48,7 @@ namespace exp
 
 namespace py = pybind11;
 namespace td_ns = dpctl::tensor::type_dispatch;
+namespace mu_ns = dpctl::tensor::math_utils;
 
 using dpctl::tensor::type_utils::is_complex;
 
@@ -71,15 +73,15 @@ template <typename argT, typename resT> struct ExpFunctor
 
             const realT x = std::real(in);
             const realT y = std::imag(in);
-            if (std::isfinite(x)) {
-                if (std::isfinite(y)) {
+            if (mu_ns::isfinite(x)) {
+                if (mu_ns::isfinite(y)) {
                     return std::exp(in);
                 }
                 else {
                     return resT{q_nan, q_nan};
                 }
             }
-            else if (std::isnan(x)) {
+            else if (mu_ns::isnan(x)) {
                 /* x is nan */
                 if (y == realT(0)) {
                     return resT{in};
@@ -93,7 +95,7 @@ template <typename argT, typename resT> struct ExpFunctor
                     if (y == realT(0)) {
                         return resT{x, y};
                     }
-                    else if (std::isfinite(y)) {
+                    else if (mu_ns::isfinite(y)) {
                         return resT{x * std::cos(y), x * std::sin(y)};
                     }
                     else {
@@ -102,7 +104,7 @@ template <typename argT, typename resT> struct ExpFunctor
                     }
                 }
                 else { /* x is -inf */
-                    if (std::isfinite(y)) {
+                    if (mu_ns::isfinite(y)) {
                         realT exp_x = std::exp(x);
                         return resT{exp_x * std::cos(y), exp_x * std::sin(y)};
                     }
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/expm1.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/expm1.hpp
index e1c23113c6..76100b249d 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/expm1.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/expm1.hpp
@@ -33,6 +33,7 @@
 
 #include "kernels/elementwise_functions/common.hpp"
 
+#include "utils/math_utils.hpp"
 #include "utils/offset_utils.hpp"
 #include "utils/type_dispatch.hpp"
 #include "utils/type_utils.hpp"
@@ -49,6 +50,7 @@ namespace expm1
 
 namespace py = pybind11;
 namespace td_ns = dpctl::tensor::type_dispatch;
+namespace mu_ns = dpctl::tensor::math_utils;
 
 using dpctl::tensor::type_utils::is_complex;
 
@@ -75,10 +77,10 @@ template <typename argT, typename resT> struct Expm1Functor
             const realT y = std::imag(in);
 
             // special cases
-            if (std::isinf(x)) {
+            if (mu_ns::isinf(x)) {
                 if (x > realT(0)) {
                     // positive infinity cases
-                    if (!std::isfinite(y)) {
+                    if (!mu_ns::isfinite(y)) {
                         return resT{x, std::numeric_limits<realT>::quiet_NaN()};
                     }
                     else if (y == realT(0)) {
@@ -91,7 +93,7 @@ template <typename argT, typename resT> struct Expm1Functor
                 }
                 else {
                     // negative infinity cases
-                    if (!std::isfinite(y)) {
+                    if (!mu_ns::isfinite(y)) {
                         // copy sign of y to guarantee
                         // conj(expm1(x)) == expm1(conj(x))
                         return resT{realT(-1), std::copysign(realT(0), y)};
@@ -103,7 +105,7 @@ template <typename argT, typename resT> struct Expm1Functor
                 }
             }
 
-            if (std::isnan(x)) {
+            if (mu_ns::isnan(x)) {
                 if (y == realT(0)) {
                     return in;
                 }
@@ -119,6 +121,7 @@ template <typename argT, typename resT> struct Expm1Functor
                 sycl::access::address_space::private_space,
                 sycl::access::decorated::yes>(&cosY_val);
             const realT sinY_val = sycl::sincos(y, cosY_val_multi_ptr);
+
             const realT sinhalfY_val = std::sin(y / 2);
 
             const realT res_re =
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/isfinite.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/isfinite.hpp
index cc91bdbb97..aedf78829f 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/isfinite.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/isfinite.hpp
@@ -30,6 +30,7 @@
 #include <cstdint>
 #include <type_traits>
 
+#include "utils/math_utils.hpp"
 #include "utils/offset_utils.hpp"
 #include "utils/type_dispatch.hpp"
 #include "utils/type_utils.hpp"
@@ -46,6 +47,7 @@ namespace isfinite
 
 namespace py = pybind11;
 namespace td_ns = dpctl::tensor::type_dispatch;
+namespace mu_ns = dpctl::tensor::math_utils;
 
 using dpctl::tensor::type_utils::is_complex;
 using dpctl::tensor::type_utils::vec_cast;
@@ -68,8 +70,8 @@ template <typename argT, typename resT> struct IsFiniteFunctor
     resT operator()(const argT &in) const
     {
         if constexpr (is_complex<argT>::value) {
-            const bool real_isfinite = std::isfinite(std::real(in));
-            const bool imag_isfinite = std::isfinite(std::imag(in));
+            const bool real_isfinite = mu_ns::isfinite(std::real(in));
+            const bool imag_isfinite = mu_ns::isfinite(std::imag(in));
             return (real_isfinite && imag_isfinite);
         }
         else if constexpr (std::is_same<argT, bool>::value ||
@@ -77,24 +79,10 @@ template <typename argT, typename resT> struct IsFiniteFunctor
         {
             return constant_value;
         }
-        else if constexpr (std::is_same_v<argT, sycl::half>) {
-            return sycl::isfinite(in);
-        }
         else {
-            return std::isfinite(in);
+            return mu_ns::isfinite(in);
         }
     }
-
-    template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT, vec_sz> &in)
-    {
-        auto const &res_vec = sycl::isfinite(in);
-
-        using deducedT = typename std::remove_cv_t<
-            std::remove_reference_t<decltype(res_vec)>>::element_type;
-
-        return vec_cast<bool, deducedT, vec_sz>(res_vec);
-    }
 };
 
 template <typename argT,
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/isinf.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/isinf.hpp
index 8cb93ea456..320f2bce16 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/isinf.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/isinf.hpp
@@ -30,6 +30,7 @@
 #include <cstdint>
 #include <type_traits>
 
+#include "utils/math_utils.hpp"
 #include "utils/offset_utils.hpp"
 #include "utils/type_dispatch.hpp"
 #include "utils/type_utils.hpp"
@@ -46,6 +47,7 @@ namespace isinf
 
 namespace py = pybind11;
 namespace td_ns = dpctl::tensor::type_dispatch;
+namespace mu_ns = dpctl::tensor::math_utils;
 
 using dpctl::tensor::type_utils::is_complex;
 using dpctl::tensor::type_utils::vec_cast;
@@ -66,8 +68,8 @@ template <typename argT, typename resT> struct IsInfFunctor
     resT operator()(const argT &in) const
     {
         if constexpr (is_complex<argT>::value) {
-            const bool real_isinf = std::isinf(std::real(in));
-            const bool imag_isinf = std::isinf(std::imag(in));
+            const bool real_isinf = mu_ns::isinf(std::real(in));
+            const bool imag_isinf = mu_ns::isinf(std::imag(in));
             return (real_isinf || imag_isinf);
         }
         else if constexpr (std::is_same<argT, bool>::value ||
@@ -79,7 +81,7 @@ template <typename argT, typename resT> struct IsInfFunctor
             return sycl::isinf(in);
         }
         else {
-            return std::isinf(in);
+            return mu_ns::isinf(in);
         }
     }
 
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/isnan.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/isnan.hpp
index 5d628437ed..ee1d31b52d 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/isnan.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/isnan.hpp
@@ -29,6 +29,7 @@
 #include <cstdint>
 #include <type_traits>
 
+#include "utils/math_utils.hpp"
 #include "utils/offset_utils.hpp"
 #include "utils/type_dispatch.hpp"
 #include "utils/type_utils.hpp"
@@ -45,6 +46,7 @@ namespace isnan
 
 namespace py = pybind11;
 namespace td_ns = dpctl::tensor::type_dispatch;
+namespace mu_ns = dpctl::tensor::math_utils;
 
 using dpctl::tensor::type_utils::is_complex;
 using dpctl::tensor::type_utils::vec_cast;
@@ -77,7 +79,7 @@ template <typename argT, typename resT> struct IsNanFunctor
             return constant_value;
         }
         else {
-            return sycl::isnan(in);
+            return mu_ns::isnan(in);
         }
     }
 
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/logaddexp.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/logaddexp.hpp
index b718a5f991..85dc44c0aa 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/logaddexp.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/logaddexp.hpp
@@ -31,6 +31,7 @@
 #include <limits>
 #include <type_traits>
 
+#include "utils/math_utils.hpp"
 #include "utils/offset_utils.hpp"
 #include "utils/type_dispatch.hpp"
 #include "utils/type_utils.hpp"
@@ -48,6 +49,7 @@ namespace logaddexp
 {
 
 namespace py = pybind11;
+namespace mu_ns = dpctl::tensor::math_utils;
 namespace td_ns = dpctl::tensor::type_dispatch;
 namespace tu_ns = dpctl::tensor::type_utils;
 
@@ -73,7 +75,7 @@ template <typename argT1, typename argT2, typename resT> struct LogAddExpFunctor
 
 #pragma unroll
         for (int i = 0; i < vec_sz; ++i) {
-            if (std::isfinite(diff[i])) {
+            if (mu_ns::isfinite(diff[i])) {
                 res[i] = std::max<resT>(in1[i], in2[i]) +
                          impl_finite<resT>(-std::abs(diff[i]));
             }
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/maximum.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/maximum.hpp
index 60afc2fc15..fc820c90b9 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/maximum.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/maximum.hpp
@@ -49,6 +49,7 @@ namespace maximum
 namespace py = pybind11;
 namespace td_ns = dpctl::tensor::type_dispatch;
 namespace tu_ns = dpctl::tensor::type_utils;
+namespace mu_ns = dpctl::tensor::math_utils;
 
 template <typename argT1, typename argT2, typename resT> struct MaximumFunctor
 {
@@ -66,12 +67,11 @@ template <typename argT1, typename argT2, typename resT> struct MaximumFunctor
                       tu_ns::is_complex<argT2>::value)
         {
             static_assert(std::is_same_v<argT1, argT2>);
-            using dpctl::tensor::math_utils::max_complex;
-            return max_complex<argT1>(in1, in2);
+            return mu_ns::max_complex<argT1>(in1, in2);
         }
         else if constexpr (std::is_floating_point_v<argT1> ||
                            std::is_same_v<argT1, sycl::half>)
-            return (std::isnan(in1) || in1 > in2) ? in1 : in2;
+            return (mu_ns::isnan(in1) || in1 > in2) ? in1 : in2;
         else
             return (in1 > in2) ? in1 : in2;
     }
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/minimum.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/minimum.hpp
index d19c829a6c..99ad45ce43 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/minimum.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/minimum.hpp
@@ -49,6 +49,7 @@ namespace minimum
 namespace py = pybind11;
 namespace td_ns = dpctl::tensor::type_dispatch;
 namespace tu_ns = dpctl::tensor::type_utils;
+namespace mu_ns = dpctl::tensor::math_utils;
 
 template <typename argT1, typename argT2, typename resT> struct MinimumFunctor
 {
@@ -66,12 +67,11 @@ template <typename argT1, typename argT2, typename resT> struct MinimumFunctor
                       tu_ns::is_complex<argT2>::value)
         {
             static_assert(std::is_same_v<argT1, argT2>);
-            using dpctl::tensor::math_utils::min_complex;
-            return min_complex<argT1>(in1, in2);
+            return mu_ns::min_complex<argT1>(in1, in2);
         }
         else if constexpr (std::is_floating_point_v<argT1> ||
                            std::is_same_v<argT1, sycl::half>)
-            return (std::isnan(in1) || in1 < in2) ? in1 : in2;
+            return (mu_ns::isnan(in1) || in1 < in2) ? in1 : in2;
         else
             return (in1 < in2) ? in1 : in2;
     }
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/proj.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/proj.hpp
index 6fbc3ae012..6b6a5864a5 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/proj.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/proj.hpp
@@ -34,6 +34,7 @@
 
 #include "kernels/elementwise_functions/common.hpp"
 
+#include "utils/math_utils.hpp"
 #include "utils/offset_utils.hpp"
 #include "utils/type_dispatch.hpp"
 #include "utils/type_utils.hpp"
@@ -50,6 +51,7 @@ namespace proj
 
 namespace py = pybind11;
 namespace td_ns = dpctl::tensor::type_dispatch;
+namespace mu_ns = dpctl::tensor::math_utils;
 
 using dpctl::tensor::type_utils::is_complex;
 
@@ -71,7 +73,7 @@ template <typename argT, typename resT> struct ProjFunctor
         const realT x = std::real(in);
         const realT y = std::imag(in);
 
-        if (std::isinf(x) || std::isinf(y)) {
+        if (mu_ns::isinf(x) || mu_ns::isinf(y)) {
             const realT res_im = std::copysign(realT(0), y);
             return resT{std::numeric_limits<realT>::infinity(), res_im};
         }
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/sign.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/sign.hpp
index 8f0488746e..c5e970e9f0 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/sign.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/sign.hpp
@@ -32,6 +32,7 @@
 
 #include "kernels/elementwise_functions/common.hpp"
 
+#include "utils/math_utils.hpp"
 #include "utils/offset_utils.hpp"
 #include "utils/type_dispatch.hpp"
 #include "utils/type_utils.hpp"
@@ -48,6 +49,7 @@ namespace sign
 
 namespace py = pybind11;
 namespace td_ns = dpctl::tensor::type_dispatch;
+namespace mu_ns = dpctl::tensor::math_utils;
 
 using dpctl::tensor::type_utils::is_complex;
 using dpctl::tensor::type_utils::vec_cast;
@@ -81,7 +83,7 @@ template <typename argT, typename resT> struct SignFunctor
                 }
             }
             else {
-                if (std::isnan(x)) {
+                if (mu_ns::isnan(x)) {
                     return std::numeric_limits<resT>::quiet_NaN();
                 }
                 else {
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/sin.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/sin.hpp
index 74cab65ec4..40954e722c 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/sin.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/sin.hpp
@@ -31,6 +31,7 @@
 
 #include "kernels/elementwise_functions/common.hpp"
 
+#include "utils/math_utils.hpp"
 #include "utils/offset_utils.hpp"
 #include "utils/type_dispatch.hpp"
 #include "utils/type_utils.hpp"
@@ -46,6 +47,7 @@ namespace sin
 {
 
 namespace py = pybind11;
+namespace mu_ns = dpctl::tensor::math_utils;
 namespace td_ns = dpctl::tensor::type_dispatch;
 
 using dpctl::tensor::type_utils::is_complex;
@@ -72,8 +74,8 @@ template <typename argT, typename resT> struct SinFunctor
             realT const &in_re = std::real(in);
             realT const &in_im = std::imag(in);
 
-            const bool in_re_finite = std::isfinite(in_re);
-            const bool in_im_finite = std::isfinite(in_im);
+            const bool in_re_finite = mu_ns::isfinite(in_re);
+            const bool in_im_finite = mu_ns::isfinite(in_im);
             /*
              * Handle the nearly-non-exceptional cases where
              * real and imaginary parts of input are finite.
@@ -112,7 +114,7 @@ template <typename argT, typename resT> struct SinFunctor
              * sinh(NaN +- I 0)   = d(NaN) + I +-0.
              */
             if (y == realT(0) && !xfinite) {
-                if (std::isnan(x)) {
+                if (mu_ns::isnan(x)) {
                     const realT sinh_re = x;
                     const realT sinh_im = y;
                     return resT{sinh_im, -sinh_re};
@@ -145,7 +147,7 @@ template <typename argT, typename resT> struct SinFunctor
              *
              * sinh(+-Inf + I y)   = +-Inf cos(y) + I Inf sin(y)
              */
-            if (std::isinf(x)) {
+            if (mu_ns::isinf(x)) {
                 if (!yfinite) {
                     const realT sinh_re = -x * x;
                     const realT sinh_im = x * (y - y);
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/sinh.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/sinh.hpp
index f535f3e946..089b61e7d2 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/sinh.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/sinh.hpp
@@ -31,6 +31,7 @@
 
 #include "kernels/elementwise_functions/common.hpp"
 
+#include "utils/math_utils.hpp"
 #include "utils/offset_utils.hpp"
 #include "utils/type_dispatch.hpp"
 #include "utils/type_utils.hpp"
@@ -46,6 +47,7 @@ namespace sinh
 {
 
 namespace py = pybind11;
+namespace mu_ns = dpctl::tensor::math_utils;
 namespace td_ns = dpctl::tensor::type_dispatch;
 
 using dpctl::tensor::type_utils::is_complex;
@@ -71,8 +73,8 @@ template <typename argT, typename resT> struct SinhFunctor
             const realT x = std::real(in);
             const realT y = std::imag(in);
 
-            const bool xfinite = std::isfinite(x);
-            const bool yfinite = std::isfinite(y);
+            const bool xfinite = mu_ns::isfinite(x);
+            const bool yfinite = mu_ns::isfinite(y);
 
             /*
              * Handle the nearly-non-exceptional cases where
@@ -101,7 +103,7 @@ template <typename argT, typename resT> struct SinhFunctor
              * sinh(NaN +- I 0)   = d(NaN) + I +-0.
              */
             if (y == realT(0) && !xfinite) {
-                if (std::isnan(x)) {
+                if (mu_ns::isnan(x)) {
                     return resT{x, y};
                 }
                 const realT res_im = std::copysign(realT(0), y);
@@ -127,7 +129,7 @@ template <typename argT, typename resT> struct SinhFunctor
              *
              * sinh(+-Inf + I y)   = +-Inf cos(y) + I Inf sin(y)
              */
-            if (!xfinite && !std::isnan(x)) {
+            if (!xfinite && !mu_ns::isnan(x)) {
                 if (!yfinite) {
                     return resT{x * x, x * (y - y)};
                 }
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/sqrt.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/sqrt.hpp
index 9b9fa95061..cbf7db4c75 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/sqrt.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/sqrt.hpp
@@ -34,6 +34,7 @@
 
 #include "kernels/elementwise_functions/common.hpp"
 
+#include "utils/math_utils.hpp"
 #include "utils/offset_utils.hpp"
 #include "utils/type_dispatch.hpp"
 #include "utils/type_utils.hpp"
@@ -49,6 +50,7 @@ namespace sqrt
 {
 
 namespace py = pybind11;
+namespace mu_ns = dpctl::tensor::math_utils;
 namespace td_ns = dpctl::tensor::type_dispatch;
 
 using dpctl::tensor::type_utils::is_complex;
@@ -108,23 +110,25 @@ template <typename argT, typename resT> struct SqrtFunctor
         realT x = std::real(z);
         realT y = std::imag(z);
 
-        if (std::isinf(y)) {
+        if (mu_ns::isinf(y)) {
             return {p_inf, y};
         }
-        else if (std::isnan(x)) {
+        else if (mu_ns::isnan(x)) {
             return {x, q_nan};
         }
-        else if (std::isinf(x)) { // x is an infinity
+        else if (mu_ns::isinf(x)) { // x is an infinity
             // y is either finite, or nan
             if (std::signbit(x)) { // x == -inf
-                return {(std::isfinite(y) ? zero : y), std::copysign(p_inf, y)};
+                return {(mu_ns::isfinite(y) ? zero : y),
+                        std::copysign(p_inf, y)};
             }
             else {
-                return {p_inf, (std::isfinite(y) ? std::copysign(zero, y) : y)};
+                return {p_inf,
+                        (mu_ns::isfinite(y) ? std::copysign(zero, y) : y)};
             }
         }
         else { // x is finite
-            if (std::isfinite(y)) {
+            if (mu_ns::isfinite(y)) {
 #ifdef USE_STD_SQRT_FOR_COMPLEX_TYPES
                 return std::sqrt(z);
 #else
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/tan.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/tan.hpp
index 3f11b0e721..3c4a53f9e8 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/tan.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/tan.hpp
@@ -32,6 +32,7 @@
 
 #include "kernels/elementwise_functions/common.hpp"
 
+#include "utils/math_utils.hpp"
 #include "utils/offset_utils.hpp"
 #include "utils/type_dispatch.hpp"
 #include "utils/type_utils.hpp"
@@ -47,6 +48,7 @@ namespace tan
 {
 
 namespace py = pybind11;
+namespace mu_ns = dpctl::tensor::math_utils;
 namespace td_ns = dpctl::tensor::type_dispatch;
 
 using dpctl::tensor::type_utils::is_complex;
@@ -94,15 +96,15 @@ template <typename argT, typename resT> struct TanFunctor
              * case is only needed to avoid a spurious invalid exception when
              * y is infinite.
              */
-            if (!std::isfinite(x)) {
-                if (std::isnan(x)) {
+            if (!mu_ns::isfinite(x)) {
+                if (mu_ns::isnan(x)) {
                     const realT tanh_re = x;
                     const realT tanh_im = (y == realT(0) ? y : x * y);
                     return resT{tanh_im, -tanh_re};
                 }
                 const realT tanh_re = std::copysign(realT(1), x);
                 const realT tanh_im = std::copysign(
-                    realT(0), std::isinf(y) ? y : std::sin(y) * std::cos(y));
+                    realT(0), mu_ns::isinf(y) ? y : std::sin(y) * std::cos(y));
                 return resT{tanh_im, -tanh_re};
             }
             /*
@@ -111,7 +113,7 @@ template <typename argT, typename resT> struct TanFunctor
              * tanh(0 + i NAN) = 0 + i NaN
              * tanh(0 +- i Inf) = 0 + i NaN
              */
-            if (!std::isfinite(y)) {
+            if (!mu_ns::isfinite(y)) {
                 if (x == realT(0)) {
                     return resT{q_nan, x};
                 }
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/tanh.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/tanh.hpp
index 4ee18eccb4..19bcb4af37 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/tanh.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/tanh.hpp
@@ -33,6 +33,7 @@
 
 #include "kernels/elementwise_functions/common.hpp"
 
+#include "utils/math_utils.hpp"
 #include "utils/offset_utils.hpp"
 #include "utils/type_dispatch.hpp"
 #include "utils/type_utils.hpp"
@@ -48,6 +49,7 @@ namespace tanh
 {
 
 namespace py = pybind11;
+namespace mu_ns = dpctl::tensor::math_utils;
 namespace td_ns = dpctl::tensor::type_dispatch;
 
 using dpctl::tensor::type_utils::is_complex;
@@ -90,13 +92,13 @@ template <typename argT, typename resT> struct TanhFunctor
              * case is only needed to avoid a spurious invalid exception when
              * y is infinite.
              */
-            if (!std::isfinite(x)) {
-                if (std::isnan(x)) {
+            if (!mu_ns::isfinite(x)) {
+                if (mu_ns::isnan(x)) {
                     return resT{q_nan, (y == realT(0) ? y : q_nan)};
                 }
                 const realT res_re = std::copysign(realT(1), x);
                 const realT res_im = std::copysign(
-                    realT(0), std::isinf(y) ? y : std::sin(y) * std::cos(y));
+                    realT(0), mu_ns::isinf(y) ? y : std::sin(y) * std::cos(y));
                 return resT{res_re, res_im};
             }
             /*
@@ -105,7 +107,7 @@ template <typename argT, typename resT> struct TanhFunctor
              * tanh(0 + i NAN) = 0 + i NaN
              * tanh(0 +- i Inf) = 0 + i NaN
              */
-            if (!std::isfinite(y)) {
+            if (!mu_ns::isfinite(y)) {
                 if (x == realT(0)) {
                     return resT{x, q_nan};
                 }
diff --git a/dpctl/tensor/libtensor/include/utils/math_utils.hpp b/dpctl/tensor/libtensor/include/utils/math_utils.hpp
index d724e03e35..8c218f8f32 100644
--- a/dpctl/tensor/libtensor/include/utils/math_utils.hpp
+++ b/dpctl/tensor/libtensor/include/utils/math_utils.hpp
@@ -23,8 +23,13 @@
 //===----------------------------------------------------------------------===//
 
 #pragma once
+#include <CL/sycl.hpp>
 #include <cmath>
 #include <complex>
+#include <cstddef>
+#include <cstdint>
+
+#include <utils/type_utils.hpp>
 
 namespace dpctl
 {
@@ -33,86 +38,199 @@ namespace tensor
 namespace math_utils
 {
 
+namespace detail
+{
+namespace dt_tu = dpctl::tensor::type_utils;
+
+template <typename T> constexpr bool is_complex_v = dt_tu::is_complex<T>::value;
+
+template <typename T,
+          typename bitseqT,
+          typename scaleT,
+          bitseqT significand_bits,
+          bitseqT exponent_mask>
+bool is_finite(const T &v)
+{
+    const scaleT scale = static_cast<scaleT>(
+        (sycl::bit_cast<bitseqT>(v) >> significand_bits) & exponent_mask);
+
+    return static_cast<bool>(scale ^ static_cast<scaleT>(exponent_mask));
+}
+
+template <typename T,
+          typename bitseqT,
+          typename scaleT,
+          bitseqT significand_bits,
+          bitseqT exponent_mask>
+bool is_inf(const T &v)
+{
+    constexpr bitseqT significand_mask = ((bitseqT(1) << significand_bits) - 1);
+
+    const bitseqT bits = sycl::bit_cast<bitseqT>(v);
+    const scaleT scale =
+        static_cast<scaleT>((bits >> significand_bits) & exponent_mask);
+
+    const bitseqT significand = bits & significand_mask;
+
+    return (!static_cast<bool>(scale ^ static_cast<scaleT>(exponent_mask)) &&
+            !static_cast<bool>(significand));
+}
+
+template <typename T,
+          typename bitseqT,
+          typename scaleT,
+          bitseqT significand_bits,
+          bitseqT exponent_mask>
+bool is_nan(const T &v)
+{
+    constexpr bitseqT significand_mask = ((bitseqT(1) << significand_bits) - 1);
+
+    const bitseqT bits = sycl::bit_cast<bitseqT>(v);
+    const scaleT scale =
+        static_cast<scaleT>((bits >> significand_bits) & exponent_mask);
+
+    const bitseqT significand = bits & significand_mask;
+
+    return (!static_cast<bool>(scale ^ static_cast<scaleT>(exponent_mask)) &&
+            static_cast<bool>(significand));
+}
+
+} // namespace detail
+
+// Work-arounds for bug in bug in SYCLOS
+template <typename T> bool isfinite(const T &x)
+{
+    static_assert(std::is_same_v<T, float> || std::is_same_v<T, double> ||
+                  std::is_same_v<T, sycl::half>);
+    if constexpr (std::is_same_v<T, float>) {
+        return detail::is_finite<T, std::uint32_t, std::uint32_t, 23, 0xff>(x);
+    }
+    else if constexpr (std::is_same_v<T, double>) {
+        return detail::is_finite<T, std::uint64_t, std::uint32_t, 52, 0x7ff>(x);
+    }
+    else if constexpr (std::is_same_v<T, sycl::half>) {
+        return detail::is_finite<T, std::uint16_t, std::uint16_t, 10, 0x1f>(x);
+    }
+}
+
+template <typename T> bool isinf(const T &x)
+{
+    static_assert(std::is_same_v<T, float> || std::is_same_v<T, double> ||
+                  std::is_same_v<T, sycl::half>);
+    if constexpr (std::is_same_v<T, float>) {
+        return detail::is_inf<T, std::uint32_t, std::uint32_t, 23, 0xff>(x);
+    }
+    else if constexpr (std::is_same_v<T, double>) {
+        return detail::is_inf<T, std::uint64_t, std::uint32_t, 52, 0x7ff>(x);
+    }
+    else if constexpr (std::is_same_v<T, sycl::half>) {
+        return detail::is_inf<T, std::uint16_t, std::uint16_t, 10, 0x1f>(x);
+    }
+}
+
+template <typename T> bool isnan(const T &x)
+{
+    static_assert(std::is_same_v<T, float> || std::is_same_v<T, double> ||
+                  std::is_same_v<T, sycl::half>);
+    if constexpr (std::is_same_v<T, float>) {
+        return detail::is_nan<T, std::uint32_t, std::uint32_t, 23, 0xff>(x);
+    }
+    else if constexpr (std::is_same_v<T, double>) {
+        return detail::is_nan<T, std::uint64_t, std::uint32_t, 52, 0x7ff>(x);
+    }
+    else if constexpr (std::is_same_v<T, sycl::half>) {
+        return detail::is_nan<T, std::uint16_t, std::uint16_t, 10, 0x1f>(x);
+    }
+}
+
 template <typename T> bool less_complex(const T &x1, const T &x2)
 {
+    static_assert(detail::is_complex_v<T>);
+
     using realT = typename T::value_type;
     realT real1 = std::real(x1);
     realT real2 = std::real(x2);
     realT imag1 = std::imag(x1);
     realT imag2 = std::imag(x2);
 
-    return (real1 == real2)
-               ? (imag1 < imag2)
-               : (real1 < real2 && !std::isnan(imag1) && !std::isnan(imag2));
+    return (real1 == real2) ? (imag1 < imag2)
+                            : (real1 < real2 && !isnan(imag1) && !isnan(imag2));
 }
 
 template <typename T> bool greater_complex(const T &x1, const T &x2)
 {
+    static_assert(detail::is_complex_v<T>);
+
     using realT = typename T::value_type;
     realT real1 = std::real(x1);
     realT real2 = std::real(x2);
     realT imag1 = std::imag(x1);
     realT imag2 = std::imag(x2);
 
-    return (real1 == real2)
-               ? (imag1 > imag2)
-               : (real1 > real2 && !std::isnan(imag1) && !std::isnan(imag2));
+    return (real1 == real2) ? (imag1 > imag2)
+                            : (real1 > real2 && !isnan(imag1) && !isnan(imag2));
 }
 
 template <typename T> bool less_equal_complex(const T &x1, const T &x2)
 {
+    static_assert(detail::is_complex_v<T>);
+
     using realT = typename T::value_type;
     realT real1 = std::real(x1);
     realT real2 = std::real(x2);
     realT imag1 = std::imag(x1);
     realT imag2 = std::imag(x2);
 
-    return (real1 == real2)
-               ? (imag1 <= imag2)
-               : (real1 < real2 && !std::isnan(imag1) && !std::isnan(imag2));
+    return (real1 == real2) ? (imag1 <= imag2)
+                            : (real1 < real2 && !isnan(imag1) && !isnan(imag2));
 }
 
 template <typename T> bool greater_equal_complex(const T &x1, const T &x2)
 {
+    static_assert(detail::is_complex_v<T>);
+
     using realT = typename T::value_type;
     realT real1 = std::real(x1);
     realT real2 = std::real(x2);
     realT imag1 = std::imag(x1);
     realT imag2 = std::imag(x2);
 
-    return (real1 == real2)
-               ? (imag1 >= imag2)
-               : (real1 > real2 && !std::isnan(imag1) && !std::isnan(imag2));
+    return (real1 == real2) ? (imag1 >= imag2)
+                            : (real1 > real2 && !isnan(imag1) && !isnan(imag2));
 }
 
 template <typename T> T max_complex(const T &x1, const T &x2)
 {
+    static_assert(detail::is_complex_v<T>);
+
     using realT = typename T::value_type;
     realT real1 = std::real(x1);
     realT real2 = std::real(x2);
     realT imag1 = std::imag(x1);
     realT imag2 = std::imag(x2);
 
-    bool isnan_imag1 = std::isnan(imag1);
+    bool isnan_imag1 = isnan(imag1);
     bool gt = (real1 == real2)
                   ? (imag1 > imag2)
-                  : (real1 > real2 && !isnan_imag1 && !std::isnan(imag2));
-    return (std::isnan(real1) || isnan_imag1 || gt) ? x1 : x2;
+                  : (real1 > real2 && !isnan_imag1 && !isnan(imag2));
+    return (isnan(real1) || isnan_imag1 || gt) ? x1 : x2;
 }
 
 template <typename T> T min_complex(const T &x1, const T &x2)
 {
+    static_assert(detail::is_complex_v<T>);
+
     using realT = typename T::value_type;
     realT real1 = std::real(x1);
     realT real2 = std::real(x2);
     realT imag1 = std::imag(x1);
     realT imag2 = std::imag(x2);
 
-    bool isnan_imag1 = std::isnan(imag1);
+    bool isnan_imag1 = isnan(imag1);
     bool lt = (real1 == real2)
                   ? (imag1 < imag2)
-                  : (real1 < real2 && !isnan_imag1 && !std::isnan(imag2));
-    return (std::isnan(real1) || isnan_imag1 || lt) ? x1 : x2;
+                  : (real1 < real2 && !isnan_imag1 && !isnan(imag2));
+    return (isnan(real1) || isnan_imag1 || lt) ? x1 : x2;
 }
 
 } // namespace math_utils