2.0

19 files changed, 3164 insertions, 0 deletions

diff --git a/.venv/lib/python3.12/site-packages/numpy/fft/__init__.py b/.venv/lib/python3.12/site-packages/numpy/fft/__init__.py
new file mode 100644
index 0000000..55f7320
--- /dev/null
+++ b/.venv/lib/python3.12/site-packages/numpy/fft/__init__.py
@@ -0,0 +1,215 @@
+"""
+Discrete Fourier Transform
+==========================
+
+.. currentmodule:: numpy.fft
+
+The SciPy module `scipy.fft` is a more comprehensive superset
+of `numpy.fft`, which includes only a basic set of routines.
+
+Standard FFTs
+-------------
+
+.. autosummary::
+   :toctree: generated/
+
+   fft       Discrete Fourier transform.
+   ifft      Inverse discrete Fourier transform.
+   fft2      Discrete Fourier transform in two dimensions.
+   ifft2     Inverse discrete Fourier transform in two dimensions.
+   fftn      Discrete Fourier transform in N-dimensions.
+   ifftn     Inverse discrete Fourier transform in N dimensions.
+
+Real FFTs
+---------
+
+.. autosummary::
+   :toctree: generated/
+
+   rfft      Real discrete Fourier transform.
+   irfft     Inverse real discrete Fourier transform.
+   rfft2     Real discrete Fourier transform in two dimensions.
+   irfft2    Inverse real discrete Fourier transform in two dimensions.
+   rfftn     Real discrete Fourier transform in N dimensions.
+   irfftn    Inverse real discrete Fourier transform in N dimensions.
+
+Hermitian FFTs
+--------------
+
+.. autosummary::
+   :toctree: generated/
+
+   hfft      Hermitian discrete Fourier transform.
+   ihfft     Inverse Hermitian discrete Fourier transform.
+
+Helper routines
+---------------
+
+.. autosummary::
+   :toctree: generated/
+
+   fftfreq   Discrete Fourier Transform sample frequencies.
+   rfftfreq  DFT sample frequencies (for usage with rfft, irfft).
+   fftshift  Shift zero-frequency component to center of spectrum.
+   ifftshift Inverse of fftshift.
+
+
+Background information
+----------------------
+
+Fourier analysis is fundamentally a method for expressing a function as a
+sum of periodic components, and for recovering the function from those
+components.  When both the function and its Fourier transform are
+replaced with discretized counterparts, it is called the discrete Fourier
+transform (DFT).  The DFT has become a mainstay of numerical computing in
+part because of a very fast algorithm for computing it, called the Fast
+Fourier Transform (FFT), which was known to Gauss (1805) and was brought
+to light in its current form by Cooley and Tukey [CT]_.  Press et al. [NR]_
+provide an accessible introduction to Fourier analysis and its
+applications.
+
+Because the discrete Fourier transform separates its input into
+components that contribute at discrete frequencies, it has a great number
+of applications in digital signal processing, e.g., for filtering, and in
+this context the discretized input to the transform is customarily
+referred to as a *signal*, which exists in the *time domain*.  The output
+is called a *spectrum* or *transform* and exists in the *frequency
+domain*.
+
+Implementation details
+----------------------
+
+There are many ways to define the DFT, varying in the sign of the
+exponent, normalization, etc.  In this implementation, the DFT is defined
+as
+
+.. math::
+   A_k =  \\sum_{m=0}^{n-1} a_m \\exp\\left\\{-2\\pi i{mk \\over n}\\right\\}
+   \\qquad k = 0,\\ldots,n-1.
+
+The DFT is in general defined for complex inputs and outputs, and a
+single-frequency component at linear frequency :math:`f` is
+represented by a complex exponential
+:math:`a_m = \\exp\\{2\\pi i\\,f m\\Delta t\\}`, where :math:`\\Delta t`
+is the sampling interval.
+
+The values in the result follow so-called "standard" order: If ``A =
+fft(a, n)``, then ``A[0]`` contains the zero-frequency term (the sum of
+the signal), which is always purely real for real inputs. Then ``A[1:n/2]``
+contains the positive-frequency terms, and ``A[n/2+1:]`` contains the
+negative-frequency terms, in order of decreasingly negative frequency.
+For an even number of input points, ``A[n/2]`` represents both positive and
+negative Nyquist frequency, and is also purely real for real input.  For
+an odd number of input points, ``A[(n-1)/2]`` contains the largest positive
+frequency, while ``A[(n+1)/2]`` contains the largest negative frequency.
+The routine ``np.fft.fftfreq(n)`` returns an array giving the frequencies
+of corresponding elements in the output.  The routine
+``np.fft.fftshift(A)`` shifts transforms and their frequencies to put the
+zero-frequency components in the middle, and ``np.fft.ifftshift(A)`` undoes
+that shift.
+
+When the input `a` is a time-domain signal and ``A = fft(a)``, ``np.abs(A)``
+is its amplitude spectrum and ``np.abs(A)**2`` is its power spectrum.
+The phase spectrum is obtained by ``np.angle(A)``.
+
+The inverse DFT is defined as
+
+.. math::
+   a_m = \\frac{1}{n}\\sum_{k=0}^{n-1}A_k\\exp\\left\\{2\\pi i{mk\\over n}\\right\\}
+   \\qquad m = 0,\\ldots,n-1.
+
+It differs from the forward transform by the sign of the exponential
+argument and the default normalization by :math:`1/n`.
+
+Type Promotion
+--------------
+
+`numpy.fft` promotes ``float32`` and ``complex64`` arrays to ``float64`` and
+``complex128`` arrays respectively. For an FFT implementation that does not
+promote input arrays, see `scipy.fftpack`.
+
+Normalization
+-------------
+
+The argument ``norm`` indicates which direction of the pair of direct/inverse
+transforms is scaled and with what normalization factor.
+The default normalization (``"backward"``) has the direct (forward) transforms
+unscaled and the inverse (backward) transforms scaled by :math:`1/n`. It is
+possible to obtain unitary transforms by setting the keyword argument ``norm``
+to ``"ortho"`` so that both direct and inverse transforms are scaled by
+:math:`1/\\sqrt{n}`. Finally, setting the keyword argument ``norm`` to
+``"forward"`` has the direct transforms scaled by :math:`1/n` and the inverse
+transforms unscaled (i.e. exactly opposite to the default ``"backward"``).
+`None` is an alias of the default option ``"backward"`` for backward
+compatibility.
+
+Real and Hermitian transforms
+-----------------------------
+
+When the input is purely real, its transform is Hermitian, i.e., the
+component at frequency :math:`f_k` is the complex conjugate of the
+component at frequency :math:`-f_k`, which means that for real
+inputs there is no information in the negative frequency components that
+is not already available from the positive frequency components.
+The family of `rfft` functions is
+designed to operate on real inputs, and exploits this symmetry by
+computing only the positive frequency components, up to and including the
+Nyquist frequency.  Thus, ``n`` input points produce ``n/2+1`` complex
+output points.  The inverses of this family assumes the same symmetry of
+its input, and for an output of ``n`` points uses ``n/2+1`` input points.
+
+Correspondingly, when the spectrum is purely real, the signal is
+Hermitian.  The `hfft` family of functions exploits this symmetry by
+using ``n/2+1`` complex points in the input (time) domain for ``n`` real
+points in the frequency domain.
+
+In higher dimensions, FFTs are used, e.g., for image analysis and
+filtering.  The computational efficiency of the FFT means that it can
+also be a faster way to compute large convolutions, using the property
+that a convolution in the time domain is equivalent to a point-by-point
+multiplication in the frequency domain.
+
+Higher dimensions
+-----------------
+
+In two dimensions, the DFT is defined as
+
+.. math::
+   A_{kl} =  \\sum_{m=0}^{M-1} \\sum_{n=0}^{N-1}
+   a_{mn}\\exp\\left\\{-2\\pi i \\left({mk\\over M}+{nl\\over N}\\right)\\right\\}
+   \\qquad k = 0, \\ldots, M-1;\\quad l = 0, \\ldots, N-1,
+
+which extends in the obvious way to higher dimensions, and the inverses
+in higher dimensions also extend in the same way.
+
+References
+----------
+
+.. [CT] Cooley, James W., and John W. Tukey, 1965, "An algorithm for the
+        machine calculation of complex Fourier series," *Math. Comput.*
+        19: 297-301.
+
+.. [NR] Press, W., Teukolsky, S., Vetterline, W.T., and Flannery, B.P.,
+        2007, *Numerical Recipes: The Art of Scientific Computing*, ch.
+        12-13.  Cambridge Univ. Press, Cambridge, UK.
+
+Examples
+--------
+
+For examples, see the various functions.
+
+"""
+
+# TODO: `numpy.fft.helper`` was deprecated in NumPy 2.0. It should
+# be deleted once downstream libraries move to `numpy.fft`.
+from . import _helper, _pocketfft, helper
+from ._helper import *
+from ._pocketfft import *
+
+__all__ = _pocketfft.__all__.copy()  # noqa: PLE0605
+__all__ += _helper.__all__
+
+from numpy._pytesttester import PytestTester
+
+test = PytestTester(__name__)
+del PytestTester
diff --git a/.venv/lib/python3.12/site-packages/numpy/fft/__init__.pyi b/.venv/lib/python3.12/site-packages/numpy/fft/__init__.pyi
new file mode 100644
index 0000000..54d0ea8
--- /dev/null
+++ b/.venv/lib/python3.12/site-packages/numpy/fft/__init__.pyi
@@ -0,0 +1,43 @@
+from ._helper import (
+    fftfreq,
+    fftshift,
+    ifftshift,
+    rfftfreq,
+)
+from ._pocketfft import (
+    fft,
+    fft2,
+    fftn,
+    hfft,
+    ifft,
+    ifft2,
+    ifftn,
+    ihfft,
+    irfft,
+    irfft2,
+    irfftn,
+    rfft,
+    rfft2,
+    rfftn,
+)
+
+__all__ = [
+    "fft",
+    "ifft",
+    "rfft",
+    "irfft",
+    "hfft",
+    "ihfft",
+    "rfftn",
+    "irfftn",
+    "rfft2",
+    "irfft2",
+    "fft2",
+    "ifft2",
+    "fftn",
+    "ifftn",
+    "fftshift",
+    "ifftshift",
+    "fftfreq",
+    "rfftfreq",
+]
diff --git a/.venv/lib/python3.12/site-packages/numpy/fft/__pycache__/__init__.cpython-312.pyc b/.venv/lib/python3.12/site-packages/numpy/fft/__pycache__/__init__.cpython-312.pyc
new file mode 100644
index 0000000..83e5e34
--- /dev/null
+++ b/.venv/lib/python3.12/site-packages/numpy/fft/__pycache__/__init__.cpython-312.pyc
diff --git a/.venv/lib/python3.12/site-packages/numpy/fft/__pycache__/_helper.cpython-312.pyc b/.venv/lib/python3.12/site-packages/numpy/fft/__pycache__/_helper.cpython-312.pyc
new file mode 100644
index 0000000..d39f4fe
--- /dev/null
+++ b/.venv/lib/python3.12/site-packages/numpy/fft/__pycache__/_helper.cpython-312.pyc
diff --git a/.venv/lib/python3.12/site-packages/numpy/fft/__pycache__/_pocketfft.cpython-312.pyc b/.venv/lib/python3.12/site-packages/numpy/fft/__pycache__/_pocketfft.cpython-312.pyc
new file mode 100644
index 0000000..0705dbb
--- /dev/null
+++ b/.venv/lib/python3.12/site-packages/numpy/fft/__pycache__/_pocketfft.cpython-312.pyc
diff --git a/.venv/lib/python3.12/site-packages/numpy/fft/__pycache__/helper.cpython-312.pyc b/.venv/lib/python3.12/site-packages/numpy/fft/__pycache__/helper.cpython-312.pyc
new file mode 100644
index 0000000..ec32c53
--- /dev/null
+++ b/.venv/lib/python3.12/site-packages/numpy/fft/__pycache__/helper.cpython-312.pyc
diff --git a/.venv/lib/python3.12/site-packages/numpy/fft/_helper.py b/.venv/lib/python3.12/site-packages/numpy/fft/_helper.py
new file mode 100644
index 0000000..77adeac
--- /dev/null
+++ b/.venv/lib/python3.12/site-packages/numpy/fft/_helper.py
@@ -0,0 +1,235 @@
+"""
+Discrete Fourier Transforms - _helper.py
+
+"""
+from numpy._core import arange, asarray, empty, integer, roll
+from numpy._core.overrides import array_function_dispatch, set_module
+
+# Created by Pearu Peterson, September 2002
+
+__all__ = ['fftshift', 'ifftshift', 'fftfreq', 'rfftfreq']
+
+integer_types = (int, integer)
+
+
+def _fftshift_dispatcher(x, axes=None):
+    return (x,)
+
+
+@array_function_dispatch(_fftshift_dispatcher, module='numpy.fft')
+def fftshift(x, axes=None):
+    """
+    Shift the zero-frequency component to the center of the spectrum.
+
+    This function swaps half-spaces for all axes listed (defaults to all).
+    Note that ``y[0]`` is the Nyquist component only if ``len(x)`` is even.
+
+    Parameters
+    ----------
+    x : array_like
+        Input array.
+    axes : int or shape tuple, optional
+        Axes over which to shift.  Default is None, which shifts all axes.
+
+    Returns
+    -------
+    y : ndarray
+        The shifted array.
+
+    See Also
+    --------
+    ifftshift : The inverse of `fftshift`.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> freqs = np.fft.fftfreq(10, 0.1)
+    >>> freqs
+    array([ 0.,  1.,  2., ..., -3., -2., -1.])
+    >>> np.fft.fftshift(freqs)
+    array([-5., -4., -3., -2., -1.,  0.,  1.,  2.,  3.,  4.])
+
+    Shift the zero-frequency component only along the second axis:
+
+    >>> freqs = np.fft.fftfreq(9, d=1./9).reshape(3, 3)
+    >>> freqs
+    array([[ 0.,  1.,  2.],
+           [ 3.,  4., -4.],
+           [-3., -2., -1.]])
+    >>> np.fft.fftshift(freqs, axes=(1,))
+    array([[ 2.,  0.,  1.],
+           [-4.,  3.,  4.],
+           [-1., -3., -2.]])
+
+    """
+    x = asarray(x)
+    if axes is None:
+        axes = tuple(range(x.ndim))
+        shift = [dim // 2 for dim in x.shape]
+    elif isinstance(axes, integer_types):
+        shift = x.shape[axes] // 2
+    else:
+        shift = [x.shape[ax] // 2 for ax in axes]
+
+    return roll(x, shift, axes)
+
+
+@array_function_dispatch(_fftshift_dispatcher, module='numpy.fft')
+def ifftshift(x, axes=None):
+    """
+    The inverse of `fftshift`. Although identical for even-length `x`, the
+    functions differ by one sample for odd-length `x`.
+
+    Parameters
+    ----------
+    x : array_like
+        Input array.
+    axes : int or shape tuple, optional
+        Axes over which to calculate.  Defaults to None, which shifts all axes.
+
+    Returns
+    -------
+    y : ndarray
+        The shifted array.
+
+    See Also
+    --------
+    fftshift : Shift zero-frequency component to the center of the spectrum.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> freqs = np.fft.fftfreq(9, d=1./9).reshape(3, 3)
+    >>> freqs
+    array([[ 0.,  1.,  2.],
+           [ 3.,  4., -4.],
+           [-3., -2., -1.]])
+    >>> np.fft.ifftshift(np.fft.fftshift(freqs))
+    array([[ 0.,  1.,  2.],
+           [ 3.,  4., -4.],
+           [-3., -2., -1.]])
+
+    """
+    x = asarray(x)
+    if axes is None:
+        axes = tuple(range(x.ndim))
+        shift = [-(dim // 2) for dim in x.shape]
+    elif isinstance(axes, integer_types):
+        shift = -(x.shape[axes] // 2)
+    else:
+        shift = [-(x.shape[ax] // 2) for ax in axes]
+
+    return roll(x, shift, axes)
+
+
+@set_module('numpy.fft')
+def fftfreq(n, d=1.0, device=None):
+    """
+    Return the Discrete Fourier Transform sample frequencies.
+
+    The returned float array `f` contains the frequency bin centers in cycles
+    per unit of the sample spacing (with zero at the start).  For instance, if
+    the sample spacing is in seconds, then the frequency unit is cycles/second.
+
+    Given a window length `n` and a sample spacing `d`::
+
+      f = [0, 1, ...,   n/2-1,     -n/2, ..., -1] / (d*n)   if n is even
+      f = [0, 1, ..., (n-1)/2, -(n-1)/2, ..., -1] / (d*n)   if n is odd
+
+    Parameters
+    ----------
+    n : int
+        Window length.
+    d : scalar, optional
+        Sample spacing (inverse of the sampling rate). Defaults to 1.
+    device : str, optional
+        The device on which to place the created array. Default: ``None``.
+        For Array-API interoperability only, so must be ``"cpu"`` if passed.
+
+        .. versionadded:: 2.0.0
+
+    Returns
+    -------
+    f : ndarray
+        Array of length `n` containing the sample frequencies.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> signal = np.array([-2, 8, 6, 4, 1, 0, 3, 5], dtype=float)
+    >>> fourier = np.fft.fft(signal)
+    >>> n = signal.size
+    >>> timestep = 0.1
+    >>> freq = np.fft.fftfreq(n, d=timestep)
+    >>> freq
+    array([ 0.  ,  1.25,  2.5 , ..., -3.75, -2.5 , -1.25])
+
+    """
+    if not isinstance(n, integer_types):
+        raise ValueError("n should be an integer")
+    val = 1.0 / (n * d)
+    results = empty(n, int, device=device)
+    N = (n - 1) // 2 + 1
+    p1 = arange(0, N, dtype=int, device=device)
+    results[:N] = p1
+    p2 = arange(-(n // 2), 0, dtype=int, device=device)
+    results[N:] = p2
+    return results * val
+
+
+@set_module('numpy.fft')
+def rfftfreq(n, d=1.0, device=None):
+    """
+    Return the Discrete Fourier Transform sample frequencies
+    (for usage with rfft, irfft).
+
+    The returned float array `f` contains the frequency bin centers in cycles
+    per unit of the sample spacing (with zero at the start).  For instance, if
+    the sample spacing is in seconds, then the frequency unit is cycles/second.
+
+    Given a window length `n` and a sample spacing `d`::
+
+      f = [0, 1, ...,     n/2-1,     n/2] / (d*n)   if n is even
+      f = [0, 1, ..., (n-1)/2-1, (n-1)/2] / (d*n)   if n is odd
+
+    Unlike `fftfreq` (but like `scipy.fftpack.rfftfreq`)
+    the Nyquist frequency component is considered to be positive.
+
+    Parameters
+    ----------
+    n : int
+        Window length.
+    d : scalar, optional
+        Sample spacing (inverse of the sampling rate). Defaults to 1.
+    device : str, optional
+        The device on which to place the created array. Default: ``None``.
+        For Array-API interoperability only, so must be ``"cpu"`` if passed.
+
+        .. versionadded:: 2.0.0
+
+    Returns
+    -------
+    f : ndarray
+        Array of length ``n//2 + 1`` containing the sample frequencies.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> signal = np.array([-2, 8, 6, 4, 1, 0, 3, 5, -3, 4], dtype=float)
+    >>> fourier = np.fft.rfft(signal)
+    >>> n = signal.size
+    >>> sample_rate = 100
+    >>> freq = np.fft.fftfreq(n, d=1./sample_rate)
+    >>> freq
+    array([  0.,  10.,  20., ..., -30., -20., -10.])
+    >>> freq = np.fft.rfftfreq(n, d=1./sample_rate)
+    >>> freq
+    array([  0.,  10.,  20.,  30.,  40.,  50.])
+
+    """
+    if not isinstance(n, integer_types):
+        raise ValueError("n should be an integer")
+    val = 1.0 / (n * d)
+    N = n // 2 + 1
+    results = arange(0, N, dtype=int, device=device)
+    return results * val
diff --git a/.venv/lib/python3.12/site-packages/numpy/fft/_helper.pyi b/.venv/lib/python3.12/site-packages/numpy/fft/_helper.pyi
new file mode 100644
index 0000000..d06bda7
--- /dev/null
+++ b/.venv/lib/python3.12/site-packages/numpy/fft/_helper.pyi
@@ -0,0 +1,45 @@
+from typing import Any, Final, TypeVar, overload
+from typing import Literal as L
+
+from numpy import complexfloating, floating, generic, integer
+from numpy._typing import (
+    ArrayLike,
+    NDArray,
+    _ArrayLike,
+    _ArrayLikeComplex_co,
+    _ArrayLikeFloat_co,
+    _ShapeLike,
+)
+
+__all__ = ["fftfreq", "fftshift", "ifftshift", "rfftfreq"]
+
+_ScalarT = TypeVar("_ScalarT", bound=generic)
+
+###
+
+integer_types: Final[tuple[type[int], type[integer]]] = ...
+
+###
+
+@overload
+def fftshift(x: _ArrayLike[_ScalarT], axes: _ShapeLike | None = None) -> NDArray[_ScalarT]: ...
+@overload
+def fftshift(x: ArrayLike, axes: _ShapeLike | None = None) -> NDArray[Any]: ...
+
+#
+@overload
+def ifftshift(x: _ArrayLike[_ScalarT], axes: _ShapeLike | None = None) -> NDArray[_ScalarT]: ...
+@overload
+def ifftshift(x: ArrayLike, axes: _ShapeLike | None = None) -> NDArray[Any]: ...
+
+#
+@overload
+def fftfreq(n: int | integer, d: _ArrayLikeFloat_co = 1.0, device: L["cpu"] | None = None) -> NDArray[floating]: ...
+@overload
+def fftfreq(n: int | integer, d: _ArrayLikeComplex_co = 1.0, device: L["cpu"] | None = None) -> NDArray[complexfloating]: ...
+
+#
+@overload
+def rfftfreq(n: int | integer, d: _ArrayLikeFloat_co = 1.0, device: L["cpu"] | None = None) -> NDArray[floating]: ...
+@overload
+def rfftfreq(n: int | integer, d: _ArrayLikeComplex_co = 1.0, device: L["cpu"] | None = None) -> NDArray[complexfloating]: ...
diff --git a/.venv/lib/python3.12/site-packages/numpy/fft/_pocketfft.py b/.venv/lib/python3.12/site-packages/numpy/fft/_pocketfft.py
new file mode 100644
index 0000000..c7f2f6a
--- /dev/null
+++ b/.venv/lib/python3.12/site-packages/numpy/fft/_pocketfft.py
@@ -0,0 +1,1693 @@
+"""
+Discrete Fourier Transforms
+
+Routines in this module:
+
+fft(a, n=None, axis=-1, norm="backward")
+ifft(a, n=None, axis=-1, norm="backward")
+rfft(a, n=None, axis=-1, norm="backward")
+irfft(a, n=None, axis=-1, norm="backward")
+hfft(a, n=None, axis=-1, norm="backward")
+ihfft(a, n=None, axis=-1, norm="backward")
+fftn(a, s=None, axes=None, norm="backward")
+ifftn(a, s=None, axes=None, norm="backward")
+rfftn(a, s=None, axes=None, norm="backward")
+irfftn(a, s=None, axes=None, norm="backward")
+fft2(a, s=None, axes=(-2,-1), norm="backward")
+ifft2(a, s=None, axes=(-2, -1), norm="backward")
+rfft2(a, s=None, axes=(-2,-1), norm="backward")
+irfft2(a, s=None, axes=(-2, -1), norm="backward")
+
+i = inverse transform
+r = transform of purely real data
+h = Hermite transform
+n = n-dimensional transform
+2 = 2-dimensional transform
+(Note: 2D routines are just nD routines with different default
+behavior.)
+
+"""
+__all__ = ['fft', 'ifft', 'rfft', 'irfft', 'hfft', 'ihfft', 'rfftn',
+           'irfftn', 'rfft2', 'irfft2', 'fft2', 'ifft2', 'fftn', 'ifftn']
+
+import functools
+import warnings
+
+from numpy._core import (
+    asarray,
+    conjugate,
+    empty_like,
+    overrides,
+    reciprocal,
+    result_type,
+    sqrt,
+    take,
+)
+from numpy.lib.array_utils import normalize_axis_index
+
+from . import _pocketfft_umath as pfu
+
+array_function_dispatch = functools.partial(
+    overrides.array_function_dispatch, module='numpy.fft')
+
+
+# `inv_norm` is a float by which the result of the transform needs to be
+# divided. This replaces the original, more intuitive 'fct` parameter to avoid
+# divisions by zero (or alternatively additional checks) in the case of
+# zero-length axes during its computation.
+def _raw_fft(a, n, axis, is_real, is_forward, norm, out=None):
+    if n < 1:
+        raise ValueError(f"Invalid number of FFT data points ({n}) specified.")
+
+    # Calculate the normalization factor, passing in the array dtype to
+    # avoid precision loss in the possible sqrt or reciprocal.
+    if not is_forward:
+        norm = _swap_direction(norm)
+
+    real_dtype = result_type(a.real.dtype, 1.0)
+    if norm is None or norm == "backward":
+        fct = 1
+    elif norm == "ortho":
+        fct = reciprocal(sqrt(n, dtype=real_dtype))
+    elif norm == "forward":
+        fct = reciprocal(n, dtype=real_dtype)
+    else:
+        raise ValueError(f'Invalid norm value {norm}; should be "backward",'
+                         '"ortho" or "forward".')
+
+    n_out = n
+    if is_real:
+        if is_forward:
+            ufunc = pfu.rfft_n_even if n % 2 == 0 else pfu.rfft_n_odd
+            n_out = n // 2 + 1
+        else:
+            ufunc = pfu.irfft
+    else:
+        ufunc = pfu.fft if is_forward else pfu.ifft
+
+    axis = normalize_axis_index(axis, a.ndim)
+
+    if out is None:
+        if is_real and not is_forward:  # irfft, complex in, real output.
+            out_dtype = real_dtype
+        else:  # Others, complex output.
+            out_dtype = result_type(a.dtype, 1j)
+        out = empty_like(a, shape=a.shape[:axis] + (n_out,) + a.shape[axis + 1:],
+                         dtype=out_dtype)
+    elif ((shape := getattr(out, "shape", None)) is not None
+          and (len(shape) != a.ndim or shape[axis] != n_out)):
+        raise ValueError("output array has wrong shape.")
+
+    return ufunc(a, fct, axes=[(axis,), (), (axis,)], out=out)
+
+
+_SWAP_DIRECTION_MAP = {"backward": "forward", None: "forward",
+                       "ortho": "ortho", "forward": "backward"}
+
+
+def _swap_direction(norm):
+    try:
+        return _SWAP_DIRECTION_MAP[norm]
+    except KeyError:
+        raise ValueError(f'Invalid norm value {norm}; should be "backward", '
+                         '"ortho" or "forward".') from None
+
+
+def _fft_dispatcher(a, n=None, axis=None, norm=None, out=None):
+    return (a, out)
+
+
+@array_function_dispatch(_fft_dispatcher)
+def fft(a, n=None, axis=-1, norm=None, out=None):
+    """
+    Compute the one-dimensional discrete Fourier Transform.
+
+    This function computes the one-dimensional *n*-point discrete Fourier
+    Transform (DFT) with the efficient Fast Fourier Transform (FFT)
+    algorithm [CT].
+
+    Parameters
+    ----------
+    a : array_like
+        Input array, can be complex.
+    n : int, optional
+        Length of the transformed axis of the output.
+        If `n` is smaller than the length of the input, the input is cropped.
+        If it is larger, the input is padded with zeros.  If `n` is not given,
+        the length of the input along the axis specified by `axis` is used.
+    axis : int, optional
+        Axis over which to compute the FFT.  If not given, the last axis is
+        used.
+    norm : {"backward", "ortho", "forward"}, optional
+        Normalization mode (see `numpy.fft`). Default is "backward".
+        Indicates which direction of the forward/backward pair of transforms
+        is scaled and with what normalization factor.
+
+        .. versionadded:: 1.20.0
+
+            The "backward", "forward" values were added.
+    out : complex ndarray, optional
+        If provided, the result will be placed in this array. It should be
+        of the appropriate shape and dtype.
+
+        .. versionadded:: 2.0.0
+
+    Returns
+    -------
+    out : complex ndarray
+        The truncated or zero-padded input, transformed along the axis
+        indicated by `axis`, or the last one if `axis` is not specified.
+
+    Raises
+    ------
+    IndexError
+        If `axis` is not a valid axis of `a`.
+
+    See Also
+    --------
+    numpy.fft : for definition of the DFT and conventions used.
+    ifft : The inverse of `fft`.
+    fft2 : The two-dimensional FFT.
+    fftn : The *n*-dimensional FFT.
+    rfftn : The *n*-dimensional FFT of real input.
+    fftfreq : Frequency bins for given FFT parameters.
+
+    Notes
+    -----
+    FFT (Fast Fourier Transform) refers to a way the discrete Fourier
+    Transform (DFT) can be calculated efficiently, by using symmetries in the
+    calculated terms.  The symmetry is highest when `n` is a power of 2, and
+    the transform is therefore most efficient for these sizes.
+
+    The DFT is defined, with the conventions used in this implementation, in
+    the documentation for the `numpy.fft` module.
+
+    References
+    ----------
+    .. [CT] Cooley, James W., and John W. Tukey, 1965, "An algorithm for the
+            machine calculation of complex Fourier series," *Math. Comput.*
+            19: 297-301.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> np.fft.fft(np.exp(2j * np.pi * np.arange(8) / 8))
+    array([-2.33486982e-16+1.14423775e-17j,  8.00000000e+00-1.25557246e-15j,
+            2.33486982e-16+2.33486982e-16j,  0.00000000e+00+1.22464680e-16j,
+           -1.14423775e-17+2.33486982e-16j,  0.00000000e+00+5.20784380e-16j,
+            1.14423775e-17+1.14423775e-17j,  0.00000000e+00+1.22464680e-16j])
+
+    In this example, real input has an FFT which is Hermitian, i.e., symmetric
+    in the real part and anti-symmetric in the imaginary part, as described in
+    the `numpy.fft` documentation:
+
+    >>> import matplotlib.pyplot as plt
+    >>> t = np.arange(256)
+    >>> sp = np.fft.fft(np.sin(t))
+    >>> freq = np.fft.fftfreq(t.shape[-1])
+    >>> _ = plt.plot(freq, sp.real, freq, sp.imag)
+    >>> plt.show()
+
+    """
+    a = asarray(a)
+    if n is None:
+        n = a.shape[axis]
+    output = _raw_fft(a, n, axis, False, True, norm, out)
+    return output
+
+
+@array_function_dispatch(_fft_dispatcher)
+def ifft(a, n=None, axis=-1, norm=None, out=None):
+    """
+    Compute the one-dimensional inverse discrete Fourier Transform.
+
+    This function computes the inverse of the one-dimensional *n*-point
+    discrete Fourier transform computed by `fft`.  In other words,
+    ``ifft(fft(a)) == a`` to within numerical accuracy.
+    For a general description of the algorithm and definitions,
+    see `numpy.fft`.
+
+    The input should be ordered in the same way as is returned by `fft`,
+    i.e.,
+
+    * ``a[0]`` should contain the zero frequency term,
+    * ``a[1:n//2]`` should contain the positive-frequency terms,
+    * ``a[n//2 + 1:]`` should contain the negative-frequency terms, in
+      increasing order starting from the most negative frequency.
+
+    For an even number of input points, ``A[n//2]`` represents the sum of
+    the values at the positive and negative Nyquist frequencies, as the two
+    are aliased together. See `numpy.fft` for details.
+
+    Parameters
+    ----------
+    a : array_like
+        Input array, can be complex.
+    n : int, optional
+        Length of the transformed axis of the output.
+        If `n` is smaller than the length of the input, the input is cropped.
+        If it is larger, the input is padded with zeros.  If `n` is not given,
+        the length of the input along the axis specified by `axis` is used.
+        See notes about padding issues.
+    axis : int, optional
+        Axis over which to compute the inverse DFT.  If not given, the last
+        axis is used.
+    norm : {"backward", "ortho", "forward"}, optional
+        Normalization mode (see `numpy.fft`). Default is "backward".
+        Indicates which direction of the forward/backward pair of transforms
+        is scaled and with what normalization factor.
+
+        .. versionadded:: 1.20.0
+
+            The "backward", "forward" values were added.
+
+    out : complex ndarray, optional
+        If provided, the result will be placed in this array. It should be
+        of the appropriate shape and dtype.
+
+        .. versionadded:: 2.0.0
+
+    Returns
+    -------
+    out : complex ndarray
+        The truncated or zero-padded input, transformed along the axis
+        indicated by `axis`, or the last one if `axis` is not specified.
+
+    Raises
+    ------
+    IndexError
+        If `axis` is not a valid axis of `a`.
+
+    See Also
+    --------
+    numpy.fft : An introduction, with definitions and general explanations.
+    fft : The one-dimensional (forward) FFT, of which `ifft` is the inverse
+    ifft2 : The two-dimensional inverse FFT.
+    ifftn : The n-dimensional inverse FFT.
+
+    Notes
+    -----
+    If the input parameter `n` is larger than the size of the input, the input
+    is padded by appending zeros at the end.  Even though this is the common
+    approach, it might lead to surprising results.  If a different padding is
+    desired, it must be performed before calling `ifft`.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> np.fft.ifft([0, 4, 0, 0])
+    array([ 1.+0.j,  0.+1.j, -1.+0.j,  0.-1.j]) # may vary
+
+    Create and plot a band-limited signal with random phases:
+
+    >>> import matplotlib.pyplot as plt
+    >>> t = np.arange(400)
+    >>> n = np.zeros((400,), dtype=complex)
+    >>> n[40:60] = np.exp(1j*np.random.uniform(0, 2*np.pi, (20,)))
+    >>> s = np.fft.ifft(n)
+    >>> plt.plot(t, s.real, label='real')
+    [<matplotlib.lines.Line2D object at ...>]
+    >>> plt.plot(t, s.imag, '--', label='imaginary')
+    [<matplotlib.lines.Line2D object at ...>]
+    >>> plt.legend()
+    <matplotlib.legend.Legend object at ...>
+    >>> plt.show()
+
+    """
+    a = asarray(a)
+    if n is None:
+        n = a.shape[axis]
+    output = _raw_fft(a, n, axis, False, False, norm, out=out)
+    return output
+
+
+@array_function_dispatch(_fft_dispatcher)
+def rfft(a, n=None, axis=-1, norm=None, out=None):
+    """
+    Compute the one-dimensional discrete Fourier Transform for real input.
+
+    This function computes the one-dimensional *n*-point discrete Fourier
+    Transform (DFT) of a real-valued array by means of an efficient algorithm
+    called the Fast Fourier Transform (FFT).
+
+    Parameters
+    ----------
+    a : array_like
+        Input array
+    n : int, optional
+        Number of points along transformation axis in the input to use.
+        If `n` is smaller than the length of the input, the input is cropped.
+        If it is larger, the input is padded with zeros. If `n` is not given,
+        the length of the input along the axis specified by `axis` is used.
+    axis : int, optional
+        Axis over which to compute the FFT. If not given, the last axis is
+        used.
+    norm : {"backward", "ortho", "forward"}, optional
+        Normalization mode (see `numpy.fft`). Default is "backward".
+        Indicates which direction of the forward/backward pair of transforms
+        is scaled and with what normalization factor.
+
+        .. versionadded:: 1.20.0
+
+            The "backward", "forward" values were added.
+
+    out : complex ndarray, optional
+        If provided, the result will be placed in this array. It should be
+        of the appropriate shape and dtype.
+
+        .. versionadded:: 2.0.0
+
+    Returns
+    -------
+    out : complex ndarray
+        The truncated or zero-padded input, transformed along the axis
+        indicated by `axis`, or the last one if `axis` is not specified.
+        If `n` is even, the length of the transformed axis is ``(n/2)+1``.
+        If `n` is odd, the length is ``(n+1)/2``.
+
+    Raises
+    ------
+    IndexError
+        If `axis` is not a valid axis of `a`.
+
+    See Also
+    --------
+    numpy.fft : For definition of the DFT and conventions used.
+    irfft : The inverse of `rfft`.
+    fft : The one-dimensional FFT of general (complex) input.
+    fftn : The *n*-dimensional FFT.
+    rfftn : The *n*-dimensional FFT of real input.
+
+    Notes
+    -----
+    When the DFT is computed for purely real input, the output is
+    Hermitian-symmetric, i.e. the negative frequency terms are just the complex
+    conjugates of the corresponding positive-frequency terms, and the
+    negative-frequency terms are therefore redundant.  This function does not
+    compute the negative frequency terms, and the length of the transformed
+    axis of the output is therefore ``n//2 + 1``.
+
+    When ``A = rfft(a)`` and fs is the sampling frequency, ``A[0]`` contains
+    the zero-frequency term 0*fs, which is real due to Hermitian symmetry.
+
+    If `n` is even, ``A[-1]`` contains the term representing both positive
+    and negative Nyquist frequency (+fs/2 and -fs/2), and must also be purely
+    real. If `n` is odd, there is no term at fs/2; ``A[-1]`` contains
+    the largest positive frequency (fs/2*(n-1)/n), and is complex in the
+    general case.
+
+    If the input `a` contains an imaginary part, it is silently discarded.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> np.fft.fft([0, 1, 0, 0])
+    array([ 1.+0.j,  0.-1.j, -1.+0.j,  0.+1.j]) # may vary
+    >>> np.fft.rfft([0, 1, 0, 0])
+    array([ 1.+0.j,  0.-1.j, -1.+0.j]) # may vary
+
+    Notice how the final element of the `fft` output is the complex conjugate
+    of the second element, for real input. For `rfft`, this symmetry is
+    exploited to compute only the non-negative frequency terms.
+
+    """
+    a = asarray(a)
+    if n is None:
+        n = a.shape[axis]
+    output = _raw_fft(a, n, axis, True, True, norm, out=out)
+    return output
+
+
+@array_function_dispatch(_fft_dispatcher)
+def irfft(a, n=None, axis=-1, norm=None, out=None):
+    """
+    Computes the inverse of `rfft`.
+
+    This function computes the inverse of the one-dimensional *n*-point
+    discrete Fourier Transform of real input computed by `rfft`.
+    In other words, ``irfft(rfft(a), len(a)) == a`` to within numerical
+    accuracy. (See Notes below for why ``len(a)`` is necessary here.)
+
+    The input is expected to be in the form returned by `rfft`, i.e. the
+    real zero-frequency term followed by the complex positive frequency terms
+    in order of increasing frequency.  Since the discrete Fourier Transform of
+    real input is Hermitian-symmetric, the negative frequency terms are taken
+    to be the complex conjugates of the corresponding positive frequency terms.
+
+    Parameters
+    ----------
+    a : array_like
+        The input array.
+    n : int, optional
+        Length of the transformed axis of the output.
+        For `n` output points, ``n//2+1`` input points are necessary.  If the
+        input is longer than this, it is cropped.  If it is shorter than this,
+        it is padded with zeros.  If `n` is not given, it is taken to be
+        ``2*(m-1)`` where ``m`` is the length of the input along the axis
+        specified by `axis`.
+    axis : int, optional
+        Axis over which to compute the inverse FFT. If not given, the last
+        axis is used.
+    norm : {"backward", "ortho", "forward"}, optional
+        Normalization mode (see `numpy.fft`). Default is "backward".
+        Indicates which direction of the forward/backward pair of transforms
+        is scaled and with what normalization factor.
+
+        .. versionadded:: 1.20.0
+
+            The "backward", "forward" values were added.
+
+    out : ndarray, optional
+        If provided, the result will be placed in this array. It should be
+        of the appropriate shape and dtype.
+
+        .. versionadded:: 2.0.0
+
+    Returns
+    -------
+    out : ndarray
+        The truncated or zero-padded input, transformed along the axis
+        indicated by `axis`, or the last one if `axis` is not specified.
+        The length of the transformed axis is `n`, or, if `n` is not given,
+        ``2*(m-1)`` where ``m`` is the length of the transformed axis of the
+        input. To get an odd number of output points, `n` must be specified.
+
+    Raises
+    ------
+    IndexError
+        If `axis` is not a valid axis of `a`.
+
+    See Also
+    --------
+    numpy.fft : For definition of the DFT and conventions used.
+    rfft : The one-dimensional FFT of real input, of which `irfft` is inverse.
+    fft : The one-dimensional FFT.
+    irfft2 : The inverse of the two-dimensional FFT of real input.
+    irfftn : The inverse of the *n*-dimensional FFT of real input.
+
+    Notes
+    -----
+    Returns the real valued `n`-point inverse discrete Fourier transform
+    of `a`, where `a` contains the non-negative frequency terms of a
+    Hermitian-symmetric sequence. `n` is the length of the result, not the
+    input.
+
+    If you specify an `n` such that `a` must be zero-padded or truncated, the
+    extra/removed values will be added/removed at high frequencies. One can
+    thus resample a series to `m` points via Fourier interpolation by:
+    ``a_resamp = irfft(rfft(a), m)``.
+
+    The correct interpretation of the hermitian input depends on the length of
+    the original data, as given by `n`. This is because each input shape could
+    correspond to either an odd or even length signal. By default, `irfft`
+    assumes an even output length which puts the last entry at the Nyquist
+    frequency; aliasing with its symmetric counterpart. By Hermitian symmetry,
+    the value is thus treated as purely real. To avoid losing information, the
+    correct length of the real input **must** be given.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> np.fft.ifft([1, -1j, -1, 1j])
+    array([0.+0.j,  1.+0.j,  0.+0.j,  0.+0.j]) # may vary
+    >>> np.fft.irfft([1, -1j, -1])
+    array([0.,  1.,  0.,  0.])
+
+    Notice how the last term in the input to the ordinary `ifft` is the
+    complex conjugate of the second term, and the output has zero imaginary
+    part everywhere.  When calling `irfft`, the negative frequencies are not
+    specified, and the output array is purely real.
+
+    """
+    a = asarray(a)
+    if n is None:
+        n = (a.shape[axis] - 1) * 2
+    output = _raw_fft(a, n, axis, True, False, norm, out=out)
+    return output
+
+
+@array_function_dispatch(_fft_dispatcher)
+def hfft(a, n=None, axis=-1, norm=None, out=None):
+    """
+    Compute the FFT of a signal that has Hermitian symmetry, i.e., a real
+    spectrum.
+
+    Parameters
+    ----------
+    a : array_like
+        The input array.
+    n : int, optional
+        Length of the transformed axis of the output. For `n` output
+        points, ``n//2 + 1`` input points are necessary.  If the input is
+        longer than this, it is cropped.  If it is shorter than this, it is
+        padded with zeros.  If `n` is not given, it is taken to be ``2*(m-1)``
+        where ``m`` is the length of the input along the axis specified by
+        `axis`.
+    axis : int, optional
+        Axis over which to compute the FFT. If not given, the last
+        axis is used.
+    norm : {"backward", "ortho", "forward"}, optional
+        Normalization mode (see `numpy.fft`). Default is "backward".
+        Indicates which direction of the forward/backward pair of transforms
+        is scaled and with what normalization factor.
+
+        .. versionadded:: 1.20.0
+
+            The "backward", "forward" values were added.
+
+    out : ndarray, optional
+        If provided, the result will be placed in this array. It should be
+        of the appropriate shape and dtype.
+
+        .. versionadded:: 2.0.0
+
+    Returns
+    -------
+    out : ndarray
+        The truncated or zero-padded input, transformed along the axis
+        indicated by `axis`, or the last one if `axis` is not specified.
+        The length of the transformed axis is `n`, or, if `n` is not given,
+        ``2*m - 2`` where ``m`` is the length of the transformed axis of
+        the input. To get an odd number of output points, `n` must be
+        specified, for instance as ``2*m - 1`` in the typical case,
+
+    Raises
+    ------
+    IndexError
+        If `axis` is not a valid axis of `a`.
+
+    See also
+    --------
+    rfft : Compute the one-dimensional FFT for real input.
+    ihfft : The inverse of `hfft`.
+
+    Notes
+    -----
+    `hfft`/`ihfft` are a pair analogous to `rfft`/`irfft`, but for the
+    opposite case: here the signal has Hermitian symmetry in the time
+    domain and is real in the frequency domain. So here it's `hfft` for
+    which you must supply the length of the result if it is to be odd.
+
+    * even: ``ihfft(hfft(a, 2*len(a) - 2)) == a``, within roundoff error,
+    * odd: ``ihfft(hfft(a, 2*len(a) - 1)) == a``, within roundoff error.
+
+    The correct interpretation of the hermitian input depends on the length of
+    the original data, as given by `n`. This is because each input shape could
+    correspond to either an odd or even length signal. By default, `hfft`
+    assumes an even output length which puts the last entry at the Nyquist
+    frequency; aliasing with its symmetric counterpart. By Hermitian symmetry,
+    the value is thus treated as purely real. To avoid losing information, the
+    shape of the full signal **must** be given.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> signal = np.array([1, 2, 3, 4, 3, 2])
+    >>> np.fft.fft(signal)
+    array([15.+0.j,  -4.+0.j,   0.+0.j,  -1.-0.j,   0.+0.j,  -4.+0.j]) # may vary
+    >>> np.fft.hfft(signal[:4]) # Input first half of signal
+    array([15.,  -4.,   0.,  -1.,   0.,  -4.])
+    >>> np.fft.hfft(signal, 6)  # Input entire signal and truncate
+    array([15.,  -4.,   0.,  -1.,   0.,  -4.])
+
+
+    >>> signal = np.array([[1, 1.j], [-1.j, 2]])
+    >>> np.conj(signal.T) - signal   # check Hermitian symmetry
+    array([[ 0.-0.j,  -0.+0.j], # may vary
+           [ 0.+0.j,  0.-0.j]])
+    >>> freq_spectrum = np.fft.hfft(signal)
+    >>> freq_spectrum
+    array([[ 1.,  1.],
+           [ 2., -2.]])
+
+    """
+    a = asarray(a)
+    if n is None:
+        n = (a.shape[axis] - 1) * 2
+    new_norm = _swap_direction(norm)
+    output = irfft(conjugate(a), n, axis, norm=new_norm, out=None)
+    return output
+
+
+@array_function_dispatch(_fft_dispatcher)
+def ihfft(a, n=None, axis=-1, norm=None, out=None):
+    """
+    Compute the inverse FFT of a signal that has Hermitian symmetry.
+
+    Parameters
+    ----------
+    a : array_like
+        Input array.
+    n : int, optional
+        Length of the inverse FFT, the number of points along
+        transformation axis in the input to use.  If `n` is smaller than
+        the length of the input, the input is cropped.  If it is larger,
+        the input is padded with zeros. If `n` is not given, the length of
+        the input along the axis specified by `axis` is used.
+    axis : int, optional
+        Axis over which to compute the inverse FFT. If not given, the last
+        axis is used.
+    norm : {"backward", "ortho", "forward"}, optional
+        Normalization mode (see `numpy.fft`). Default is "backward".
+        Indicates which direction of the forward/backward pair of transforms
+        is scaled and with what normalization factor.
+
+        .. versionadded:: 1.20.0
+
+            The "backward", "forward" values were added.
+
+    out : complex ndarray, optional
+        If provided, the result will be placed in this array. It should be
+        of the appropriate shape and dtype.
+
+        .. versionadded:: 2.0.0
+
+    Returns
+    -------
+    out : complex ndarray
+        The truncated or zero-padded input, transformed along the axis
+        indicated by `axis`, or the last one if `axis` is not specified.
+        The length of the transformed axis is ``n//2 + 1``.
+
+    See also
+    --------
+    hfft, irfft
+
+    Notes
+    -----
+    `hfft`/`ihfft` are a pair analogous to `rfft`/`irfft`, but for the
+    opposite case: here the signal has Hermitian symmetry in the time
+    domain and is real in the frequency domain. So here it's `hfft` for
+    which you must supply the length of the result if it is to be odd:
+
+    * even: ``ihfft(hfft(a, 2*len(a) - 2)) == a``, within roundoff error,
+    * odd: ``ihfft(hfft(a, 2*len(a) - 1)) == a``, within roundoff error.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> spectrum = np.array([ 15, -4, 0, -1, 0, -4])
+    >>> np.fft.ifft(spectrum)
+    array([1.+0.j,  2.+0.j,  3.+0.j,  4.+0.j,  3.+0.j,  2.+0.j]) # may vary
+    >>> np.fft.ihfft(spectrum)
+    array([ 1.-0.j,  2.-0.j,  3.-0.j,  4.-0.j]) # may vary
+
+    """
+    a = asarray(a)
+    if n is None:
+        n = a.shape[axis]
+    new_norm = _swap_direction(norm)
+    out = rfft(a, n, axis, norm=new_norm, out=out)
+    return conjugate(out, out=out)
+
+
+def _cook_nd_args(a, s=None, axes=None, invreal=0):
+    if s is None:
+        shapeless = True
+        if axes is None:
+            s = list(a.shape)
+        else:
+            s = take(a.shape, axes)
+    else:
+        shapeless = False
+    s = list(s)
+    if axes is None:
+        if not shapeless:
+            msg = ("`axes` should not be `None` if `s` is not `None` "
+                   "(Deprecated in NumPy 2.0). In a future version of NumPy, "
+                   "this will raise an error and `s[i]` will correspond to "
+                   "the size along the transformed axis specified by "
+                   "`axes[i]`. To retain current behaviour, pass a sequence "
+                   "[0, ..., k-1] to `axes` for an array of dimension k.")
+            warnings.warn(msg, DeprecationWarning, stacklevel=3)
+        axes = list(range(-len(s), 0))
+    if len(s) != len(axes):
+        raise ValueError("Shape and axes have different lengths.")
+    if invreal and shapeless:
+        s[-1] = (a.shape[axes[-1]] - 1) * 2
+    if None in s:
+        msg = ("Passing an array containing `None` values to `s` is "
+               "deprecated in NumPy 2.0 and will raise an error in "
+               "a future version of NumPy. To use the default behaviour "
+               "of the corresponding 1-D transform, pass the value matching "
+               "the default for its `n` parameter. To use the default "
+               "behaviour for every axis, the `s` argument can be omitted.")
+        warnings.warn(msg, DeprecationWarning, stacklevel=3)
+    # use the whole input array along axis `i` if `s[i] == -1`
+    s = [a.shape[_a] if _s == -1 else _s for _s, _a in zip(s, axes)]
+    return s, axes
+
+
+def _raw_fftnd(a, s=None, axes=None, function=fft, norm=None, out=None):
+    a = asarray(a)
+    s, axes = _cook_nd_args(a, s, axes)
+    itl = list(range(len(axes)))
+    itl.reverse()
+    for ii in itl:
+        a = function(a, n=s[ii], axis=axes[ii], norm=norm, out=out)
+    return a
+
+
+def _fftn_dispatcher(a, s=None, axes=None, norm=None, out=None):
+    return (a, out)
+
+
+@array_function_dispatch(_fftn_dispatcher)
+def fftn(a, s=None, axes=None, norm=None, out=None):
+    """
+    Compute the N-dimensional discrete Fourier Transform.
+
+    This function computes the *N*-dimensional discrete Fourier Transform over
+    any number of axes in an *M*-dimensional array by means of the Fast Fourier
+    Transform (FFT).
+
+    Parameters
+    ----------
+    a : array_like
+        Input array, can be complex.
+    s : sequence of ints, optional
+        Shape (length of each transformed axis) of the output
+        (``s[0]`` refers to axis 0, ``s[1]`` to axis 1, etc.).
+        This corresponds to ``n`` for ``fft(x, n)``.
+        Along any axis, if the given shape is smaller than that of the input,
+        the input is cropped. If it is larger, the input is padded with zeros.
+
+        .. versionchanged:: 2.0
+
+            If it is ``-1``, the whole input is used (no padding/trimming).
+
+        If `s` is not given, the shape of the input along the axes specified
+        by `axes` is used.
+
+        .. deprecated:: 2.0
+
+            If `s` is not ``None``, `axes` must not be ``None`` either.
+
+        .. deprecated:: 2.0
+
+            `s` must contain only ``int`` s, not ``None`` values. ``None``
+            values currently mean that the default value for ``n`` is used
+            in the corresponding 1-D transform, but this behaviour is
+            deprecated.
+
+    axes : sequence of ints, optional
+        Axes over which to compute the FFT.  If not given, the last ``len(s)``
+        axes are used, or all axes if `s` is also not specified.
+        Repeated indices in `axes` means that the transform over that axis is
+        performed multiple times.
+
+        .. deprecated:: 2.0
+
+            If `s` is specified, the corresponding `axes` to be transformed
+            must be explicitly specified too.
+
+    norm : {"backward", "ortho", "forward"}, optional
+        Normalization mode (see `numpy.fft`). Default is "backward".
+        Indicates which direction of the forward/backward pair of transforms
+        is scaled and with what normalization factor.
+
+        .. versionadded:: 1.20.0
+
+            The "backward", "forward" values were added.
+
+    out : complex ndarray, optional
+        If provided, the result will be placed in this array. It should be
+        of the appropriate shape and dtype for all axes (and hence is
+        incompatible with passing in all but the trivial ``s``).
+
+        .. versionadded:: 2.0.0
+
+    Returns
+    -------
+    out : complex ndarray
+        The truncated or zero-padded input, transformed along the axes
+        indicated by `axes`, or by a combination of `s` and `a`,
+        as explained in the parameters section above.
+
+    Raises
+    ------
+    ValueError
+        If `s` and `axes` have different length.
+    IndexError
+        If an element of `axes` is larger than than the number of axes of `a`.
+
+    See Also
+    --------
+    numpy.fft : Overall view of discrete Fourier transforms, with definitions
+        and conventions used.
+    ifftn : The inverse of `fftn`, the inverse *n*-dimensional FFT.
+    fft : The one-dimensional FFT, with definitions and conventions used.
+    rfftn : The *n*-dimensional FFT of real input.
+    fft2 : The two-dimensional FFT.
+    fftshift : Shifts zero-frequency terms to centre of array
+
+    Notes
+    -----
+    The output, analogously to `fft`, contains the term for zero frequency in
+    the low-order corner of all axes, the positive frequency terms in the
+    first half of all axes, the term for the Nyquist frequency in the middle
+    of all axes and the negative frequency terms in the second half of all
+    axes, in order of decreasingly negative frequency.
+
+    See `numpy.fft` for details, definitions and conventions used.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> a = np.mgrid[:3, :3, :3][0]
+    >>> np.fft.fftn(a, axes=(1, 2))
+    array([[[ 0.+0.j,   0.+0.j,   0.+0.j], # may vary
+            [ 0.+0.j,   0.+0.j,   0.+0.j],
+            [ 0.+0.j,   0.+0.j,   0.+0.j]],
+           [[ 9.+0.j,   0.+0.j,   0.+0.j],
+            [ 0.+0.j,   0.+0.j,   0.+0.j],
+            [ 0.+0.j,   0.+0.j,   0.+0.j]],
+           [[18.+0.j,   0.+0.j,   0.+0.j],
+            [ 0.+0.j,   0.+0.j,   0.+0.j],
+            [ 0.+0.j,   0.+0.j,   0.+0.j]]])
+    >>> np.fft.fftn(a, (2, 2), axes=(0, 1))
+    array([[[ 2.+0.j,  2.+0.j,  2.+0.j], # may vary
+            [ 0.+0.j,  0.+0.j,  0.+0.j]],
+           [[-2.+0.j, -2.+0.j, -2.+0.j],
+            [ 0.+0.j,  0.+0.j,  0.+0.j]]])
+
+    >>> import matplotlib.pyplot as plt
+    >>> [X, Y] = np.meshgrid(2 * np.pi * np.arange(200) / 12,
+    ...                      2 * np.pi * np.arange(200) / 34)
+    >>> S = np.sin(X) + np.cos(Y) + np.random.uniform(0, 1, X.shape)
+    >>> FS = np.fft.fftn(S)
+    >>> plt.imshow(np.log(np.abs(np.fft.fftshift(FS))**2))
+    <matplotlib.image.AxesImage object at 0x...>
+    >>> plt.show()
+
+    """
+    return _raw_fftnd(a, s, axes, fft, norm, out=out)
+
+
+@array_function_dispatch(_fftn_dispatcher)
+def ifftn(a, s=None, axes=None, norm=None, out=None):
+    """
+    Compute the N-dimensional inverse discrete Fourier Transform.
+
+    This function computes the inverse of the N-dimensional discrete
+    Fourier Transform over any number of axes in an M-dimensional array by
+    means of the Fast Fourier Transform (FFT).  In other words,
+    ``ifftn(fftn(a)) == a`` to within numerical accuracy.
+    For a description of the definitions and conventions used, see `numpy.fft`.
+
+    The input, analogously to `ifft`, should be ordered in the same way as is
+    returned by `fftn`, i.e. it should have the term for zero frequency
+    in all axes in the low-order corner, the positive frequency terms in the
+    first half of all axes, the term for the Nyquist frequency in the middle
+    of all axes and the negative frequency terms in the second half of all
+    axes, in order of decreasingly negative frequency.
+
+    Parameters
+    ----------
+    a : array_like
+        Input array, can be complex.
+    s : sequence of ints, optional
+        Shape (length of each transformed axis) of the output
+        (``s[0]`` refers to axis 0, ``s[1]`` to axis 1, etc.).
+        This corresponds to ``n`` for ``ifft(x, n)``.
+        Along any axis, if the given shape is smaller than that of the input,
+        the input is cropped. If it is larger, the input is padded with zeros.
+
+        .. versionchanged:: 2.0
+
+            If it is ``-1``, the whole input is used (no padding/trimming).
+
+        If `s` is not given, the shape of the input along the axes specified
+        by `axes` is used. See notes for issue on `ifft` zero padding.
+
+        .. deprecated:: 2.0
+
+            If `s` is not ``None``, `axes` must not be ``None`` either.
+
+        .. deprecated:: 2.0
+
+            `s` must contain only ``int`` s, not ``None`` values. ``None``
+            values currently mean that the default value for ``n`` is used
+            in the corresponding 1-D transform, but this behaviour is
+            deprecated.
+
+    axes : sequence of ints, optional
+        Axes over which to compute the IFFT.  If not given, the last ``len(s)``
+        axes are used, or all axes if `s` is also not specified.
+        Repeated indices in `axes` means that the inverse transform over that
+        axis is performed multiple times.
+
+        .. deprecated:: 2.0
+
+            If `s` is specified, the corresponding `axes` to be transformed
+            must be explicitly specified too.
+
+    norm : {"backward", "ortho", "forward"}, optional
+        Normalization mode (see `numpy.fft`). Default is "backward".
+        Indicates which direction of the forward/backward pair of transforms
+        is scaled and with what normalization factor.
+
+        .. versionadded:: 1.20.0
+
+            The "backward", "forward" values were added.
+
+    out : complex ndarray, optional
+        If provided, the result will be placed in this array. It should be
+        of the appropriate shape and dtype for all axes (and hence is
+        incompatible with passing in all but the trivial ``s``).
+
+        .. versionadded:: 2.0.0
+
+    Returns
+    -------
+    out : complex ndarray
+        The truncated or zero-padded input, transformed along the axes
+        indicated by `axes`, or by a combination of `s` or `a`,
+        as explained in the parameters section above.
+
+    Raises
+    ------
+    ValueError
+        If `s` and `axes` have different length.
+    IndexError
+        If an element of `axes` is larger than than the number of axes of `a`.
+
+    See Also
+    --------
+    numpy.fft : Overall view of discrete Fourier transforms, with definitions
+         and conventions used.
+    fftn : The forward *n*-dimensional FFT, of which `ifftn` is the inverse.
+    ifft : The one-dimensional inverse FFT.
+    ifft2 : The two-dimensional inverse FFT.
+    ifftshift : Undoes `fftshift`, shifts zero-frequency terms to beginning
+        of array.
+
+    Notes
+    -----
+    See `numpy.fft` for definitions and conventions used.
+
+    Zero-padding, analogously with `ifft`, is performed by appending zeros to
+    the input along the specified dimension.  Although this is the common
+    approach, it might lead to surprising results.  If another form of zero
+    padding is desired, it must be performed before `ifftn` is called.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> a = np.eye(4)
+    >>> np.fft.ifftn(np.fft.fftn(a, axes=(0,)), axes=(1,))
+    array([[1.+0.j,  0.+0.j,  0.+0.j,  0.+0.j], # may vary
+           [0.+0.j,  1.+0.j,  0.+0.j,  0.+0.j],
+           [0.+0.j,  0.+0.j,  1.+0.j,  0.+0.j],
+           [0.+0.j,  0.+0.j,  0.+0.j,  1.+0.j]])
+
+
+    Create and plot an image with band-limited frequency content:
+
+    >>> import matplotlib.pyplot as plt
+    >>> n = np.zeros((200,200), dtype=complex)
+    >>> n[60:80, 20:40] = np.exp(1j*np.random.uniform(0, 2*np.pi, (20, 20)))
+    >>> im = np.fft.ifftn(n).real
+    >>> plt.imshow(im)
+    <matplotlib.image.AxesImage object at 0x...>
+    >>> plt.show()
+
+    """
+    return _raw_fftnd(a, s, axes, ifft, norm, out=out)
+
+
+@array_function_dispatch(_fftn_dispatcher)
+def fft2(a, s=None, axes=(-2, -1), norm=None, out=None):
+    """
+    Compute the 2-dimensional discrete Fourier Transform.
+
+    This function computes the *n*-dimensional discrete Fourier Transform
+    over any axes in an *M*-dimensional array by means of the
+    Fast Fourier Transform (FFT).  By default, the transform is computed over
+    the last two axes of the input array, i.e., a 2-dimensional FFT.
+
+    Parameters
+    ----------
+    a : array_like
+        Input array, can be complex
+    s : sequence of ints, optional
+        Shape (length of each transformed axis) of the output
+        (``s[0]`` refers to axis 0, ``s[1]`` to axis 1, etc.).
+        This corresponds to ``n`` for ``fft(x, n)``.
+        Along each axis, if the given shape is smaller than that of the input,
+        the input is cropped. If it is larger, the input is padded with zeros.
+
+        .. versionchanged:: 2.0
+
+            If it is ``-1``, the whole input is used (no padding/trimming).
+
+        If `s` is not given, the shape of the input along the axes specified
+        by `axes` is used.
+
+        .. deprecated:: 2.0
+
+            If `s` is not ``None``, `axes` must not be ``None`` either.
+
+        .. deprecated:: 2.0
+
+            `s` must contain only ``int`` s, not ``None`` values. ``None``
+            values currently mean that the default value for ``n`` is used
+            in the corresponding 1-D transform, but this behaviour is
+            deprecated.
+
+    axes : sequence of ints, optional
+        Axes over which to compute the FFT.  If not given, the last two
+        axes are used.  A repeated index in `axes` means the transform over
+        that axis is performed multiple times.  A one-element sequence means
+        that a one-dimensional FFT is performed. Default: ``(-2, -1)``.
+
+        .. deprecated:: 2.0
+
+            If `s` is specified, the corresponding `axes` to be transformed
+            must not be ``None``.
+
+    norm : {"backward", "ortho", "forward"}, optional
+        Normalization mode (see `numpy.fft`). Default is "backward".
+        Indicates which direction of the forward/backward pair of transforms
+        is scaled and with what normalization factor.
+
+        .. versionadded:: 1.20.0
+
+            The "backward", "forward" values were added.
+
+    out : complex ndarray, optional
+        If provided, the result will be placed in this array. It should be
+        of the appropriate shape and dtype for all axes (and hence only the
+        last axis can have ``s`` not equal to the shape at that axis).
+
+        .. versionadded:: 2.0.0
+
+    Returns
+    -------
+    out : complex ndarray
+        The truncated or zero-padded input, transformed along the axes
+        indicated by `axes`, or the last two axes if `axes` is not given.
+
+    Raises
+    ------
+    ValueError
+        If `s` and `axes` have different length, or `axes` not given and
+        ``len(s) != 2``.
+    IndexError
+        If an element of `axes` is larger than than the number of axes of `a`.
+
+    See Also
+    --------
+    numpy.fft : Overall view of discrete Fourier transforms, with definitions
+         and conventions used.
+    ifft2 : The inverse two-dimensional FFT.
+    fft : The one-dimensional FFT.
+    fftn : The *n*-dimensional FFT.
+    fftshift : Shifts zero-frequency terms to the center of the array.
+        For two-dimensional input, swaps first and third quadrants, and second
+        and fourth quadrants.
+
+    Notes
+    -----
+    `fft2` is just `fftn` with a different default for `axes`.
+
+    The output, analogously to `fft`, contains the term for zero frequency in
+    the low-order corner of the transformed axes, the positive frequency terms
+    in the first half of these axes, the term for the Nyquist frequency in the
+    middle of the axes and the negative frequency terms in the second half of
+    the axes, in order of decreasingly negative frequency.
+
+    See `fftn` for details and a plotting example, and `numpy.fft` for
+    definitions and conventions used.
+
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> a = np.mgrid[:5, :5][0]
+    >>> np.fft.fft2(a)
+    array([[ 50.  +0.j        ,   0.  +0.j        ,   0.  +0.j        , # may vary
+              0.  +0.j        ,   0.  +0.j        ],
+           [-12.5+17.20477401j,   0.  +0.j        ,   0.  +0.j        ,
+              0.  +0.j        ,   0.  +0.j        ],
+           [-12.5 +4.0614962j ,   0.  +0.j        ,   0.  +0.j        ,
+              0.  +0.j        ,   0.  +0.j        ],
+           [-12.5 -4.0614962j ,   0.  +0.j        ,   0.  +0.j        ,
+              0.  +0.j        ,   0.  +0.j        ],
+           [-12.5-17.20477401j,   0.  +0.j        ,   0.  +0.j        ,
+              0.  +0.j        ,   0.  +0.j        ]])
+
+    """
+    return _raw_fftnd(a, s, axes, fft, norm, out=out)
+
+
+@array_function_dispatch(_fftn_dispatcher)
+def ifft2(a, s=None, axes=(-2, -1), norm=None, out=None):
+    """
+    Compute the 2-dimensional inverse discrete Fourier Transform.
+
+    This function computes the inverse of the 2-dimensional discrete Fourier
+    Transform over any number of axes in an M-dimensional array by means of
+    the Fast Fourier Transform (FFT).  In other words, ``ifft2(fft2(a)) == a``
+    to within numerical accuracy.  By default, the inverse transform is
+    computed over the last two axes of the input array.
+
+    The input, analogously to `ifft`, should be ordered in the same way as is
+    returned by `fft2`, i.e. it should have the term for zero frequency
+    in the low-order corner of the two axes, the positive frequency terms in
+    the first half of these axes, the term for the Nyquist frequency in the
+    middle of the axes and the negative frequency terms in the second half of
+    both axes, in order of decreasingly negative frequency.
+
+    Parameters
+    ----------
+    a : array_like
+        Input array, can be complex.
+    s : sequence of ints, optional
+        Shape (length of each axis) of the output (``s[0]`` refers to axis 0,
+        ``s[1]`` to axis 1, etc.).  This corresponds to `n` for ``ifft(x, n)``.
+        Along each axis, if the given shape is smaller than that of the input,
+        the input is cropped. If it is larger, the input is padded with zeros.
+
+        .. versionchanged:: 2.0
+
+            If it is ``-1``, the whole input is used (no padding/trimming).
+
+        If `s` is not given, the shape of the input along the axes specified
+        by `axes` is used.  See notes for issue on `ifft` zero padding.
+
+        .. deprecated:: 2.0
+
+            If `s` is not ``None``, `axes` must not be ``None`` either.
+
+        .. deprecated:: 2.0
+
+            `s` must contain only ``int`` s, not ``None`` values. ``None``
+            values currently mean that the default value for ``n`` is used
+            in the corresponding 1-D transform, but this behaviour is
+            deprecated.
+
+    axes : sequence of ints, optional
+        Axes over which to compute the FFT.  If not given, the last two
+        axes are used.  A repeated index in `axes` means the transform over
+        that axis is performed multiple times.  A one-element sequence means
+        that a one-dimensional FFT is performed. Default: ``(-2, -1)``.
+
+        .. deprecated:: 2.0
+
+            If `s` is specified, the corresponding `axes` to be transformed
+            must not be ``None``.
+
+    norm : {"backward", "ortho", "forward"}, optional
+        Normalization mode (see `numpy.fft`). Default is "backward".
+        Indicates which direction of the forward/backward pair of transforms
+        is scaled and with what normalization factor.
+
+        .. versionadded:: 1.20.0
+
+            The "backward", "forward" values were added.
+
+    out : complex ndarray, optional
+        If provided, the result will be placed in this array. It should be
+        of the appropriate shape and dtype for all axes (and hence is
+        incompatible with passing in all but the trivial ``s``).
+
+        .. versionadded:: 2.0.0
+
+    Returns
+    -------
+    out : complex ndarray
+        The truncated or zero-padded input, transformed along the axes
+        indicated by `axes`, or the last two axes if `axes` is not given.
+
+    Raises
+    ------
+    ValueError
+        If `s` and `axes` have different length, or `axes` not given and
+        ``len(s) != 2``.
+    IndexError
+        If an element of `axes` is larger than than the number of axes of `a`.
+
+    See Also
+    --------
+    numpy.fft : Overall view of discrete Fourier transforms, with definitions
+         and conventions used.
+    fft2 : The forward 2-dimensional FFT, of which `ifft2` is the inverse.
+    ifftn : The inverse of the *n*-dimensional FFT.
+    fft : The one-dimensional FFT.
+    ifft : The one-dimensional inverse FFT.
+
+    Notes
+    -----
+    `ifft2` is just `ifftn` with a different default for `axes`.
+
+    See `ifftn` for details and a plotting example, and `numpy.fft` for
+    definition and conventions used.
+
+    Zero-padding, analogously with `ifft`, is performed by appending zeros to
+    the input along the specified dimension.  Although this is the common
+    approach, it might lead to surprising results.  If another form of zero
+    padding is desired, it must be performed before `ifft2` is called.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> a = 4 * np.eye(4)
+    >>> np.fft.ifft2(a)
+    array([[1.+0.j,  0.+0.j,  0.+0.j,  0.+0.j], # may vary
+           [0.+0.j,  0.+0.j,  0.+0.j,  1.+0.j],
+           [0.+0.j,  0.+0.j,  1.+0.j,  0.+0.j],
+           [0.+0.j,  1.+0.j,  0.+0.j,  0.+0.j]])
+
+    """
+    return _raw_fftnd(a, s, axes, ifft, norm, out=None)
+
+
+@array_function_dispatch(_fftn_dispatcher)
+def rfftn(a, s=None, axes=None, norm=None, out=None):
+    """
+    Compute the N-dimensional discrete Fourier Transform for real input.
+
+    This function computes the N-dimensional discrete Fourier Transform over
+    any number of axes in an M-dimensional real array by means of the Fast
+    Fourier Transform (FFT).  By default, all axes are transformed, with the
+    real transform performed over the last axis, while the remaining
+    transforms are complex.
+
+    Parameters
+    ----------
+    a : array_like
+        Input array, taken to be real.
+    s : sequence of ints, optional
+        Shape (length along each transformed axis) to use from the input.
+        (``s[0]`` refers to axis 0, ``s[1]`` to axis 1, etc.).
+        The final element of `s` corresponds to `n` for ``rfft(x, n)``, while
+        for the remaining axes, it corresponds to `n` for ``fft(x, n)``.
+        Along any axis, if the given shape is smaller than that of the input,
+        the input is cropped. If it is larger, the input is padded with zeros.
+
+        .. versionchanged:: 2.0
+
+            If it is ``-1``, the whole input is used (no padding/trimming).
+
+        If `s` is not given, the shape of the input along the axes specified
+        by `axes` is used.
+
+        .. deprecated:: 2.0
+
+            If `s` is not ``None``, `axes` must not be ``None`` either.
+
+        .. deprecated:: 2.0
+
+            `s` must contain only ``int`` s, not ``None`` values. ``None``
+            values currently mean that the default value for ``n`` is used
+            in the corresponding 1-D transform, but this behaviour is
+            deprecated.
+
+    axes : sequence of ints, optional
+        Axes over which to compute the FFT.  If not given, the last ``len(s)``
+        axes are used, or all axes if `s` is also not specified.
+
+        .. deprecated:: 2.0
+
+            If `s` is specified, the corresponding `axes` to be transformed
+            must be explicitly specified too.
+
+    norm : {"backward", "ortho", "forward"}, optional
+        Normalization mode (see `numpy.fft`). Default is "backward".
+        Indicates which direction of the forward/backward pair of transforms
+        is scaled and with what normalization factor.
+
+        .. versionadded:: 1.20.0
+
+            The "backward", "forward" values were added.
+
+    out : complex ndarray, optional
+        If provided, the result will be placed in this array. It should be
+        of the appropriate shape and dtype for all axes (and hence is
+        incompatible with passing in all but the trivial ``s``).
+
+        .. versionadded:: 2.0.0
+
+    Returns
+    -------
+    out : complex ndarray
+        The truncated or zero-padded input, transformed along the axes
+        indicated by `axes`, or by a combination of `s` and `a`,
+        as explained in the parameters section above.
+        The length of the last axis transformed will be ``s[-1]//2+1``,
+        while the remaining transformed axes will have lengths according to
+        `s`, or unchanged from the input.
+
+    Raises
+    ------
+    ValueError
+        If `s` and `axes` have different length.
+    IndexError
+        If an element of `axes` is larger than than the number of axes of `a`.
+
+    See Also
+    --------
+    irfftn : The inverse of `rfftn`, i.e. the inverse of the n-dimensional FFT
+         of real input.
+    fft : The one-dimensional FFT, with definitions and conventions used.
+    rfft : The one-dimensional FFT of real input.
+    fftn : The n-dimensional FFT.
+    rfft2 : The two-dimensional FFT of real input.
+
+    Notes
+    -----
+    The transform for real input is performed over the last transformation
+    axis, as by `rfft`, then the transform over the remaining axes is
+    performed as by `fftn`.  The order of the output is as for `rfft` for the
+    final transformation axis, and as for `fftn` for the remaining
+    transformation axes.
+
+    See `fft` for details, definitions and conventions used.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> a = np.ones((2, 2, 2))
+    >>> np.fft.rfftn(a)
+    array([[[8.+0.j,  0.+0.j], # may vary
+            [0.+0.j,  0.+0.j]],
+           [[0.+0.j,  0.+0.j],
+            [0.+0.j,  0.+0.j]]])
+
+    >>> np.fft.rfftn(a, axes=(2, 0))
+    array([[[4.+0.j,  0.+0.j], # may vary
+            [4.+0.j,  0.+0.j]],
+           [[0.+0.j,  0.+0.j],
+            [0.+0.j,  0.+0.j]]])
+
+    """
+    a = asarray(a)
+    s, axes = _cook_nd_args(a, s, axes)
+    a = rfft(a, s[-1], axes[-1], norm, out=out)
+    for ii in range(len(axes) - 2, -1, -1):
+        a = fft(a, s[ii], axes[ii], norm, out=out)
+    return a
+
+
+@array_function_dispatch(_fftn_dispatcher)
+def rfft2(a, s=None, axes=(-2, -1), norm=None, out=None):
+    """
+    Compute the 2-dimensional FFT of a real array.
+
+    Parameters
+    ----------
+    a : array
+        Input array, taken to be real.
+    s : sequence of ints, optional
+        Shape of the FFT.
+
+        .. versionchanged:: 2.0
+
+            If it is ``-1``, the whole input is used (no padding/trimming).
+
+        .. deprecated:: 2.0
+
+            If `s` is not ``None``, `axes` must not be ``None`` either.
+
+        .. deprecated:: 2.0
+
+            `s` must contain only ``int`` s, not ``None`` values. ``None``
+            values currently mean that the default value for ``n`` is used
+            in the corresponding 1-D transform, but this behaviour is
+            deprecated.
+
+    axes : sequence of ints, optional
+        Axes over which to compute the FFT. Default: ``(-2, -1)``.
+
+        .. deprecated:: 2.0
+
+            If `s` is specified, the corresponding `axes` to be transformed
+            must not be ``None``.
+
+    norm : {"backward", "ortho", "forward"}, optional
+        Normalization mode (see `numpy.fft`). Default is "backward".
+        Indicates which direction of the forward/backward pair of transforms
+        is scaled and with what normalization factor.
+
+        .. versionadded:: 1.20.0
+
+            The "backward", "forward" values were added.
+
+    out : complex ndarray, optional
+        If provided, the result will be placed in this array. It should be
+        of the appropriate shape and dtype for the last inverse transform.
+        incompatible with passing in all but the trivial ``s``).
+
+        .. versionadded:: 2.0.0
+
+    Returns
+    -------
+    out : ndarray
+        The result of the real 2-D FFT.
+
+    See Also
+    --------
+    rfftn : Compute the N-dimensional discrete Fourier Transform for real
+            input.
+
+    Notes
+    -----
+    This is really just `rfftn` with different default behavior.
+    For more details see `rfftn`.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> a = np.mgrid[:5, :5][0]
+    >>> np.fft.rfft2(a)
+    array([[ 50.  +0.j        ,   0.  +0.j        ,   0.  +0.j        ],
+           [-12.5+17.20477401j,   0.  +0.j        ,   0.  +0.j        ],
+           [-12.5 +4.0614962j ,   0.  +0.j        ,   0.  +0.j        ],
+           [-12.5 -4.0614962j ,   0.  +0.j        ,   0.  +0.j        ],
+           [-12.5-17.20477401j,   0.  +0.j        ,   0.  +0.j        ]])
+    """
+    return rfftn(a, s, axes, norm, out=out)
+
+
+@array_function_dispatch(_fftn_dispatcher)
+def irfftn(a, s=None, axes=None, norm=None, out=None):
+    """
+    Computes the inverse of `rfftn`.
+
+    This function computes the inverse of the N-dimensional discrete
+    Fourier Transform for real input over any number of axes in an
+    M-dimensional array by means of the Fast Fourier Transform (FFT).  In
+    other words, ``irfftn(rfftn(a), a.shape) == a`` to within numerical
+    accuracy. (The ``a.shape`` is necessary like ``len(a)`` is for `irfft`,
+    and for the same reason.)
+
+    The input should be ordered in the same way as is returned by `rfftn`,
+    i.e. as for `irfft` for the final transformation axis, and as for `ifftn`
+    along all the other axes.
+
+    Parameters
+    ----------
+    a : array_like
+        Input array.
+    s : sequence of ints, optional
+        Shape (length of each transformed axis) of the output
+        (``s[0]`` refers to axis 0, ``s[1]`` to axis 1, etc.). `s` is also the
+        number of input points used along this axis, except for the last axis,
+        where ``s[-1]//2+1`` points of the input are used.
+        Along any axis, if the shape indicated by `s` is smaller than that of
+        the input, the input is cropped.  If it is larger, the input is padded
+        with zeros.
+
+        .. versionchanged:: 2.0
+
+            If it is ``-1``, the whole input is used (no padding/trimming).
+
+        If `s` is not given, the shape of the input along the axes
+        specified by axes is used. Except for the last axis which is taken to
+        be ``2*(m-1)`` where ``m`` is the length of the input along that axis.
+
+        .. deprecated:: 2.0
+
+            If `s` is not ``None``, `axes` must not be ``None`` either.
+
+        .. deprecated:: 2.0
+
+            `s` must contain only ``int`` s, not ``None`` values. ``None``
+            values currently mean that the default value for ``n`` is used
+            in the corresponding 1-D transform, but this behaviour is
+            deprecated.
+
+    axes : sequence of ints, optional
+        Axes over which to compute the inverse FFT. If not given, the last
+        `len(s)` axes are used, or all axes if `s` is also not specified.
+        Repeated indices in `axes` means that the inverse transform over that
+        axis is performed multiple times.
+
+        .. deprecated:: 2.0
+
+            If `s` is specified, the corresponding `axes` to be transformed
+            must be explicitly specified too.
+
+    norm : {"backward", "ortho", "forward"}, optional
+        Normalization mode (see `numpy.fft`). Default is "backward".
+        Indicates which direction of the forward/backward pair of transforms
+        is scaled and with what normalization factor.
+
+        .. versionadded:: 1.20.0
+
+            The "backward", "forward" values were added.
+
+    out : ndarray, optional
+        If provided, the result will be placed in this array. It should be
+        of the appropriate shape and dtype for the last transformation.
+
+        .. versionadded:: 2.0.0
+
+    Returns
+    -------
+    out : ndarray
+        The truncated or zero-padded input, transformed along the axes
+        indicated by `axes`, or by a combination of `s` or `a`,
+        as explained in the parameters section above.
+        The length of each transformed axis is as given by the corresponding
+        element of `s`, or the length of the input in every axis except for the
+        last one if `s` is not given.  In the final transformed axis the length
+        of the output when `s` is not given is ``2*(m-1)`` where ``m`` is the
+        length of the final transformed axis of the input.  To get an odd
+        number of output points in the final axis, `s` must be specified.
+
+    Raises
+    ------
+    ValueError
+        If `s` and `axes` have different length.
+    IndexError
+        If an element of `axes` is larger than than the number of axes of `a`.
+
+    See Also
+    --------
+    rfftn : The forward n-dimensional FFT of real input,
+            of which `ifftn` is the inverse.
+    fft : The one-dimensional FFT, with definitions and conventions used.
+    irfft : The inverse of the one-dimensional FFT of real input.
+    irfft2 : The inverse of the two-dimensional FFT of real input.
+
+    Notes
+    -----
+    See `fft` for definitions and conventions used.
+
+    See `rfft` for definitions and conventions used for real input.
+
+    The correct interpretation of the hermitian input depends on the shape of
+    the original data, as given by `s`. This is because each input shape could
+    correspond to either an odd or even length signal. By default, `irfftn`
+    assumes an even output length which puts the last entry at the Nyquist
+    frequency; aliasing with its symmetric counterpart. When performing the
+    final complex to real transform, the last value is thus treated as purely
+    real. To avoid losing information, the correct shape of the real input
+    **must** be given.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> a = np.zeros((3, 2, 2))
+    >>> a[0, 0, 0] = 3 * 2 * 2
+    >>> np.fft.irfftn(a)
+    array([[[1.,  1.],
+            [1.,  1.]],
+           [[1.,  1.],
+            [1.,  1.]],
+           [[1.,  1.],
+            [1.,  1.]]])
+
+    """
+    a = asarray(a)
+    s, axes = _cook_nd_args(a, s, axes, invreal=1)
+    for ii in range(len(axes) - 1):
+        a = ifft(a, s[ii], axes[ii], norm)
+    a = irfft(a, s[-1], axes[-1], norm, out=out)
+    return a
+
+
+@array_function_dispatch(_fftn_dispatcher)
+def irfft2(a, s=None, axes=(-2, -1), norm=None, out=None):
+    """
+    Computes the inverse of `rfft2`.
+
+    Parameters
+    ----------
+    a : array_like
+        The input array
+    s : sequence of ints, optional
+        Shape of the real output to the inverse FFT.
+
+        .. versionchanged:: 2.0
+
+            If it is ``-1``, the whole input is used (no padding/trimming).
+
+        .. deprecated:: 2.0
+
+            If `s` is not ``None``, `axes` must not be ``None`` either.
+
+        .. deprecated:: 2.0
+
+            `s` must contain only ``int`` s, not ``None`` values. ``None``
+            values currently mean that the default value for ``n`` is used
+            in the corresponding 1-D transform, but this behaviour is
+            deprecated.
+
+    axes : sequence of ints, optional
+        The axes over which to compute the inverse fft.
+        Default: ``(-2, -1)``, the last two axes.
+
+        .. deprecated:: 2.0
+
+            If `s` is specified, the corresponding `axes` to be transformed
+            must not be ``None``.
+
+    norm : {"backward", "ortho", "forward"}, optional
+        Normalization mode (see `numpy.fft`). Default is "backward".
+        Indicates which direction of the forward/backward pair of transforms
+        is scaled and with what normalization factor.
+
+        .. versionadded:: 1.20.0
+
+            The "backward", "forward" values were added.
+
+    out : ndarray, optional
+        If provided, the result will be placed in this array. It should be
+        of the appropriate shape and dtype for the last transformation.
+
+        .. versionadded:: 2.0.0
+
+    Returns
+    -------
+    out : ndarray
+        The result of the inverse real 2-D FFT.
+
+    See Also
+    --------
+    rfft2 : The forward two-dimensional FFT of real input,
+            of which `irfft2` is the inverse.
+    rfft : The one-dimensional FFT for real input.
+    irfft : The inverse of the one-dimensional FFT of real input.
+    irfftn : Compute the inverse of the N-dimensional FFT of real input.
+
+    Notes
+    -----
+    This is really `irfftn` with different defaults.
+    For more details see `irfftn`.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> a = np.mgrid[:5, :5][0]
+    >>> A = np.fft.rfft2(a)
+    >>> np.fft.irfft2(A, s=a.shape)
+    array([[0., 0., 0., 0., 0.],
+           [1., 1., 1., 1., 1.],
+           [2., 2., 2., 2., 2.],
+           [3., 3., 3., 3., 3.],
+           [4., 4., 4., 4., 4.]])
+    """
+    return irfftn(a, s, axes, norm, out=None)
diff --git a/.venv/lib/python3.12/site-packages/numpy/fft/_pocketfft.pyi b/.venv/lib/python3.12/site-packages/numpy/fft/_pocketfft.pyi
new file mode 100644
index 0000000..215cf14
--- /dev/null
+++ b/.venv/lib/python3.12/site-packages/numpy/fft/_pocketfft.pyi
@@ -0,0 +1,138 @@
+from collections.abc import Sequence
+from typing import Literal as L
+from typing import TypeAlias
+
+from numpy import complex128, float64
+from numpy._typing import ArrayLike, NDArray, _ArrayLikeNumber_co
+
+__all__ = [
+    "fft",
+    "ifft",
+    "rfft",
+    "irfft",
+    "hfft",
+    "ihfft",
+    "rfftn",
+    "irfftn",
+    "rfft2",
+    "irfft2",
+    "fft2",
+    "ifft2",
+    "fftn",
+    "ifftn",
+]
+
+_NormKind: TypeAlias = L["backward", "ortho", "forward"] | None
+
+def fft(
+    a: ArrayLike,
+    n: int | None = ...,
+    axis: int = ...,
+    norm: _NormKind = ...,
+    out: NDArray[complex128] | None = ...,
+) -> NDArray[complex128]: ...
+
+def ifft(
+    a: ArrayLike,
+    n: int | None = ...,
+    axis: int = ...,
+    norm: _NormKind = ...,
+    out: NDArray[complex128] | None = ...,
+) -> NDArray[complex128]: ...
+
+def rfft(
+    a: ArrayLike,
+    n: int | None = ...,
+    axis: int = ...,
+    norm: _NormKind = ...,
+    out: NDArray[complex128] | None = ...,
+) -> NDArray[complex128]: ...
+
+def irfft(
+    a: ArrayLike,
+    n: int | None = ...,
+    axis: int = ...,
+    norm: _NormKind = ...,
+    out: NDArray[float64] | None = ...,
+) -> NDArray[float64]: ...
+
+# Input array must be compatible with `np.conjugate`
+def hfft(
+    a: _ArrayLikeNumber_co,
+    n: int | None = ...,
+    axis: int = ...,
+    norm: _NormKind = ...,
+    out: NDArray[float64] | None = ...,
+) -> NDArray[float64]: ...
+
+def ihfft(
+    a: ArrayLike,
+    n: int | None = ...,
+    axis: int = ...,
+    norm: _NormKind = ...,
+    out: NDArray[complex128] | None = ...,
+) -> NDArray[complex128]: ...
+
+def fftn(
+    a: ArrayLike,
+    s: Sequence[int] | None = ...,
+    axes: Sequence[int] | None = ...,
+    norm: _NormKind = ...,
+    out: NDArray[complex128] | None = ...,
+) -> NDArray[complex128]: ...
+
+def ifftn(
+    a: ArrayLike,
+    s: Sequence[int] | None = ...,
+    axes: Sequence[int] | None = ...,
+    norm: _NormKind = ...,
+    out: NDArray[complex128] | None = ...,
+) -> NDArray[complex128]: ...
+
+def rfftn(
+    a: ArrayLike,
+    s: Sequence[int] | None = ...,
+    axes: Sequence[int] | None = ...,
+    norm: _NormKind = ...,
+    out: NDArray[complex128] | None = ...,
+) -> NDArray[complex128]: ...
+
+def irfftn(
+    a: ArrayLike,
+    s: Sequence[int] | None = ...,
+    axes: Sequence[int] | None = ...,
+    norm: _NormKind = ...,
+    out: NDArray[float64] | None = ...,
+) -> NDArray[float64]: ...
+
+def fft2(
+    a: ArrayLike,
+    s: Sequence[int] | None = ...,
+    axes: Sequence[int] | None = ...,
+    norm: _NormKind = ...,
+    out: NDArray[complex128] | None = ...,
+) -> NDArray[complex128]: ...
+
+def ifft2(
+    a: ArrayLike,
+    s: Sequence[int] | None = ...,
+    axes: Sequence[int] | None = ...,
+    norm: _NormKind = ...,
+    out: NDArray[complex128] | None = ...,
+) -> NDArray[complex128]: ...
+
+def rfft2(
+    a: ArrayLike,
+    s: Sequence[int] | None = ...,
+    axes: Sequence[int] | None = ...,
+    norm: _NormKind = ...,
+    out: NDArray[complex128] | None = ...,
+) -> NDArray[complex128]: ...
+
+def irfft2(
+    a: ArrayLike,
+    s: Sequence[int] | None = ...,
+    axes: Sequence[int] | None = ...,
+    norm: _NormKind = ...,
+    out: NDArray[float64] | None = ...,
+) -> NDArray[float64]: ...
diff --git a/.venv/lib/python3.12/site-packages/numpy/fft/_pocketfft_umath.cpython-312-x86_64-linux-gnu.so b/.venv/lib/python3.12/site-packages/numpy/fft/_pocketfft_umath.cpython-312-x86_64-linux-gnu.so
new file mode 100755
index 0000000..df57850
--- /dev/null
+++ b/.venv/lib/python3.12/site-packages/numpy/fft/_pocketfft_umath.cpython-312-x86_64-linux-gnu.so
diff --git a/.venv/lib/python3.12/site-packages/numpy/fft/helper.py b/.venv/lib/python3.12/site-packages/numpy/fft/helper.py
new file mode 100644
index 0000000..08d5662
--- /dev/null
+++ b/.venv/lib/python3.12/site-packages/numpy/fft/helper.py
@@ -0,0 +1,17 @@
+def __getattr__(attr_name):
+    import warnings
+
+    from numpy.fft import _helper
+    ret = getattr(_helper, attr_name, None)
+    if ret is None:
+        raise AttributeError(
+            f"module 'numpy.fft.helper' has no attribute {attr_name}")
+    warnings.warn(
+        "The numpy.fft.helper has been made private and renamed to "
+        "numpy.fft._helper. All four functions exported by it (i.e. fftshift, "
+        "ifftshift, fftfreq, rfftfreq) are available from numpy.fft. "
+        f"Please use numpy.fft.{attr_name} instead.",
+        DeprecationWarning,
+        stacklevel=3
+    )
+    return ret
diff --git a/.venv/lib/python3.12/site-packages/numpy/fft/helper.pyi b/.venv/lib/python3.12/site-packages/numpy/fft/helper.pyi
new file mode 100644
index 0000000..887cbe7
--- /dev/null
+++ b/.venv/lib/python3.12/site-packages/numpy/fft/helper.pyi
@@ -0,0 +1,22 @@
+from typing import Any
+from typing import Literal as L
+
+from typing_extensions import deprecated
+
+import numpy as np
+from numpy._typing import ArrayLike, NDArray, _ShapeLike
+
+from ._helper import integer_types as integer_types
+
+__all__ = ["fftfreq", "fftshift", "ifftshift", "rfftfreq"]
+
+###
+
+@deprecated("Please use `numpy.fft.fftshift` instead.")
+def fftshift(x: ArrayLike, axes: _ShapeLike | None = None) -> NDArray[Any]: ...
+@deprecated("Please use `numpy.fft.ifftshift` instead.")
+def ifftshift(x: ArrayLike, axes: _ShapeLike | None = None) -> NDArray[Any]: ...
+@deprecated("Please use `numpy.fft.fftfreq` instead.")
+def fftfreq(n: int | np.integer, d: ArrayLike = 1.0, device: L["cpu"] | None = None) -> NDArray[Any]: ...
+@deprecated("Please use `numpy.fft.rfftfreq` instead.")
+def rfftfreq(n: int | np.integer, d: ArrayLike = 1.0, device: L["cpu"] | None = None) -> NDArray[Any]: ...
diff --git a/.venv/lib/python3.12/site-packages/numpy/fft/tests/__init__.py b/.venv/lib/python3.12/site-packages/numpy/fft/tests/__init__.py
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/.venv/lib/python3.12/site-packages/numpy/fft/tests/__init__.py
diff --git a/.venv/lib/python3.12/site-packages/numpy/fft/tests/__pycache__/__init__.cpython-312.pyc b/.venv/lib/python3.12/site-packages/numpy/fft/tests/__pycache__/__init__.cpython-312.pyc
new file mode 100644
index 0000000..f2fb847
--- /dev/null
+++ b/.venv/lib/python3.12/site-packages/numpy/fft/tests/__pycache__/__init__.cpython-312.pyc
diff --git a/.venv/lib/python3.12/site-packages/numpy/fft/tests/__pycache__/test_helper.cpython-312.pyc b/.venv/lib/python3.12/site-packages/numpy/fft/tests/__pycache__/test_helper.cpython-312.pyc
new file mode 100644
index 0000000..a1e4182
--- /dev/null
+++ b/.venv/lib/python3.12/site-packages/numpy/fft/tests/__pycache__/test_helper.cpython-312.pyc
diff --git a/.venv/lib/python3.12/site-packages/numpy/fft/tests/__pycache__/test_pocketfft.cpython-312.pyc b/.venv/lib/python3.12/site-packages/numpy/fft/tests/__pycache__/test_pocketfft.cpython-312.pyc
new file mode 100644
index 0000000..853bb4e
--- /dev/null
+++ b/.venv/lib/python3.12/site-packages/numpy/fft/tests/__pycache__/test_pocketfft.cpython-312.pyc
diff --git a/.venv/lib/python3.12/site-packages/numpy/fft/tests/test_helper.py b/.venv/lib/python3.12/site-packages/numpy/fft/tests/test_helper.py
new file mode 100644
index 0000000..c02a736
--- /dev/null
+++ b/.venv/lib/python3.12/site-packages/numpy/fft/tests/test_helper.py
@@ -0,0 +1,167 @@
+"""Test functions for fftpack.helper module
+
+Copied from fftpack.helper by Pearu Peterson, October 2005
+
+"""
+import numpy as np
+from numpy import fft, pi
+from numpy.testing import assert_array_almost_equal
+
+
+class TestFFTShift:
+
+    def test_definition(self):
+        x = [0, 1, 2, 3, 4, -4, -3, -2, -1]
+        y = [-4, -3, -2, -1, 0, 1, 2, 3, 4]
+        assert_array_almost_equal(fft.fftshift(x), y)
+        assert_array_almost_equal(fft.ifftshift(y), x)
+        x = [0, 1, 2, 3, 4, -5, -4, -3, -2, -1]
+        y = [-5, -4, -3, -2, -1, 0, 1, 2, 3, 4]
+        assert_array_almost_equal(fft.fftshift(x), y)
+        assert_array_almost_equal(fft.ifftshift(y), x)
+
+    def test_inverse(self):
+        for n in [1, 4, 9, 100, 211]:
+            x = np.random.random((n,))
+            assert_array_almost_equal(fft.ifftshift(fft.fftshift(x)), x)
+
+    def test_axes_keyword(self):
+        freqs = [[0, 1, 2], [3, 4, -4], [-3, -2, -1]]
+        shifted = [[-1, -3, -2], [2, 0, 1], [-4, 3, 4]]
+        assert_array_almost_equal(fft.fftshift(freqs, axes=(0, 1)), shifted)
+        assert_array_almost_equal(fft.fftshift(freqs, axes=0),
+                                  fft.fftshift(freqs, axes=(0,)))
+        assert_array_almost_equal(fft.ifftshift(shifted, axes=(0, 1)), freqs)
+        assert_array_almost_equal(fft.ifftshift(shifted, axes=0),
+                                  fft.ifftshift(shifted, axes=(0,)))
+
+        assert_array_almost_equal(fft.fftshift(freqs), shifted)
+        assert_array_almost_equal(fft.ifftshift(shifted), freqs)
+
+    def test_uneven_dims(self):
+        """ Test 2D input, which has uneven dimension sizes """
+        freqs = [
+            [0, 1],
+            [2, 3],
+            [4, 5]
+        ]
+
+        # shift in dimension 0
+        shift_dim0 = [
+            [4, 5],
+            [0, 1],
+            [2, 3]
+        ]
+        assert_array_almost_equal(fft.fftshift(freqs, axes=0), shift_dim0)
+        assert_array_almost_equal(fft.ifftshift(shift_dim0, axes=0), freqs)
+        assert_array_almost_equal(fft.fftshift(freqs, axes=(0,)), shift_dim0)
+        assert_array_almost_equal(fft.ifftshift(shift_dim0, axes=[0]), freqs)
+
+        # shift in dimension 1
+        shift_dim1 = [
+            [1, 0],
+            [3, 2],
+            [5, 4]
+        ]
+        assert_array_almost_equal(fft.fftshift(freqs, axes=1), shift_dim1)
+        assert_array_almost_equal(fft.ifftshift(shift_dim1, axes=1), freqs)
+
+        # shift in both dimensions
+        shift_dim_both = [
+            [5, 4],
+            [1, 0],
+            [3, 2]
+        ]
+        assert_array_almost_equal(fft.fftshift(freqs, axes=(0, 1)), shift_dim_both)
+        assert_array_almost_equal(fft.ifftshift(shift_dim_both, axes=(0, 1)), freqs)
+        assert_array_almost_equal(fft.fftshift(freqs, axes=[0, 1]), shift_dim_both)
+        assert_array_almost_equal(fft.ifftshift(shift_dim_both, axes=[0, 1]), freqs)
+
+        # axes=None (default) shift in all dimensions
+        assert_array_almost_equal(fft.fftshift(freqs, axes=None), shift_dim_both)
+        assert_array_almost_equal(fft.ifftshift(shift_dim_both, axes=None), freqs)
+        assert_array_almost_equal(fft.fftshift(freqs), shift_dim_both)
+        assert_array_almost_equal(fft.ifftshift(shift_dim_both), freqs)
+
+    def test_equal_to_original(self):
+        """ Test the new (>=v1.15) and old implementations are equal (see #10073) """
+        from numpy._core import arange, asarray, concatenate, take
+
+        def original_fftshift(x, axes=None):
+            """ How fftshift was implemented in v1.14"""
+            tmp = asarray(x)
+            ndim = tmp.ndim
+            if axes is None:
+                axes = list(range(ndim))
+            elif isinstance(axes, int):
+                axes = (axes,)
+            y = tmp
+            for k in axes:
+                n = tmp.shape[k]
+                p2 = (n + 1) // 2
+                mylist = concatenate((arange(p2, n), arange(p2)))
+                y = take(y, mylist, k)
+            return y
+
+        def original_ifftshift(x, axes=None):
+            """ How ifftshift was implemented in v1.14 """
+            tmp = asarray(x)
+            ndim = tmp.ndim
+            if axes is None:
+                axes = list(range(ndim))
+            elif isinstance(axes, int):
+                axes = (axes,)
+            y = tmp
+            for k in axes:
+                n = tmp.shape[k]
+                p2 = n - (n + 1) // 2
+                mylist = concatenate((arange(p2, n), arange(p2)))
+                y = take(y, mylist, k)
+            return y
+
+        # create possible 2d array combinations and try all possible keywords
+        # compare output to original functions
+        for i in range(16):
+            for j in range(16):
+                for axes_keyword in [0, 1, None, (0,), (0, 1)]:
+                    inp = np.random.rand(i, j)
+
+                    assert_array_almost_equal(fft.fftshift(inp, axes_keyword),
+                                              original_fftshift(inp, axes_keyword))
+
+                    assert_array_almost_equal(fft.ifftshift(inp, axes_keyword),
+                                              original_ifftshift(inp, axes_keyword))
+
+
+class TestFFTFreq:
+
+    def test_definition(self):
+        x = [0, 1, 2, 3, 4, -4, -3, -2, -1]
+        assert_array_almost_equal(9 * fft.fftfreq(9), x)
+        assert_array_almost_equal(9 * pi * fft.fftfreq(9, pi), x)
+        x = [0, 1, 2, 3, 4, -5, -4, -3, -2, -1]
+        assert_array_almost_equal(10 * fft.fftfreq(10), x)
+        assert_array_almost_equal(10 * pi * fft.fftfreq(10, pi), x)
+
+
+class TestRFFTFreq:
+
+    def test_definition(self):
+        x = [0, 1, 2, 3, 4]
+        assert_array_almost_equal(9 * fft.rfftfreq(9), x)
+        assert_array_almost_equal(9 * pi * fft.rfftfreq(9, pi), x)
+        x = [0, 1, 2, 3, 4, 5]
+        assert_array_almost_equal(10 * fft.rfftfreq(10), x)
+        assert_array_almost_equal(10 * pi * fft.rfftfreq(10, pi), x)
+
+
+class TestIRFFTN:
+
+    def test_not_last_axis_success(self):
+        ar, ai = np.random.random((2, 16, 8, 32))
+        a = ar + 1j * ai
+
+        axes = (-2,)
+
+        # Should not raise error
+        fft.irfftn(a, axes=axes)
diff --git a/.venv/lib/python3.12/site-packages/numpy/fft/tests/test_pocketfft.py b/.venv/lib/python3.12/site-packages/numpy/fft/tests/test_pocketfft.py
new file mode 100644
index 0000000..0211818
--- /dev/null
+++ b/.venv/lib/python3.12/site-packages/numpy/fft/tests/test_pocketfft.py
@@ -0,0 +1,589 @@
+import queue
+import threading
+
+import pytest
+
+import numpy as np
+from numpy.random import random
+from numpy.testing import IS_WASM, assert_allclose, assert_array_equal, assert_raises
+
+
+def fft1(x):
+    L = len(x)
+    phase = -2j * np.pi * (np.arange(L) / L)
+    phase = np.arange(L).reshape(-1, 1) * phase
+    return np.sum(x * np.exp(phase), axis=1)
+
+
+class TestFFTShift:
+
+    def test_fft_n(self):
+        assert_raises(ValueError, np.fft.fft, [1, 2, 3], 0)
+
+
+class TestFFT1D:
+
+    def test_identity(self):
+        maxlen = 512
+        x = random(maxlen) + 1j * random(maxlen)
+        xr = random(maxlen)
+        for i in range(1, maxlen):
+            assert_allclose(np.fft.ifft(np.fft.fft(x[0:i])), x[0:i],
+                            atol=1e-12)
+            assert_allclose(np.fft.irfft(np.fft.rfft(xr[0:i]), i),
+                            xr[0:i], atol=1e-12)
+
+    @pytest.mark.parametrize("dtype", [np.single, np.double, np.longdouble])
+    def test_identity_long_short(self, dtype):
+        # Test with explicitly given number of points, both for n
+        # smaller and for n larger than the input size.
+        maxlen = 16
+        atol = 5 * np.spacing(np.array(1., dtype=dtype))
+        x = random(maxlen).astype(dtype) + 1j * random(maxlen).astype(dtype)
+        xx = np.concatenate([x, np.zeros_like(x)])
+        xr = random(maxlen).astype(dtype)
+        xxr = np.concatenate([xr, np.zeros_like(xr)])
+        for i in range(1, maxlen * 2):
+            check_c = np.fft.ifft(np.fft.fft(x, n=i), n=i)
+            assert check_c.real.dtype == dtype
+            assert_allclose(check_c, xx[0:i], atol=atol, rtol=0)
+            check_r = np.fft.irfft(np.fft.rfft(xr, n=i), n=i)
+            assert check_r.dtype == dtype
+            assert_allclose(check_r, xxr[0:i], atol=atol, rtol=0)
+
+    @pytest.mark.parametrize("dtype", [np.single, np.double, np.longdouble])
+    def test_identity_long_short_reversed(self, dtype):
+        # Also test explicitly given number of points in reversed order.
+        maxlen = 16
+        atol = 5 * np.spacing(np.array(1., dtype=dtype))
+        x = random(maxlen).astype(dtype) + 1j * random(maxlen).astype(dtype)
+        xx = np.concatenate([x, np.zeros_like(x)])
+        for i in range(1, maxlen * 2):
+            check_via_c = np.fft.fft(np.fft.ifft(x, n=i), n=i)
+            assert check_via_c.dtype == x.dtype
+            assert_allclose(check_via_c, xx[0:i], atol=atol, rtol=0)
+            # For irfft, we can neither recover the imaginary part of
+            # the first element, nor the imaginary part of the last
+            # element if npts is even.  So, set to 0 for the comparison.
+            y = x.copy()
+            n = i // 2 + 1
+            y.imag[0] = 0
+            if i % 2 == 0:
+                y.imag[n - 1:] = 0
+            yy = np.concatenate([y, np.zeros_like(y)])
+            check_via_r = np.fft.rfft(np.fft.irfft(x, n=i), n=i)
+            assert check_via_r.dtype == x.dtype
+            assert_allclose(check_via_r, yy[0:n], atol=atol, rtol=0)
+
+    def test_fft(self):
+        x = random(30) + 1j * random(30)
+        assert_allclose(fft1(x), np.fft.fft(x), atol=1e-6)
+        assert_allclose(fft1(x), np.fft.fft(x, norm="backward"), atol=1e-6)
+        assert_allclose(fft1(x) / np.sqrt(30),
+                        np.fft.fft(x, norm="ortho"), atol=1e-6)
+        assert_allclose(fft1(x) / 30.,
+                        np.fft.fft(x, norm="forward"), atol=1e-6)
+
+    @pytest.mark.parametrize("axis", (0, 1))
+    @pytest.mark.parametrize("dtype", (complex, float))
+    @pytest.mark.parametrize("transpose", (True, False))
+    def test_fft_out_argument(self, dtype, transpose, axis):
+        def zeros_like(x):
+            if transpose:
+                return np.zeros_like(x.T).T
+            else:
+                return np.zeros_like(x)
+
+        # tests below only test the out parameter
+        if dtype is complex:
+            y = random((10, 20)) + 1j * random((10, 20))
+            fft, ifft = np.fft.fft, np.fft.ifft
+        else:
+            y = random((10, 20))
+            fft, ifft = np.fft.rfft, np.fft.irfft
+
+        expected = fft(y, axis=axis)
+        out = zeros_like(expected)
+        result = fft(y, out=out, axis=axis)
+        assert result is out
+        assert_array_equal(result, expected)
+
+        expected2 = ifft(expected, axis=axis)
+        out2 = out if dtype is complex else zeros_like(expected2)
+        result2 = ifft(out, out=out2, axis=axis)
+        assert result2 is out2
+        assert_array_equal(result2, expected2)
+
+    @pytest.mark.parametrize("axis", [0, 1])
+    def test_fft_inplace_out(self, axis):
+        # Test some weirder in-place combinations
+        y = random((20, 20)) + 1j * random((20, 20))
+        # Fully in-place.
+        y1 = y.copy()
+        expected1 = np.fft.fft(y1, axis=axis)
+        result1 = np.fft.fft(y1, axis=axis, out=y1)
+        assert result1 is y1
+        assert_array_equal(result1, expected1)
+        # In-place of part of the array; rest should be unchanged.
+        y2 = y.copy()
+        out2 = y2[:10] if axis == 0 else y2[:, :10]
+        expected2 = np.fft.fft(y2, n=10, axis=axis)
+        result2 = np.fft.fft(y2, n=10, axis=axis, out=out2)
+        assert result2 is out2
+        assert_array_equal(result2, expected2)
+        if axis == 0:
+            assert_array_equal(y2[10:], y[10:])
+        else:
+            assert_array_equal(y2[:, 10:], y[:, 10:])
+        # In-place of another part of the array.
+        y3 = y.copy()
+        y3_sel = y3[5:] if axis == 0 else y3[:, 5:]
+        out3 = y3[5:15] if axis == 0 else y3[:, 5:15]
+        expected3 = np.fft.fft(y3_sel, n=10, axis=axis)
+        result3 = np.fft.fft(y3_sel, n=10, axis=axis, out=out3)
+        assert result3 is out3
+        assert_array_equal(result3, expected3)
+        if axis == 0:
+            assert_array_equal(y3[:5], y[:5])
+            assert_array_equal(y3[15:], y[15:])
+        else:
+            assert_array_equal(y3[:, :5], y[:, :5])
+            assert_array_equal(y3[:, 15:], y[:, 15:])
+        # In-place with n > nin; rest should be unchanged.
+        y4 = y.copy()
+        y4_sel = y4[:10] if axis == 0 else y4[:, :10]
+        out4 = y4[:15] if axis == 0 else y4[:, :15]
+        expected4 = np.fft.fft(y4_sel, n=15, axis=axis)
+        result4 = np.fft.fft(y4_sel, n=15, axis=axis, out=out4)
+        assert result4 is out4
+        assert_array_equal(result4, expected4)
+        if axis == 0:
+            assert_array_equal(y4[15:], y[15:])
+        else:
+            assert_array_equal(y4[:, 15:], y[:, 15:])
+        # Overwrite in a transpose.
+        y5 = y.copy()
+        out5 = y5.T
+        result5 = np.fft.fft(y5, axis=axis, out=out5)
+        assert result5 is out5
+        assert_array_equal(result5, expected1)
+        # Reverse strides.
+        y6 = y.copy()
+        out6 = y6[::-1] if axis == 0 else y6[:, ::-1]
+        result6 = np.fft.fft(y6, axis=axis, out=out6)
+        assert result6 is out6
+        assert_array_equal(result6, expected1)
+
+    def test_fft_bad_out(self):
+        x = np.arange(30.)
+        with pytest.raises(TypeError, match="must be of ArrayType"):
+            np.fft.fft(x, out="")
+        with pytest.raises(ValueError, match="has wrong shape"):
+            np.fft.fft(x, out=np.zeros_like(x).reshape(5, -1))
+        with pytest.raises(TypeError, match="Cannot cast"):
+            np.fft.fft(x, out=np.zeros_like(x, dtype=float))
+
+    @pytest.mark.parametrize('norm', (None, 'backward', 'ortho', 'forward'))
+    def test_ifft(self, norm):
+        x = random(30) + 1j * random(30)
+        assert_allclose(
+            x, np.fft.ifft(np.fft.fft(x, norm=norm), norm=norm),
+            atol=1e-6)
+        # Ensure we get the correct error message
+        with pytest.raises(ValueError,
+                           match='Invalid number of FFT data points'):
+            np.fft.ifft([], norm=norm)
+
+    def test_fft2(self):
+        x = random((30, 20)) + 1j * random((30, 20))
+        assert_allclose(np.fft.fft(np.fft.fft(x, axis=1), axis=0),
+                        np.fft.fft2(x), atol=1e-6)
+        assert_allclose(np.fft.fft2(x),
+                        np.fft.fft2(x, norm="backward"), atol=1e-6)
+        assert_allclose(np.fft.fft2(x) / np.sqrt(30 * 20),
+                        np.fft.fft2(x, norm="ortho"), atol=1e-6)
+        assert_allclose(np.fft.fft2(x) / (30. * 20.),
+                        np.fft.fft2(x, norm="forward"), atol=1e-6)
+
+    def test_ifft2(self):
+        x = random((30, 20)) + 1j * random((30, 20))
+        assert_allclose(np.fft.ifft(np.fft.ifft(x, axis=1), axis=0),
+                        np.fft.ifft2(x), atol=1e-6)
+        assert_allclose(np.fft.ifft2(x),
+                        np.fft.ifft2(x, norm="backward"), atol=1e-6)
+        assert_allclose(np.fft.ifft2(x) * np.sqrt(30 * 20),
+                        np.fft.ifft2(x, norm="ortho"), atol=1e-6)
+        assert_allclose(np.fft.ifft2(x) * (30. * 20.),
+                        np.fft.ifft2(x, norm="forward"), atol=1e-6)
+
+    def test_fftn(self):
+        x = random((30, 20, 10)) + 1j * random((30, 20, 10))
+        assert_allclose(
+            np.fft.fft(np.fft.fft(np.fft.fft(x, axis=2), axis=1), axis=0),
+            np.fft.fftn(x), atol=1e-6)
+        assert_allclose(np.fft.fftn(x),
+                        np.fft.fftn(x, norm="backward"), atol=1e-6)
+        assert_allclose(np.fft.fftn(x) / np.sqrt(30 * 20 * 10),
+                        np.fft.fftn(x, norm="ortho"), atol=1e-6)
+        assert_allclose(np.fft.fftn(x) / (30. * 20. * 10.),
+                        np.fft.fftn(x, norm="forward"), atol=1e-6)
+
+    def test_ifftn(self):
+        x = random((30, 20, 10)) + 1j * random((30, 20, 10))
+        assert_allclose(
+            np.fft.ifft(np.fft.ifft(np.fft.ifft(x, axis=2), axis=1), axis=0),
+            np.fft.ifftn(x), atol=1e-6)
+        assert_allclose(np.fft.ifftn(x),
+                        np.fft.ifftn(x, norm="backward"), atol=1e-6)
+        assert_allclose(np.fft.ifftn(x) * np.sqrt(30 * 20 * 10),
+                        np.fft.ifftn(x, norm="ortho"), atol=1e-6)
+        assert_allclose(np.fft.ifftn(x) * (30. * 20. * 10.),
+                        np.fft.ifftn(x, norm="forward"), atol=1e-6)
+
+    def test_rfft(self):
+        x = random(30)
+        for n in [x.size, 2 * x.size]:
+            for norm in [None, 'backward', 'ortho', 'forward']:
+                assert_allclose(
+                    np.fft.fft(x, n=n, norm=norm)[:(n // 2 + 1)],
+                    np.fft.rfft(x, n=n, norm=norm), atol=1e-6)
+            assert_allclose(
+                np.fft.rfft(x, n=n),
+                np.fft.rfft(x, n=n, norm="backward"), atol=1e-6)
+            assert_allclose(
+                np.fft.rfft(x, n=n) / np.sqrt(n),
+                np.fft.rfft(x, n=n, norm="ortho"), atol=1e-6)
+            assert_allclose(
+                np.fft.rfft(x, n=n) / n,
+                np.fft.rfft(x, n=n, norm="forward"), atol=1e-6)
+
+    def test_rfft_even(self):
+        x = np.arange(8)
+        n = 4
+        y = np.fft.rfft(x, n)
+        assert_allclose(y, np.fft.fft(x[:n])[:n // 2 + 1], rtol=1e-14)
+
+    def test_rfft_odd(self):
+        x = np.array([1, 0, 2, 3, -3])
+        y = np.fft.rfft(x)
+        assert_allclose(y, np.fft.fft(x)[:3], rtol=1e-14)
+
+    def test_irfft(self):
+        x = random(30)
+        assert_allclose(x, np.fft.irfft(np.fft.rfft(x)), atol=1e-6)
+        assert_allclose(x, np.fft.irfft(np.fft.rfft(x, norm="backward"),
+                        norm="backward"), atol=1e-6)
+        assert_allclose(x, np.fft.irfft(np.fft.rfft(x, norm="ortho"),
+                        norm="ortho"), atol=1e-6)
+        assert_allclose(x, np.fft.irfft(np.fft.rfft(x, norm="forward"),
+                        norm="forward"), atol=1e-6)
+
+    def test_rfft2(self):
+        x = random((30, 20))
+        assert_allclose(np.fft.fft2(x)[:, :11], np.fft.rfft2(x), atol=1e-6)
+        assert_allclose(np.fft.rfft2(x),
+                        np.fft.rfft2(x, norm="backward"), atol=1e-6)
+        assert_allclose(np.fft.rfft2(x) / np.sqrt(30 * 20),
+                        np.fft.rfft2(x, norm="ortho"), atol=1e-6)
+        assert_allclose(np.fft.rfft2(x) / (30. * 20.),
+                        np.fft.rfft2(x, norm="forward"), atol=1e-6)
+
+    def test_irfft2(self):
+        x = random((30, 20))
+        assert_allclose(x, np.fft.irfft2(np.fft.rfft2(x)), atol=1e-6)
+        assert_allclose(x, np.fft.irfft2(np.fft.rfft2(x, norm="backward"),
+                        norm="backward"), atol=1e-6)
+        assert_allclose(x, np.fft.irfft2(np.fft.rfft2(x, norm="ortho"),
+                        norm="ortho"), atol=1e-6)
+        assert_allclose(x, np.fft.irfft2(np.fft.rfft2(x, norm="forward"),
+                        norm="forward"), atol=1e-6)
+
+    def test_rfftn(self):
+        x = random((30, 20, 10))
+        assert_allclose(np.fft.fftn(x)[:, :, :6], np.fft.rfftn(x), atol=1e-6)
+        assert_allclose(np.fft.rfftn(x),
+                        np.fft.rfftn(x, norm="backward"), atol=1e-6)
+        assert_allclose(np.fft.rfftn(x) / np.sqrt(30 * 20 * 10),
+                        np.fft.rfftn(x, norm="ortho"), atol=1e-6)
+        assert_allclose(np.fft.rfftn(x) / (30. * 20. * 10.),
+                        np.fft.rfftn(x, norm="forward"), atol=1e-6)
+        # Regression test for gh-27159
+        x = np.ones((2, 3))
+        result = np.fft.rfftn(x, axes=(0, 0, 1), s=(10, 20, 40))
+        assert result.shape == (10, 21)
+        expected = np.fft.fft(np.fft.fft(np.fft.rfft(x, axis=1, n=40),
+                            axis=0, n=20), axis=0, n=10)
+        assert expected.shape == (10, 21)
+        assert_allclose(result, expected, atol=1e-6)
+
+    def test_irfftn(self):
+        x = random((30, 20, 10))
+        assert_allclose(x, np.fft.irfftn(np.fft.rfftn(x)), atol=1e-6)
+        assert_allclose(x, np.fft.irfftn(np.fft.rfftn(x, norm="backward"),
+                        norm="backward"), atol=1e-6)
+        assert_allclose(x, np.fft.irfftn(np.fft.rfftn(x, norm="ortho"),
+                        norm="ortho"), atol=1e-6)
+        assert_allclose(x, np.fft.irfftn(np.fft.rfftn(x, norm="forward"),
+                        norm="forward"), atol=1e-6)
+
+    def test_hfft(self):
+        x = random(14) + 1j * random(14)
+        x_herm = np.concatenate((random(1), x, random(1)))
+        x = np.concatenate((x_herm, x[::-1].conj()))
+        assert_allclose(np.fft.fft(x), np.fft.hfft(x_herm), atol=1e-6)
+        assert_allclose(np.fft.hfft(x_herm),
+                        np.fft.hfft(x_herm, norm="backward"), atol=1e-6)
+        assert_allclose(np.fft.hfft(x_herm) / np.sqrt(30),
+                        np.fft.hfft(x_herm, norm="ortho"), atol=1e-6)
+        assert_allclose(np.fft.hfft(x_herm) / 30.,
+                        np.fft.hfft(x_herm, norm="forward"), atol=1e-6)
+
+    def test_ihfft(self):
+        x = random(14) + 1j * random(14)
+        x_herm = np.concatenate((random(1), x, random(1)))
+        x = np.concatenate((x_herm, x[::-1].conj()))
+        assert_allclose(x_herm, np.fft.ihfft(np.fft.hfft(x_herm)), atol=1e-6)
+        assert_allclose(x_herm, np.fft.ihfft(np.fft.hfft(x_herm,
+                        norm="backward"), norm="backward"), atol=1e-6)
+        assert_allclose(x_herm, np.fft.ihfft(np.fft.hfft(x_herm,
+                        norm="ortho"), norm="ortho"), atol=1e-6)
+        assert_allclose(x_herm, np.fft.ihfft(np.fft.hfft(x_herm,
+                        norm="forward"), norm="forward"), atol=1e-6)
+
+    @pytest.mark.parametrize("op", [np.fft.fftn, np.fft.ifftn,
+                                    np.fft.rfftn, np.fft.irfftn])
+    def test_axes(self, op):
+        x = random((30, 20, 10))
+        axes = [(0, 1, 2), (0, 2, 1), (1, 0, 2), (1, 2, 0), (2, 0, 1), (2, 1, 0)]
+        for a in axes:
+            op_tr = op(np.transpose(x, a))
+            tr_op = np.transpose(op(x, axes=a), a)
+            assert_allclose(op_tr, tr_op, atol=1e-6)
+
+    @pytest.mark.parametrize("op", [np.fft.fftn, np.fft.ifftn,
+                                    np.fft.fft2, np.fft.ifft2])
+    def test_s_negative_1(self, op):
+        x = np.arange(100).reshape(10, 10)
+        # should use the whole input array along the first axis
+        assert op(x, s=(-1, 5), axes=(0, 1)).shape == (10, 5)
+
+    @pytest.mark.parametrize("op", [np.fft.fftn, np.fft.ifftn,
+                                    np.fft.rfftn, np.fft.irfftn])
+    def test_s_axes_none(self, op):
+        x = np.arange(100).reshape(10, 10)
+        with pytest.warns(match='`axes` should not be `None` if `s`'):
+            op(x, s=(-1, 5))
+
+    @pytest.mark.parametrize("op", [np.fft.fft2, np.fft.ifft2])
+    def test_s_axes_none_2D(self, op):
+        x = np.arange(100).reshape(10, 10)
+        with pytest.warns(match='`axes` should not be `None` if `s`'):
+            op(x, s=(-1, 5), axes=None)
+
+    @pytest.mark.parametrize("op", [np.fft.fftn, np.fft.ifftn,
+                                    np.fft.rfftn, np.fft.irfftn,
+                                    np.fft.fft2, np.fft.ifft2])
+    def test_s_contains_none(self, op):
+        x = random((30, 20, 10))
+        with pytest.warns(match='array containing `None` values to `s`'):
+            op(x, s=(10, None, 10), axes=(0, 1, 2))
+
+    def test_all_1d_norm_preserving(self):
+        # verify that round-trip transforms are norm-preserving
+        x = random(30)
+        x_norm = np.linalg.norm(x)
+        n = x.size * 2
+        func_pairs = [(np.fft.fft, np.fft.ifft),
+                      (np.fft.rfft, np.fft.irfft),
+                      # hfft: order so the first function takes x.size samples
+                      #       (necessary for comparison to x_norm above)
+                      (np.fft.ihfft, np.fft.hfft),
+                      ]
+        for forw, back in func_pairs:
+            for n in [x.size, 2 * x.size]:
+                for norm in [None, 'backward', 'ortho', 'forward']:
+                    tmp = forw(x, n=n, norm=norm)
+                    tmp = back(tmp, n=n, norm=norm)
+                    assert_allclose(x_norm,
+                                    np.linalg.norm(tmp), atol=1e-6)
+
+    @pytest.mark.parametrize("axes", [(0, 1), (0, 2), None])
+    @pytest.mark.parametrize("dtype", (complex, float))
+    @pytest.mark.parametrize("transpose", (True, False))
+    def test_fftn_out_argument(self, dtype, transpose, axes):
+        def zeros_like(x):
+            if transpose:
+                return np.zeros_like(x.T).T
+            else:
+                return np.zeros_like(x)
+
+        # tests below only test the out parameter
+        if dtype is complex:
+            x = random((10, 5, 6)) + 1j * random((10, 5, 6))
+            fft, ifft = np.fft.fftn, np.fft.ifftn
+        else:
+            x = random((10, 5, 6))
+            fft, ifft = np.fft.rfftn, np.fft.irfftn
+
+        expected = fft(x, axes=axes)
+        out = zeros_like(expected)
+        result = fft(x, out=out, axes=axes)
+        assert result is out
+        assert_array_equal(result, expected)
+
+        expected2 = ifft(expected, axes=axes)
+        out2 = out if dtype is complex else zeros_like(expected2)
+        result2 = ifft(out, out=out2, axes=axes)
+        assert result2 is out2
+        assert_array_equal(result2, expected2)
+
+    @pytest.mark.parametrize("fft", [np.fft.fftn, np.fft.ifftn, np.fft.rfftn])
+    def test_fftn_out_and_s_interaction(self, fft):
+        # With s, shape varies, so generally one cannot pass in out.
+        if fft is np.fft.rfftn:
+            x = random((10, 5, 6))
+        else:
+            x = random((10, 5, 6)) + 1j * random((10, 5, 6))
+        with pytest.raises(ValueError, match="has wrong shape"):
+            fft(x, out=np.zeros_like(x), s=(3, 3, 3), axes=(0, 1, 2))
+        # Except on the first axis done (which is the last of axes).
+        s = (10, 5, 5)
+        expected = fft(x, s=s, axes=(0, 1, 2))
+        out = np.zeros_like(expected)
+        result = fft(x, s=s, axes=(0, 1, 2), out=out)
+        assert result is out
+        assert_array_equal(result, expected)
+
+    @pytest.mark.parametrize("s", [(9, 5, 5), (3, 3, 3)])
+    def test_irfftn_out_and_s_interaction(self, s):
+        # Since for irfftn, the output is real and thus cannot be used for
+        # intermediate steps, it should always work.
+        x = random((9, 5, 6, 2)) + 1j * random((9, 5, 6, 2))
+        expected = np.fft.irfftn(x, s=s, axes=(0, 1, 2))
+        out = np.zeros_like(expected)
+        result = np.fft.irfftn(x, s=s, axes=(0, 1, 2), out=out)
+        assert result is out
+        assert_array_equal(result, expected)
+
+
+@pytest.mark.parametrize(
+        "dtype",
+        [np.float32, np.float64, np.complex64, np.complex128])
+@pytest.mark.parametrize("order", ["F", 'non-contiguous'])
+@pytest.mark.parametrize(
+        "fft",
+        [np.fft.fft, np.fft.fft2, np.fft.fftn,
+         np.fft.ifft, np.fft.ifft2, np.fft.ifftn])
+def test_fft_with_order(dtype, order, fft):
+    # Check that FFT/IFFT produces identical results for C, Fortran and
+    # non contiguous arrays
+    rng = np.random.RandomState(42)
+    X = rng.rand(8, 7, 13).astype(dtype, copy=False)
+    # See discussion in pull/14178
+    _tol = 8.0 * np.sqrt(np.log2(X.size)) * np.finfo(X.dtype).eps
+    if order == 'F':
+        Y = np.asfortranarray(X)
+    else:
+        # Make a non contiguous array
+        Y = X[::-1]
+        X = np.ascontiguousarray(X[::-1])
+
+    if fft.__name__.endswith('fft'):
+        for axis in range(3):
+            X_res = fft(X, axis=axis)
+            Y_res = fft(Y, axis=axis)
+            assert_allclose(X_res, Y_res, atol=_tol, rtol=_tol)
+    elif fft.__name__.endswith(('fft2', 'fftn')):
+        axes = [(0, 1), (1, 2), (0, 2)]
+        if fft.__name__.endswith('fftn'):
+            axes.extend([(0,), (1,), (2,), None])
+        for ax in axes:
+            X_res = fft(X, axes=ax)
+            Y_res = fft(Y, axes=ax)
+            assert_allclose(X_res, Y_res, atol=_tol, rtol=_tol)
+    else:
+        raise ValueError
+
+
+@pytest.mark.parametrize("order", ["F", "C"])
+@pytest.mark.parametrize("n", [None, 7, 12])
+def test_fft_output_order(order, n):
+    rng = np.random.RandomState(42)
+    x = rng.rand(10)
+    x = np.asarray(x, dtype=np.complex64, order=order)
+    res = np.fft.fft(x, n=n)
+    assert res.flags.c_contiguous == x.flags.c_contiguous
+    assert res.flags.f_contiguous == x.flags.f_contiguous
+
+@pytest.mark.skipif(IS_WASM, reason="Cannot start thread")
+class TestFFTThreadSafe:
+    threads = 16
+    input_shape = (800, 200)
+
+    def _test_mtsame(self, func, *args):
+        def worker(args, q):
+            q.put(func(*args))
+
+        q = queue.Queue()
+        expected = func(*args)
+
+        # Spin off a bunch of threads to call the same function simultaneously
+        t = [threading.Thread(target=worker, args=(args, q))
+             for i in range(self.threads)]
+        [x.start() for x in t]
+
+        [x.join() for x in t]
+        # Make sure all threads returned the correct value
+        for i in range(self.threads):
+            assert_array_equal(q.get(timeout=5), expected,
+                'Function returned wrong value in multithreaded context')
+
+    def test_fft(self):
+        a = np.ones(self.input_shape) * 1 + 0j
+        self._test_mtsame(np.fft.fft, a)
+
+    def test_ifft(self):
+        a = np.ones(self.input_shape) * 1 + 0j
+        self._test_mtsame(np.fft.ifft, a)
+
+    def test_rfft(self):
+        a = np.ones(self.input_shape)
+        self._test_mtsame(np.fft.rfft, a)
+
+    def test_irfft(self):
+        a = np.ones(self.input_shape) * 1 + 0j
+        self._test_mtsame(np.fft.irfft, a)
+
+
+def test_irfft_with_n_1_regression():
+    # Regression test for gh-25661
+    x = np.arange(10)
+    np.fft.irfft(x, n=1)
+    np.fft.hfft(x, n=1)
+    np.fft.irfft(np.array([0], complex), n=10)
+
+
+def test_irfft_with_n_large_regression():
+    # Regression test for gh-25679
+    x = np.arange(5) * (1 + 1j)
+    result = np.fft.hfft(x, n=10)
+    expected = np.array([20., 9.91628173, -11.8819096, 7.1048486,
+                         -6.62459848, 4., -3.37540152, -0.16057669,
+                         1.8819096, -20.86055364])
+    assert_allclose(result, expected)
+
+
+@pytest.mark.parametrize("fft", [
+    np.fft.fft, np.fft.ifft, np.fft.rfft, np.fft.irfft
+])
+@pytest.mark.parametrize("data", [
+    np.array([False, True, False]),
+    np.arange(10, dtype=np.uint8),
+    np.arange(5, dtype=np.int16),
+])
+def test_fft_with_integer_or_bool_input(data, fft):
+    # Regression test for gh-25819
+    result = fft(data)
+    float_data = data.astype(np.result_type(data, 1.))
+    expected = fft(float_data)
+    assert_array_equal(result, expected)