@@ -209,6 +209,67 @@ def _init_nd_shape_and_axes(x, shape, axes):
209209
210210
211211def _iter_complementary (x , axes , func , kwargs , result ):
212+ """
213+ Apply FFT function by iterating over complementary axes.
214+
215+ This function applies an FFT operation to slices of the input array
216+ by iterating over all axes that are NOT in the `axes` parameter
217+ (the complementary axes). For each position in the complementary axes,
218+ it applies the FFT function to a slice along the specified axes.
219+
220+ Parameters
221+ ----------
222+ x : ndarray
223+ Input array.
224+ axes : int, sequence of ints, or None
225+ Axes along which to perform the FFT operation. The function iterates
226+ over the complementary axes (axes not in this parameter). If ``None``,
227+ performs direct N-D FFT without iteration.
228+ Default: None
229+ func : callable
230+ FFT function to apply to each slice. Should accept array input and
231+ return transformed output.
232+ kwargs : dict
233+ Additional keyword arguments to pass to `func`.
234+ result : ndarray
235+ Pre-allocated output array where results are stored.
236+
237+ Returns
238+ -------
239+ ndarray
240+ The transformed array (same as `result`).
241+
242+ Notes
243+ -----
244+ For complex input, the function uses in-place operations with the `out`
245+ parameter passed for better performance. For real input, `np.copyto` is
246+ used instead to avoid element ordering issues that can occur with the
247+ `out` parameter in certain FFT operations.
248+
249+ Examples
250+ --------
251+ Consider an input array with shape (3, 4, 5) and performing FFT
252+ along axis 2 only:
253+
254+ >>> x = np.random.random((3, 4, 5))
255+ >>> result = np.empty((3, 4, 5), dtype=np.complex128)
256+ >>> _iter_complementary(
257+ ... x, axes=(2,), func=_direct_fftnd,
258+ ... kwargs={'direction': 1, 'fsc': 1.0}, result=result
259+ ... )
260+
261+ The function will iterate over axes 0 and 1 (complementary axes)
262+ and apply `_direct_fftnd` to each 1-D slice along axis 2:
263+
264+ - Iteration 0: func(x[0, 0, :]) -> result[0, 0, :]
265+ - Iteration 1: func(x[0, 1, :]) -> result[0, 1, :]
266+ - ...
267+ - Iteration 11: func(x[2, 3, :]) -> result[2, 3, :]
268+
269+ Total: 3 * 4 = 12 FFT operations on arrays of shape (5,).
270+
271+ """
272+
212273 if axes is None :
213274 # s and axes are None, direct N-D FFT
214275 return func (x , ** kwargs , out = result )
@@ -233,8 +294,10 @@ def _iter_complementary(x, axes, func, kwargs, result):
233294 m_ind = _flat_to_multi (ind , sub_shape )
234295 for k1 , k2 in zip (dual_ind , m_ind ):
235296 sl [k1 ] = k2
297+ tsl = tuple (sl )
298+
236299 if np .issubdtype (x .dtype , np .complexfloating ):
237- func (x [tuple ( sl ) ], ** kwargs , out = result [tuple ( sl ) ])
300+ func (x [tsl ], ** kwargs , out = result [tsl ])
238301 else :
239302 # For c2c FFT, if the input is real, half of the output is the
240303 # complex conjugate of the other half. Instead of upcasting the
@@ -247,7 +310,7 @@ def _iter_complementary(x, axes, func, kwargs, result):
247310 # array appeared in the second half of the NumPy output array,
248311 # while the equivalent element in the NumPy array was the conjugate
249312 # of the mkl_fft output array.
250- np .copyto (result [tuple ( sl ) ], func (x [tuple ( sl ) ], ** kwargs ))
313+ np .copyto (result [tsl ], func (x [tsl ], ** kwargs ))
251314
252315 return result
253316
@@ -260,7 +323,49 @@ def _iter_fftnd(
260323 direction = + 1 ,
261324 scale_function = lambda ind : 1.0 ,
262325):
263- a = np .asarray (a )
326+ """
327+ Perform N-D FFT as a series of 1-D FFTs along specified axes.
328+
329+ This function implements N-D FFT by applying 1-D FFT iteratively along each
330+ axis. The axes are processed in reverse order to end with the first axis
331+ given.
332+
333+ Parameters
334+ ----------
335+ a : ndarray
336+ Input array.
337+ s : sequence of ints, optional
338+ Shape of the FFT output along each axis in `axes`. If not provided, the
339+ shape is inferred from the input array.
340+ Default: ``None``
341+ axes : sequence of ints, optional
342+ Axes along which to compute the FFT. If not provided, all axes are used.
343+ Default: ``None``
344+ out : ndarray, optional
345+ Output array to store the result. Used for in-place operations when
346+ possible.
347+ Default: ``None``
348+ direction : int, optional
349+ FFT direction: ``+1`` for forward FFT, ``-1`` for inverse FFT.
350+ Default: ``+1``
351+ scale_function : callable, optional
352+ Function that takes iteration index and returns the scaling factor for
353+ that step. Used to apply normalization at specific iteration steps.
354+ Default: ``lambda ind: 1.0``
355+
356+ Returns
357+ -------
358+ ndarray
359+ The transformed array.
360+
361+ Notes
362+ -----
363+ The function optimizes memory usage by performing in-place calculations
364+ when possible. In-place operations are used everywhere except when the
365+ array size changes after the first FFT along an axis.
366+
367+ """
368+
264369 s , axes = _init_nd_shape_and_axes (a , s , axes )
265370
266371 # Combine the two, but in reverse, to end with the first axis given.
@@ -412,8 +517,17 @@ def _c2c_fftnd_impl(
412517 out = out ,
413518 )
414519 else :
520+ x = np .asarray (x )
521+
522+ # Fast path: FFT over no axes is a complete identity operation.
523+ # Returns the input unchanged (same object, no copy), preserving
524+ # dtype and avoiding any FFT computation. The out parameter is
525+ # ignored to match NumPy behavior.
526+ _ , xa = _cook_nd_args (x , s , axes )
527+ if len (xa ) == 0 :
528+ return x
529+
415530 if _complementary and x .dtype in valid_dtypes :
416- x = np .asarray (x )
417531 if out is None :
418532 res = np .empty_like (x , dtype = _output_dtype (x .dtype ))
419533 else :
0 commit comments