fix an issue for real inputs of irfft

Author: Vahid Tavanashad

CHANGELOG.md

@@ -13,6 +13,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Changed
 * NumPy interface `mkl_fft.interfaces.numpy_fft` is aligned with numpy-2.* [gh-139](https://github.com/IntelPython/mkl_fft/pull/139), [gh-157](https://github.com/IntelPython/mkl_fft/pull/157)
 
+### Fixed
+* Fixed an issue for calling `mkl_fft.interfaces.numpy.fftn` with an empty axes [gh-139](https://github.com/IntelPython/mkl_fft/pull/139)
+* Fixed an issue for calling `mkl_fft.interfaces.numpy.fftn` with a zero-size array [gh-139](https://github.com/IntelPython/mkl_fft/pull/139)
+* Fixed inconsistency of input and output arrays dtype for `irfft` function [gh-180](https://github.com/IntelPython/mkl_fft/pull/180)
+
 ## [1.3.14] (04/10/2025)
 
 resolves gh-152 by adding an explicit `mkl-service` dependency to `mkl-fft` when building the wheel

mkl_fft/_pydfti.pyx

@@ -89,7 +89,7 @@ def _tls_dfti_cache_capsule():
 cdef extern from "Python.h":
     ctypedef int size_t
 
-    long PyInt_AsLong(object ob)
+    long PyLong_AsLong(object ob)
     int PyObject_HasAttrString(object, char*)
 
 
@@ -262,7 +262,7 @@ cdef cnp.ndarray _process_arguments(
     xnd[0] = cnp.PyArray_NDIM(x_arr)  # tensor-rank of the array
 
     err = 0
-    axis_[0] = PyInt_AsLong(axis)
+    axis_[0] = PyLong_AsLong(axis)
     if (axis_[0] == -1 and PyErr_Occurred()):
         PyErr_Clear()
         err = 1
@@ -278,7 +278,7 @@ cdef cnp.ndarray _process_arguments(
         n_[0] = x_arr.shape[axis_[0]]
     else:
         try:
-            n_[0] = PyInt_AsLong(n)
+            n_[0] = PyLong_AsLong(n)
         except:
             err = 1
 
@@ -334,7 +334,7 @@ cdef int _is_integral(object num):
     if num is None:
         return 0
     try:
-        n = PyInt_AsLong(num)
+        n = PyLong_AsLong(num)
         _integral = 1 if n > 0 else 0
     except:
         _integral = 0
@@ -665,13 +665,12 @@ def _r2c_fft1d_impl(
         return f_arr
 
 
-# this routine is functionally equivalent to numpy.fft.irfft
 def _c2r_fft1d_impl(
     x, n=None, axis=-1, overwrite_x=False, double fsc=1.0, out=None
 ):
     """
-    Uses MKL to perform 1D FFT on the real input array x along the given axis,
-    producing complex output, but giving only half of the harmonics.
+    Uses MKL to perform 1D FFT on the real/complex input array x along the
+    given axis, producing real output.
 
     cf. numpy.fft.irfft
     """
@@ -704,13 +703,13 @@ def _c2r_fft1d_impl(
     else:
         # we must cast the input and allocate the output,
         # so we cast to complex double and operate in place
-        try:
+        if x_type is cnp.NPY_FLOAT:
             x_arr = <cnp.ndarray> cnp.PyArray_FROM_OTF(
-                x_arr, cnp.NPY_CDOUBLE, cnp.NPY_BEHAVED)
-        except:
-            raise ValueError(
-                "First argument should be a real or "
-                "a complex sequence of single or double precision"
+                x_arr, cnp.NPY_CFLOAT, cnp.NPY_BEHAVED
+            )
+        else:
+            x_arr = <cnp.ndarray> cnp.PyArray_FROM_OTF(
+                x_arr, cnp.NPY_CDOUBLE, cnp.NPY_BEHAVED
             )
         x_type = cnp.PyArray_TYPE(x_arr)
         in_place = 1

mkl_fft/interfaces/_numpy_fft.py

@@ -295,13 +295,11 @@ def hfft(a, n=None, axis=-1, norm=None, out=None):
     """
     norm = _swap_direction(norm)
     x = _downcast_float128_array(a)
-    x = np.array(x, copy=True)
-    np.conjugate(x, out=x)
     fsc = _compute_fwd_scale(norm, n, 2 * (x.shape[axis] - 1))
 
     return _trycall(
         mkl_fft.irfft,
-        (x,),
+        (np.conjugate(x),),
         {"n": n, "axis": axis, "fwd_scale": fsc, "out": out},
     )
 
@@ -317,9 +315,9 @@ def ihfft(a, n=None, axis=-1, norm=None, out=None):
     x = _downcast_float128_array(a)
     fsc = _compute_fwd_scale(norm, n, x.shape[axis])
 
-    output = _trycall(
+    result = _trycall(
         mkl_fft.rfft, (x,), {"n": n, "axis": axis, "fwd_scale": fsc, "out": out}
     )
 
-    np.conjugate(output, out=output)
-    return output
+    np.conjugate(result, out=result)
+    return result

mkl_fft/tests/test_fft1d.py

@@ -457,3 +457,13 @@ def test_irfft_out_strided(axis):
     expected = np.fft.irfft(x, axis=axis, out=out)
 
     assert_allclose(result, expected)
+
+
+@requires_numpy_2
+@pytest.mark.parametrize("dt", ["i4", "f4", "f8", "c8", "c16"])
+def test_irfft_dtype(dt):
+    x = np.array(rnd.random((20, 20)), dtype=dt)
+    result = mkl_fft.irfft(x)
+    expected = np.fft.irfft(x)
+    assert result.dtype == expected.dtype
+    assert_allclose(result, expected, rtol=1e-7, atol=1e-7)