diff --git a/CHANGELOG.md b/CHANGELOG.md
index 15f1840317a8..9e4a741cc3f5 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -15,6 +15,7 @@ The release drops support for Python 3.9, making Python 3.10 the minimum require
 * Added implementation of `dpnp.linalg.lu_solve` for batch inputs (SciPy-compatible) [#2619](https://github.com/IntelPython/dpnp/pull/2619)
 * Added `dpnp.exceptions` submodule to aggregate the generic exceptions used by dpnp [#2616](https://github.com/IntelPython/dpnp/pull/2616)
 * Added implementation of `dpnp.scipy.special.erfcx` [#2596](https://github.com/IntelPython/dpnp/pull/2596)
+* Added implementation of `dpnp.scipy.special.erfinv` and `dpnp.scipy.special.erfcinv` [#2624](https://github.com/IntelPython/dpnp/pull/2624)
 
 ### Changed
 
diff --git a/dpnp/backend/extensions/ufunc/elementwise_functions/erf_funcs.cpp b/dpnp/backend/extensions/ufunc/elementwise_functions/erf_funcs.cpp
index fbd2ac211980..5254e50d3faf 100644
--- a/dpnp/backend/extensions/ufunc/elementwise_functions/erf_funcs.cpp
+++ b/dpnp/backend/extensions/ufunc/elementwise_functions/erf_funcs.cpp
@@ -213,6 +213,8 @@ static void populate(py::module_ m,
 MACRO_DEFINE_IMPL(erf, Erf);
 MACRO_DEFINE_IMPL(erfc, Erfc);
 MACRO_DEFINE_IMPL(erfcx, Erfcx);
+MACRO_DEFINE_IMPL(erfinv, Erfinv);
+MACRO_DEFINE_IMPL(erfcinv, Erfcinv);
 } // namespace impl
 
 void init_erf_funcs(py::module_ m)
@@ -236,5 +238,13 @@ void init_erf_funcs(py::module_ m)
     impl::populate<impl::ErfcxContigFactory, impl::ErfcxStridedFactory>(
         m, "_erfcx", "", impl::erfcx_contig_dispatch_vector,
         impl::erfcx_strided_dispatch_vector);
+
+    impl::populate<impl::ErfinvContigFactory, impl::ErfinvStridedFactory>(
+        m, "_erfinv", "", impl::erfinv_contig_dispatch_vector,
+        impl::erfinv_strided_dispatch_vector);
+
+    impl::populate<impl::ErfcinvContigFactory, impl::ErfcinvStridedFactory>(
+        m, "_erfcinv", "", impl::erfcinv_contig_dispatch_vector,
+        impl::erfcinv_strided_dispatch_vector);
 }
 } // namespace dpnp::extensions::ufunc
diff --git a/dpnp/backend/extensions/vm/erf_funcs.cpp b/dpnp/backend/extensions/vm/erf_funcs.cpp
index c443a7737750..4e84403eb061 100644
--- a/dpnp/backend/extensions/vm/erf_funcs.cpp
+++ b/dpnp/backend/extensions/vm/erf_funcs.cpp
@@ -134,6 +134,8 @@ using ew_cmn_ns::unary_contig_impl_fn_ptr_t;
 MACRO_DEFINE_IMPL(erf, Erf);
 MACRO_DEFINE_IMPL(erfc, Erfc);
 MACRO_DEFINE_IMPL(erfcx, Erfcx);
+MACRO_DEFINE_IMPL(erfinv, Erfinv);
+MACRO_DEFINE_IMPL(erfcinv, Erfcinv);
 
 template <template <typename fnT, typename T> typename factoryT>
 static void populate(py::module_ m,
@@ -194,5 +196,17 @@ void init_erf_funcs(py::module_ m)
         "Call `erfcx` function from OneMKL VM library to compute the scaled "
         "complementary error function value of vector elements",
         impl::erfcx_contig_dispatch_vector);
+
+    impl::populate<impl::ErfinvContigFactory>(
+        m, "_erfinv",
+        "Call `erfinv` function from OneMKL VM library to compute the inverse "
+        "of the error function value of vector elements",
+        impl::erfinv_contig_dispatch_vector);
+
+    impl::populate<impl::ErfcinvContigFactory>(
+        m, "_erfcinv",
+        "Call `erfcinv` function from OneMKL VM library to compute the inverse "
+        "of the complementary error function value of vector elements",
+        impl::erfcinv_contig_dispatch_vector);
 }
 } // namespace dpnp::extensions::vm
