ENH: add sinc (#20)

* ENH: add `sinc`

* TST: sinc: add tests

* DOC: sinc: tweak

* DOC: sinc: more tweaks

* TST: sinc: add 3D test

* improve `where` call

Co-authored-by: jakirkham <[email protected]>

* Apply suggestions from code review

* DOC: sinc: add `xp` param

* Update _funcs.py

Co-authored-by: jakirkham <[email protected]>

* appease formatter

---------

Co-authored-by: jakirkham <[email protected]>

docs/api-reference.md

@@ -10,4 +10,5 @@
     cov
     expand_dims
     kron
+    sinc
 ```

src/array_api_extra/__init__.py

@@ -1,7 +1,7 @@
 from __future__ import annotations
 
-from ._funcs import atleast_nd, cov, expand_dims, kron
+from ._funcs import atleast_nd, cov, expand_dims, kron, sinc
 
 __version__ = "0.1.2.dev0"
 
-__all__ = ["__version__", "atleast_nd", "cov", "expand_dims", "kron"]
+__all__ = ["__version__", "atleast_nd", "cov", "expand_dims", "kron", "sinc"]

src/array_api_extra/_funcs.py

@@ -6,7 +6,7 @@
 if TYPE_CHECKING:
     from ._typing import Array, ModuleType
 
-__all__ = ["atleast_nd", "cov", "expand_dims", "kron"]
+__all__ = ["atleast_nd", "cov", "expand_dims", "kron", "sinc"]
 
 
 def atleast_nd(x: Array, /, *, ndim: int, xp: ModuleType) -> Array:
@@ -348,3 +348,88 @@ def kron(a: Array, b: Array, /, *, xp: ModuleType) -> Array:
     a_shape = xp.asarray(a_shape)
     b_shape = xp.asarray(b_shape)
     return xp.reshape(result, tuple(xp.multiply(a_shape, b_shape)))
+
+
+def sinc(x: Array, /, *, xp: ModuleType) -> Array:
+    r"""
+    Return the normalized sinc function.
+
+    The sinc function is equal to :math:`\sin(\pi x)/(\pi x)` for any argument
+    :math:`x\ne 0`. ``sinc(0)`` takes the limit value 1, making ``sinc`` not
+    only everywhere continuous but also infinitely differentiable.
+
+    .. note::
+
+        Note the normalization factor of ``pi`` used in the definition.
+        This is the most commonly used definition in signal processing.
+        Use ``sinc(x / xp.pi)`` to obtain the unnormalized sinc function
+        :math:`\sin(x)/x` that is more common in mathematics.
+
+    Parameters
+    ----------
+    x : array
+        Array (possibly multi-dimensional) of values for which to calculate
+        ``sinc(x)``. Must have a real floating point dtype.
+    xp : array_namespace
+        The standard-compatible namespace for `x`.
+
+    Returns
+    -------
+    res : array
+        ``sinc(x)`` calculated elementwise, which has the same shape as the input.
+
+    Notes
+    -----
+    The name sinc is short for "sine cardinal" or "sinus cardinalis".
+
+    The sinc function is used in various signal processing applications,
+    including in anti-aliasing, in the construction of a Lanczos resampling
+    filter, and in interpolation.
+
+    For bandlimited interpolation of discrete-time signals, the ideal
+    interpolation kernel is proportional to the sinc function.
+
+    References
+    ----------
+    .. [1] Weisstein, Eric W. "Sinc Function." From MathWorld--A Wolfram Web
+           Resource. https://mathworld.wolfram.com/SincFunction.html
+    .. [2] Wikipedia, "Sinc function",
+           https://en.wikipedia.org/wiki/Sinc_function
+
+    Examples
+    --------
+    >>> import array_api_strict as xp
+    >>> import array_api_extra as xpx
+    >>> x = xp.linspace(-4, 4, 41)
+    >>> xpx.sinc(x, xp=xp)
+    Array([-3.89817183e-17, -4.92362781e-02,
+           -8.40918587e-02, -8.90384387e-02,
+           -5.84680802e-02,  3.89817183e-17,
+            6.68206631e-02,  1.16434881e-01,
+            1.26137788e-01,  8.50444803e-02,
+           -3.89817183e-17, -1.03943254e-01,
+           -1.89206682e-01, -2.16236208e-01,
+           -1.55914881e-01,  3.89817183e-17,
+            2.33872321e-01,  5.04551152e-01,
+            7.56826729e-01,  9.35489284e-01,
+            1.00000000e+00,  9.35489284e-01,
+            7.56826729e-01,  5.04551152e-01,
+            2.33872321e-01,  3.89817183e-17,
+           -1.55914881e-01, -2.16236208e-01,
+           -1.89206682e-01, -1.03943254e-01,
+           -3.89817183e-17,  8.50444803e-02,
+            1.26137788e-01,  1.16434881e-01,
+            6.68206631e-02,  3.89817183e-17,
+           -5.84680802e-02, -8.90384387e-02,
+           -8.40918587e-02, -4.92362781e-02,
+           -3.89817183e-17], dtype=array_api_strict.float64)
+
+    """
+    if not xp.isdtype(x.dtype, "real floating"):
+        err_msg = "`x` must have a real floating data type."
+        raise ValueError(err_msg)
+    # no scalars in `where` - array-api#807
+    y = xp.pi * xp.where(
+        x, x, xp.asarray(xp.finfo(x.dtype).smallest_normal, dtype=x.dtype)
+    )
+    return xp.sin(y) / y

tests/test_funcs.py

@@ -9,7 +9,7 @@
 import pytest
 from numpy.testing import assert_allclose, assert_array_equal, assert_equal
 
-from array_api_extra import atleast_nd, cov, expand_dims, kron
+from array_api_extra import atleast_nd, cov, expand_dims, kron, sinc
 
 if TYPE_CHECKING:
     Array = Any  # To be changed to a Protocol later (see array-api#589)
@@ -224,3 +224,22 @@ def test_positive_negative_repeated(self):
         a = xp.empty((2, 3, 4, 5))
         with pytest.raises(ValueError, match="Duplicate dimensions"):
             expand_dims(a, axis=(3, -3), xp=xp)
+
+
+class TestSinc:
+    def test_simple(self):
+        assert_array_equal(sinc(xp.asarray(0.0), xp=xp), xp.asarray(1.0))
+        w = sinc(xp.linspace(-1, 1, 100), xp=xp)
+        # check symmetry
+        assert_allclose(w, xp.flip(w, axis=0))
+
+    @pytest.mark.parametrize("x", [0, 1 + 3j])
+    def test_dtype(self, x):
+        with pytest.raises(ValueError, match="real floating data type"):
+            sinc(xp.asarray(x), xp=xp)
+
+    def test_3d(self):
+        x = xp.reshape(xp.arange(18, dtype=xp.float64), (3, 3, 2))
+        expected = xp.zeros((3, 3, 2))
+        expected[0, 0, 0] = 1.0
+        assert_allclose(sinc(x, xp=xp), expected, atol=1e-15)