BUG: behavioral change of numpy.fft.rfft with numpy 2.x

[serge-sans-paille]

Describe the issue:

With numpy 2.x:

% python -c 'import numpy as np; print(np.version.version); x = np.arange(8.) ; print(np.fft.rfft(x, 4))' 
2.1.0.dev0+git20240423.ef5f10d
[ 6.+0.j -2.+2.j -2.+4.j]

with numpy 1.26.4:

% python -c 'import numpy as np; print(np.version.version); x = np.arange(8.) ; print(np.fft.rfft(x, 4))'
1.26.4
[ 6.+0.j -2.+2.j -2.+0.j]

that's a notable output difference, maybe introduced by the vectorized form of pocketfft?

Reproduce the code example:

import numpy as np
x = np.arange(8.)
print(np.fft.rfft(x, 4))

Error message:

No response

Python and NumPy Versions:

2.1.0.dev0+git20240423.ef5f10d
3.12.2 (main, Feb 21 2024, 00:00:00) [GCC 13.2.1 20231205 (Red Hat 13.2.1-6)]

Runtime Environment:

No response

Context for the issue:

Difference found while porting pythran to numpy 2.x, as pythran mimics the 'legacy' behavior.
I'm unsure whether the new behavior is intended or not, at least it's not documented in the git log as far as I can tell (?)

[seberg]

@mhvk is maybe another one of those "FFT seems simple but is not?" off-by-one errors. I think the SciPy fft returns the same as NumPy used to here, so it seems plausible (didn't try to understand yet what is going on, maybe it makes sense).

[whym1here]

Is this intended behavior?

[WarrenWeckesser]

@whym1here, no, it is a bug. When n is even, the imaginary part of the last coefficient returned by rfft must be zero.

The equivalent calculation using fft is

In [39]: x = np.arange(8)

In [40]: n = 4

In [41]: np.fft.fft(x[:n])[:n//2 + 1]
Out[41]: array([ 6.+0.j, -2.+2.j, -2.+0.j])

[WarrenWeckesser]

I have a possible fix in #26354.

[pijyoi]

There's an over-copy of the 5th input element when only a shorter length fft was required.

>>> np.fft.rfft(np.array([0,1,2,3]), 4)
array([ 6.+0.j, -2.+2.j, -2.+0.j])

>>> np.fft.rfft(np.array([0,1,2,3,123]), 4)
array([ 6.  +0.j, -2.  +2.j, -2.+123.j])

>>> np.fft.rfft(np.array([0,1,2,3,123])[:4], 4)
array([ 6.+0.j, -2.+2.j, -2.+0.j])

[serge-sans-paille]

Merci / thanks!