diff --git a/dpnp/backend/kernels/elementwise_functions/erf.hpp b/dpnp/backend/kernels/elementwise_functions/erf.hpp
index 70a763ff396e..4ead1174a0c1 100644
--- a/dpnp/backend/kernels/elementwise_functions/erf.hpp
+++ b/dpnp/backend/kernels/elementwise_functions/erf.hpp
@@ -43,6 +43,7 @@
 #include <sycl/ext/intel/math.hpp>
 #else
 #include "erfcx.hpp"
+#include "erfinv.hpp"
 #endif
 
 namespace dpnp::kernels::erfs
@@ -85,13 +86,15 @@ struct BaseFunctor
     template <typename ArgT, typename ResT>                                    \
     using __f_name__##Functor = BaseFunctor<__f_name__##Op, ArgT, ResT>;
 
-MACRO_DEFINE_FUNCTOR(sycl::erf, Erf);
-MACRO_DEFINE_FUNCTOR(sycl::erfc, Erfc);
-MACRO_DEFINE_FUNCTOR(
 #if defined(__SYCL_EXT_INTEL_MATH_SUPPORT)
-    sycl::ext::intel::math::erfcx,
+using namespace sycl::ext::intel::math;
 #else
-    impl::erfcx,
+using namespace impl;
 #endif
-    Erfcx);
+
+MACRO_DEFINE_FUNCTOR(sycl::erf, Erf);
+MACRO_DEFINE_FUNCTOR(sycl::erfc, Erfc);
+MACRO_DEFINE_FUNCTOR(erfcx, Erfcx);
+MACRO_DEFINE_FUNCTOR(erfinv, Erfinv);
+MACRO_DEFINE_FUNCTOR(erfcinv, Erfcinv);
 } // namespace dpnp::kernels::erfs
diff --git a/dpnp/backend/kernels/elementwise_functions/erfcx.hpp b/dpnp/backend/kernels/elementwise_functions/erfcx.hpp
index 9fbb8aeac566..1dd3205d3484 100644
--- a/dpnp/backend/kernels/elementwise_functions/erfcx.hpp
+++ b/dpnp/backend/kernels/elementwise_functions/erfcx.hpp
@@ -1622,19 +1622,17 @@ For x < 0, we use the relationship erfcx(-x) = 2 exp(x^2) - erfc(x), with the
 usual checks for overflow etcetera.
 */
 template <typename Tp>
-Tp erfcx(Tp x)
+inline Tp erfcx(Tp x)
 {
     static_assert(std::is_floating_point_v<Tp>,
                   "erfcx requires a floating-point type");
 
     if (x >= 0) {
-        if (x > 50) // continued-fraction expansion is faster
-        {
+        if (x > 50) { // continued-fraction expansion is faster
             // 1/sqrt(pi)
             constexpr Tp inv_sqrtpi = 0.564189583547756286948079451560772586L;
 
-            if (x > 5e7) // 1-term expansion, important to avoid overflow
-            {
+            if (x > 5e7) { // 1-term expansion, important to avoid overflow
                 return inv_sqrtpi / x;
             }
 
diff --git a/dpnp/backend/kernels/elementwise_functions/erfinv.hpp b/dpnp/backend/kernels/elementwise_functions/erfinv.hpp
new file mode 100644
index 000000000000..80409d6c14a4
--- /dev/null
+++ b/dpnp/backend/kernels/elementwise_functions/erfinv.hpp
@@ -0,0 +1,219 @@
+//*****************************************************************************
+// Copyright (c) 2025, 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.
+//*****************************************************************************
+
+#pragma once
+
+#include <limits>
+#include <sycl/sycl.hpp>
+
+namespace dpnp::kernels::erfs::impl
+{
+template <typename Tp>
+inline Tp polevl(Tp x, const Tp *coeff, int i)
+{
+    Tp p = *coeff++;
+
+    do {
+        p = p * x + *coeff++;
+    } while (--i);
+    return p;
+}
+
+template <typename Tp>
+inline Tp p1evl(Tp x, const Tp *coeff, int i)
+{
+    Tp p = x + *coeff++;
+
+    while (--i) {
+        p = p * x + *coeff++;
+    }
+    return p;
+}
+
+template <typename Tp>
+inline Tp ndtri(Tp y0)
+{
+    Tp y;
+    int code = 1;
+
+    if (y0 == 0.0) {
+        return -std::numeric_limits<Tp>::infinity();
+    }
+    else if (y0 == 1.0) {
+        return std::numeric_limits<Tp>::infinity();
+    }
+    else if (y0 < 0.0 || y0 > 1.0) {
+        return std::numeric_limits<Tp>::quiet_NaN();
+    }
+
+    // exp(-2)
+    constexpr Tp exp_minus2 = 0.13533528323661269189399949497248L;
+    if (y0 > (1.0 - exp_minus2)) {
+        y = 1.0 - y0;
+        code = 0;
+    }
+    else {
+        y = y0;
+    }
+
+    if (y > exp_minus2) {
+        // sqrt(2*pi)
+        constexpr Tp root_2_pi = 2.50662827463100050241576528481105L;
+
+        // approximation for 0 <= |y - 0.5| <= 3/8
+        constexpr Tp p[] = {
+            -5.99633501014107895267E1, 9.80010754185999661536E1,
+            -5.66762857469070293439E1, 1.39312609387279679503E1,
+            -1.23916583867381258016E0,
+        };
+        constexpr Tp q[] = {
+            1.95448858338141759834E0, 4.67627912898881538453E0,
+            8.63602421390890590575E1, -2.25462687854119370527E2,
+            2.00260212380060660359E2, -8.20372256168333339912E1,
+            1.59056225126211695515E1, -1.18331621121330003142E0,
+        };
+
+        y -= 0.5;
+        Tp y2 = y * y;
+        Tp x = y + y * (y2 * polevl(y2, p, 4) / p1evl(y2, q, 8));
+        return x * root_2_pi;
+    }
+
+    Tp x = sycl::sqrt(-2.0 * sycl::log(y));
+    Tp x0 = x - sycl::log(x) / x;
+    Tp z = 1.0 / x;
+
+    Tp x1;
+    if (x < 8.0) {
+        // approximation for 2 <= sqrt(-2*log(y)) < 8
+        constexpr Tp p[] = {
+            4.05544892305962419923E0,   3.15251094599893866154E1,
+            5.71628192246421288162E1,   4.40805073893200834700E1,
+            1.46849561928858024014E1,   2.18663306850790267539E0,
+            -1.40256079171354495875E-1, -3.50424626827848203418E-2,
+            -8.57456785154685413611E-4,
+        };
+
+        constexpr Tp q[] = {
+            1.57799883256466749731E1,   4.53907635128879210584E1,
+            4.13172038254672030440E1,   1.50425385692907503408E1,
+            2.50464946208309415979E0,   -1.42182922854787788574E-1,
+            -3.80806407691578277194E-2, -9.33259480895457427372E-4,
+        };
+
+        x1 = z * polevl(z, p, 8) / p1evl(z, q, 8);
+    }
+    else {
+        // approximation for 8 <= sqrt(-2*log(y)) < 64
+        constexpr Tp p[] = {
+            3.23774891776946035970E0,  6.91522889068984211695E0,
+            3.93881025292474443415E0,  1.33303460815807542389E0,
+            2.01485389549179081538E-1, 1.23716634817820021358E-2,
+            3.01581553508235416007E-4, 2.65806974686737550832E-6,
+            6.23974539184983293730E-9,
+        };
+
+        constexpr Tp q[] = {
+            6.02427039364742014255E0,  3.67983563856160859403E0,
+            1.37702099489081330271E0,  2.16236993594496635890E-1,
+            1.34204006088543189037E-2, 3.28014464682127739104E-4,
+            2.89247864745380683936E-6, 6.79019408009981274425E-9,
+        };
+
+        x1 = z * polevl(z, p, 8) / p1evl(z, q, 8);
+    }
+
+    x = x0 - x1;
+    if (code != 0) {
+        x = -x;
+    }
+    return x;
+}
+
+template <typename Tp>
+inline Tp erfinv(Tp y)
+{
+    static_assert(std::is_floating_point_v<Tp>,
+                  "erfinv requires a floating-point type");
+
+    constexpr Tp lower = -1;
+    constexpr Tp upper = 1;
+
+    constexpr Tp thresh = 1e-7;
+
+    // For small arguments, use the Taylor expansion.
+    // Otherwise, y + 1 loses precision for |y| << 1.
+    if ((-thresh < y) && (y < thresh)) {
+        // 2/sqrt(pi)
+        constexpr Tp inv_sqrtpi = 1.1283791670955125738961589031215452L;
+        return y / inv_sqrtpi;
+    }
+
+    if ((lower < y) && (y < upper)) {
+        // 1/sqrt(2)
+        constexpr Tp one_div_root_2 = 0.7071067811865475244008443621048490L;
+        return ndtri(0.5 * (y + 1)) * one_div_root_2;
+    }
+    else if (y == lower) {
+        return -std::numeric_limits<Tp>::infinity();
+    }
+    else if (y == upper) {
+        return std::numeric_limits<Tp>::infinity();
+    }
+    else if (sycl::isnan(y)) {
+        return y;
+    }
+    return std::numeric_limits<Tp>::quiet_NaN();
+}
+
+template <typename Tp>
+inline Tp erfcinv(Tp y)
+{
+    static_assert(std::is_floating_point_v<Tp>,
+                  "erfcinv requires a floating-point type");
+
+    constexpr Tp lower = 0;
+    constexpr Tp upper = 2;
+
+    if ((lower < y) && (y < upper)) {
+        // 1/sqrt(2)
+        constexpr Tp one_div_root_2 = 0.7071067811865475244008443621048490L;
+        return -ndtri(0.5 * y) * one_div_root_2;
+    }
+    else if (y == lower) {
+        return std::numeric_limits<Tp>::infinity();
+    }
+    else if (y == upper) {
+        return -std::numeric_limits<Tp>::infinity();
+    }
+    else if (sycl::isnan(y)) {
+        return y;
+    }
+    return std::numeric_limits<Tp>::quiet_NaN();
+}
+} // namespace dpnp::kernels::erfs::impl
diff --git a/dpnp/scipy/special/__init__.py b/dpnp/scipy/special/__init__.py
index 50c72d5ace3d..ef8ccfe07c25 100644
--- a/dpnp/scipy/special/__init__.py
+++ b/dpnp/scipy/special/__init__.py
@@ -45,11 +45,15 @@
 from ._erf import (
     erf,
     erfc,
+    erfcinv,
     erfcx,
+    erfinv,
 )
 
 __all__ = [
     "erf",
     "erfc",
+    "erfcinv",
     "erfcx",
+    "erfinv",
 ]
diff --git a/dpnp/scipy/special/_erf.py b/dpnp/scipy/special/_erf.py
index 1e817e46a82d..fef621c9ea4a 100644
--- a/dpnp/scipy/special/_erf.py
+++ b/dpnp/scipy/special/_erf.py
@@ -45,7 +45,7 @@
 import dpnp.backend.extensions.ufunc._ufunc_impl as ufi
 from dpnp.dpnp_algo.dpnp_elementwise_common import DPNPUnaryFunc
 
-__all__ = ["erf", "erfc", "erfcx"]
+__all__ = ["erf", "erfc", "erfcinv", "erfcx", "erfinv"]
 
 
 # pylint: disable=too-few-public-methods
@@ -177,6 +177,7 @@ def __call__(self, x, out=None):  # pylint: disable=signature-differs
 
 _ERFCX_DOCSTRING = r"""
 Calculates the scaled complementary error function of a given input array.
+
 It is defined as :math:`\exp(x^2) * \operatorname{erfc}(x)`.
 
 For full documentation refer to :obj:`scipy.special.erfcx`.
@@ -218,3 +219,110 @@ def __call__(self, x, out=None):  # pylint: disable=signature-differs
     mkl_fn_to_call="_mkl_erf_to_call",
     mkl_impl_fn="_erfcx",
 )
+
+_ERFINV_DOCSTRING = r"""
+Calculates the inverse of the Gauss error function of a given input array.
+
+It is defined as :math:`\operatorname{erf}^{-1}(x)`.
+
+For full documentation refer to :obj:`scipy.special.erfinv`.
+
+Parameters
+----------
+x : {dpnp.ndarray, usm_ndarray}
+    Input array, expected to have a real-valued floating-point data type.
+    Domain: [-1, 1].
+out : {dpnp.ndarray, usm_ndarray}, optional
+    Optional output array for the function values.
+
+Returns
+-------
+out : dpnp.ndarray
+    The values of the inverse of the error function at the given points `x`.
+
+See Also
+--------
+:obj:`dpnp.scipy.special.erf` : Gauss error function.
+:obj:`dpnp.scipy.special.erfc` : Complementary error function.
+:obj:`dpnp.scipy.special.erfcx` : Scaled complementary error function.
+:obj:`dpnp.scipy.special.erfcinv` : Inverse of the complementary error function.
+
+Examples
+--------
+>>> import dpnp as np
+>>> x = np.linspace(-1.0, 1.0, num=9)
+>>> y = np.scipy.special.erfinv(x)
+>>> y
+array([       -inf, -0.81341985, -0.47693628, -0.22531206,  0.        ,
+        0.22531206,  0.47693628,  0.81341985,         inf])
+
+Verify that ``erf(erfinv(x))`` is ``x``:
+
+>>> np.scipy.special.erf(y)
+array([-1.  , -0.75, -0.5 , -0.25,  0.  ,  0.25,  0.5 ,  0.75,  1.  ])
+
+"""
+
+erfinv = DPNPErf(
+    "erfinv",
+    ufi._erf_result_type,
+    ufi._erfinv,
+    _ERFINV_DOCSTRING,
+    mkl_fn_to_call="_mkl_erf_to_call",
+    mkl_impl_fn="_erfinv",
+)
+
+_ERFCINV_DOCSTRING = r"""
+Calculates the inverse of the complementary error function of a given input
+array.
+
+It is defined as :math:`\operatorname{erfinv}(1 - x)`.
+
+For full documentation refer to :obj:`scipy.special.erfcinv`.
+
+Parameters
+----------
+x : {dpnp.ndarray, usm_ndarray}
+    Input array, expected to have a real-valued floating-point data type.
+    Domain: [0, 2].
+out : {dpnp.ndarray, usm_ndarray}, optional
+    Optional output array for the function values.
+
+Returns
+-------
+out : dpnp.ndarray
+    The values of the inverse of the complementary error function at the given
+    points `x`.
+
+See Also
+--------
+:obj:`dpnp.scipy.special.erf` : Gauss error function.
+:obj:`dpnp.scipy.special.erfc` : Complementary error function.
+:obj:`dpnp.scipy.special.erfcx` : Scaled complementary error function.
+:obj:`dpnp.scipy.special.erfinv` : Inverse of the error function.
+
+Examples
+--------
+>>> import dpnp as np
+>>> x = np.linspace(0.0, 2.0, num=11)
+>>> y = np.scipy.special.erfcinv(x)
+>>> y
+array([        inf,  0.9061938 ,  0.59511608,  0.37080716,  0.17914345,
+        0.        , -0.17914345, -0.37080716, -0.59511608, -0.9061938 ,
+              -inf])
+
+Verify that ``erfc(erfcinv(x))`` is ``x``:
+
+>>> np.scipy.special.erfc(y)
+array([0. , 0.2, 0.4, 0.6, 0.8, 1. , 1.2, 1.4, 1.6, 1.8, 2. ])
+
+"""
+
+erfcinv = DPNPErf(
+    "erfcinv",
+    ufi._erf_result_type,
+    ufi._erfcinv,
+    _ERFCINV_DOCSTRING,
+    mkl_fn_to_call="_mkl_erf_to_call",
+    mkl_impl_fn="_erfcinv",
+)
diff --git a/dpnp/tests/test_special.py b/dpnp/tests/test_special.py
index 66b62aad65ac..257ee3c5758b 100644
--- a/dpnp/tests/test_special.py
+++ b/dpnp/tests/test_special.py
@@ -63,17 +63,27 @@ def test_complex(self, func, dt):
             getattr(dpnp.scipy.special, func)(x)
 
 
+@with_requires("scipy")
+@pytest.mark.parametrize("inverse", [True, False])
 class TestConsistency:
 
     tol = 8 * dpnp.finfo(dpnp.default_float_type()).resolution
 
-    def _check_variant_func(self, func, other_func, rtol, atol=0):
-        # TODO: replace with dpnp.random.RandomState, once pareto is added
-        rng = numpy.random.RandomState(1234)
-        n = 10000
-        a = rng.pareto(0.02, n) * (2 * rng.randint(0, 2, n) - 1)
-        a = dpnp.array(a)
-        a = a[::-1]
+    def _check_variant_func(
+        self, inverse, func, other_func, rtol, atol=0, data=None
+    ):
+        if data is not None:
+            a = data
+        else:
+            # TODO: replace with dpnp.random.RandomState, once pareto is added
+            rng = numpy.random.RandomState(1234)
+            n = 10000
+            a = rng.pareto(0.02, n) * (2 * rng.randint(0, 2, n) - 1)
+            a = dpnp.array(a)
+
+        if inverse:
+            # select either ufunc or vm implementation (SYCL or oneMKL call)
+            a = a[::-1]
 
         res = other_func(a)
         mask = dpnp.isfinite(res)
@@ -84,17 +94,76 @@ def _check_variant_func(self, func, other_func, rtol, atol=0):
             # calling numpy testing func, because it's more verbose
             assert_allclose(x.asnumpy(), y.asnumpy(), rtol=rtol, atol=atol)
 
-    def test_erfc(self):
+    def test_erfc(self, inverse):
         self._check_variant_func(
+            inverse,
             dpnp.scipy.special.erfc,
             lambda z: 1 - dpnp.scipy.special.erf(z),
             rtol=self.tol,
             atol=self.tol,
         )
 
-    def test_erfcx(self):
+    def test_erfcx(self, inverse):
         self._check_variant_func(
+            inverse,
             dpnp.scipy.special.erfcx,
             lambda z: dpnp.exp(z * z) * dpnp.scipy.special.erfc(z),
             rtol=10 * self.tol,
         )
+
+    def test_erfinv(self, inverse):
+        self._check_variant_func(
+            inverse,
+            dpnp.scipy.special.erfinv,
+            lambda z: dpnp.scipy.special.erfcinv(1 - z),
+            rtol=self.tol,
+            data=dpnp.linspace(-1, 1, 101),
+        )
+
+
+@with_requires("scipy")
+@pytest.mark.parametrize("func", ["erfinv", "erfcinv"])
+class TestInverse:
+
+    @pytest.mark.parametrize(
+        "dt", get_all_dtypes(no_none=True, no_float16=False, no_complex=True)
+    )
+    def test_literal_values(self, func, dt):
+        import scipy.special
+
+        a = numpy.array(
+            [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9], dtype=dt
+        )
+        ia = dpnp.array(a)
+
+        result = getattr(dpnp.scipy.special, func)(ia)
+        expected = getattr(scipy.special, func)(a)
+        assert_dtype_allclose(result, expected)
+
+    def test_bound_values(self, func):
+        import scipy.special
+
+        a = numpy.array([-1, 0, 1, -100, 100, 0, 1, 2, -100, 100])
+        ia = dpnp.array(a)
+
+        result = getattr(dpnp.scipy.special, func)(ia)
+        expected = getattr(scipy.special, func)(a)
+        assert_dtype_allclose(result, expected)
+
+    def test_small_values(self, func):
+        import scipy.special
+
+        a = numpy.array(
+            [1e-20, 1e-15, 1e-14, 1e-10, 1e-8, 0.9e-7, 1.1e-7, 1e-6]
+        )
+        ia = dpnp.array(a)
+
+        result = getattr(dpnp.scipy.special, func)(ia)
+        expected = getattr(scipy.special, func)(a)
+        assert_dtype_allclose(result, expected)
+
+    @pytest.mark.parametrize("dt", get_complex_dtypes())
+    def test_complex(self, func, dt):
+        x = dpnp.empty(5, dtype=dt)
+        with pytest.raises(ValueError):
+            getattr(dpnp.scipy.special, func)(x)
diff --git a/dpnp/tests/test_strides.py b/dpnp/tests/test_strides.py
index 3690ff61ec02..102a20d60c04 100644
--- a/dpnp/tests/test_strides.py
+++ b/dpnp/tests/test_strides.py
@@ -167,12 +167,21 @@ def test_reduce_hypot(dtype, stride):
 
 
 @with_requires("scipy")
-@pytest.mark.parametrize("func", ["erf", "erfc", "erfcx"])
+@pytest.mark.parametrize(
+    "func, low, high",
+    [
+        ("erf", -3, 3),
+        ("erfc", -3, 3),
+        ("erfcx", -3, 3),
+        ("erfinv", -1, 1),
+        ("erfcinv", 0, 2),
+    ],
+)
 @pytest.mark.parametrize("stride", [2, -1, -3])
-def test_erf_funcs(func, stride):
+def test_erf_funcs(func, low, high, stride):
     import scipy.special
 
-    x = generate_random_numpy_array(10)
+    x = generate_random_numpy_array(10, low=low, high=high)
     a, ia = x[::stride], dpnp.array(x)[::stride]
 
     result = getattr(dpnp.scipy.special, func)(ia)
diff --git a/dpnp/tests/test_sycl_queue.py b/dpnp/tests/test_sycl_queue.py
index 3ed3ef17f96e..cf905cc94646 100644
--- a/dpnp/tests/test_sycl_queue.py
+++ b/dpnp/tests/test_sycl_queue.py
@@ -1488,10 +1488,19 @@ def test_interp(device, left, right, period):
     assert_sycl_queue_equal(result.sycl_queue, x.sycl_queue)
 
 
-@pytest.mark.parametrize("func", ["erf", "erfc", "erfcx"])
+@pytest.mark.parametrize(
+    "func, low, high",
+    [
+        ("erf", -3, 3),
+        ("erfc", -3, 3),
+        ("erfcx", -3, 3),
+        ("erfinv", -1, 1),
+        ("erfcinv", 0, 2),
+    ],
+)
 @pytest.mark.parametrize("device", valid_dev, ids=dev_ids)
-def test_erf_funcs(func, device):
-    x = dpnp.linspace(-3, 3, num=5, device=device)
+def test_erf_funcs(func, low, high, device):
+    x = dpnp.linspace(low, high, num=5, device=device)
 
     result = getattr(dpnp.scipy.special, func)(x)
     assert_sycl_queue_equal(result.sycl_queue, x.sycl_queue)
diff --git a/dpnp/tests/test_usm_type.py b/dpnp/tests/test_usm_type.py
index 97ceb0220e53..453ffea735e1 100644
--- a/dpnp/tests/test_usm_type.py
+++ b/dpnp/tests/test_usm_type.py
@@ -1296,10 +1296,19 @@ def test_choose(usm_type_x, usm_type_ind):
     assert z.usm_type == du.get_coerced_usm_type([usm_type_x, usm_type_ind])
 
 
-@pytest.mark.parametrize("func", ["erf", "erfc", "erfcx"])
+@pytest.mark.parametrize(
+    "func, low, high",
+    [
+        ("erf", -3, 3),
+        ("erfc", -3, 3),
+        ("erfcx", -3, 3),
+        ("erfinv", -1, 1),
+        ("erfcinv", 0, 2),
+    ],
+)
 @pytest.mark.parametrize("usm_type", list_of_usm_types)
-def test_erf_funcs(func, usm_type):
-    x = dpnp.linspace(-3, 3, num=5, usm_type=usm_type)
+def test_erf_funcs(func, low, high, usm_type):
+    x = dpnp.linspace(low, high, num=5, usm_type=usm_type)
     y = getattr(dpnp.scipy.special, func)(x)
     assert x.usm_type == y.usm_type == usm_type
 
diff --git a/dpnp/tests/third_party/cupyx/scipy_tests/special_tests/test_erf.py b/dpnp/tests/third_party/cupyx/scipy_tests/special_tests/test_erf.py
index 13d1523ff781..1bdee3dbe54c 100644
--- a/dpnp/tests/third_party/cupyx/scipy_tests/special_tests/test_erf.py
+++ b/dpnp/tests/third_party/cupyx/scipy_tests/special_tests/test_erf.py
@@ -31,14 +31,12 @@ def test_erfc(self):
     def test_erfcx(self):
         self.check_unary("erfcx")
 
-    @pytest.mark.skip("erfinv() is not supported yet")
     def test_erfinv(self):
         self.check_unary("erfinv")
         self.check_unary_random("erfinv", scale=2, offset=-1)
         self.check_unary_boundary("erfinv", boundary=-1)
         self.check_unary_boundary("erfinv", boundary=1)
 
-    @pytest.mark.skip("erfcinv() is not supported yet")
     def test_erfcinv(self):
         self.check_unary("erfcinv")
         self.check_unary_random("erfcinv", scale=2, offset=0)
@@ -75,42 +73,40 @@ def check_unary_boundary(self, name, xp, scp, dtype, boundary):
         a = xp.array(a, dtype=dtype)
         return getattr(scp.special, name)(a)
 
-    @pytest.mark.skip("erfinv() is not supported yet")
     @testing.with_requires("scipy>=1.4.0")
     @testing.for_dtypes(["f", "d"])
     def test_erfinv_behavior(self, dtype):
         a = cupy.empty((1,), dtype=dtype)
 
         a[:] = 1.0 + 1e-6
-        a = cupyx.scipy.special.erfinv(a)
+        a = cupy.scipy.special.erfinv(a)
         assert cupy.isnan(a)
         a[:] = -1.0 - 1e-6
-        a = cupyx.scipy.special.erfinv(a)
+        a = cupy.scipy.special.erfinv(a)
         assert cupy.isnan(a)
         a[:] = 1.0
-        a = cupyx.scipy.special.erfinv(a)
+        a = cupy.scipy.special.erfinv(a)
         assert numpy.isposinf(cupy.asnumpy(a))
         a[:] = -1.0
-        a = cupyx.scipy.special.erfinv(a)
+        a = cupy.scipy.special.erfinv(a)
         assert numpy.isneginf(cupy.asnumpy(a))
 
-    @pytest.mark.skip("erfcinv() is not supported yet")
     @testing.with_requires("scipy>=1.4.0")
     @testing.for_dtypes(["f", "d"])
     def test_erfcinv_behavior(self, dtype):
         a = cupy.empty((1,), dtype=dtype)
 
         a[:] = 2.0 + 1e-6
-        a = cupyx.scipy.special.erfcinv(a)
+        a = cupy.scipy.special.erfcinv(a)
         assert cupy.isnan(a)
         a[:] = 0.0 - 1e-6
-        a = cupyx.scipy.special.erfcinv(a)
+        a = cupy.scipy.special.erfcinv(a)
         assert cupy.isnan(a)
         a[:] = 0.0
-        a = cupyx.scipy.special.erfcinv(a)
+        a = cupy.scipy.special.erfcinv(a)
         assert numpy.isposinf(cupy.asnumpy(a))
         a[:] = 2.0
-        a = cupyx.scipy.special.erfcinv(a)
+        a = cupy.scipy.special.erfcinv(a)
         assert numpy.isneginf(cupy.asnumpy(a))