BUG: correct irfft with n=1 on larger input #25668
Conversation
faea9fb to
502a7b9
Compare
2800814 to
0b7e4d6
Compare
|
I can confirm this fixes the segfault. |
|
@seberg - would you be able to review the correction to #25536. Note that I also decided to include the missing release note that you asked for in #25536 (comment) |
|
p.s. I skipped ci on the release note addition, but you can still see that the tests passed for the actual code commit. |
There was a problem hiding this comment.
This looks right to me, but honestly I find few things more confusing than these FFT lengths (considering how simple it is in some sense).
@ngoldbaum maybe you (or someone else) also can have a quick look?
numpy/fft/_pocketfft_umath.c
Outdated
| if (i_last > 0) { | ||
| op_double[i_last * 2 - 1] = last[0]; /* real of the last point */ | ||
| } | ||
| if (npts % 2 == 1) { /* For odd n, we still imag of the last point */ |
There was a problem hiding this comment.
Comment is mangled slightly, but OK.
There was a problem hiding this comment.
I'll try to get it. And might as well add a comment about what FFTPACK format actually is, since it helps to understand what is done...
There was a problem hiding this comment.
I think of it in terms of the real transform taking advantage of the Hermitian symmetry of the spectrum to only return the lower half. The question is then how to reconstruct the full spectrum, which depends on whether the center of symmetry is at the highest frequency, or half a sample to the right. This same sort of question involving reflection symmetry turns up in many places.
EDIT: The highest frequency is required to be real in the first case, is full complex in the second, so one can also determine the symmetry by counting degrees of freedom, which is what is happening here.
b01cf67 to
1350f6f
Compare
|
Yikes, that took me long to realize I was just making things far too complicated. But now it is finally OK. |
numpy/fft/_pocketfft_umath.c
Outdated
| (char *)&op_double[1], sizeof(npy_cdouble), npts / 2 - 1, | ||
| /* | ||
| * Copy complex input to buffer in FFTpack order: | ||
| * R0,R1,I1,...Rn-1,In-1,Rn[,In] (last for npts odd) |
There was a problem hiding this comment.
I got tripped up reading this that there's no I0, maybe a comment explicitly saying that the first entry is defined to be real?
There was a problem hiding this comment.
The data is often packed a bit further in the even case, R0, RN, R1, I1 ... . Don't know what fftpack did, but the form NumPy uses differed from the form SciPy returned back in the day. Having all complex numbers in the proper order and alignment is much less confusing.
There was a problem hiding this comment.
Yes, the order @charris mentions is that used in Numerical Recipes, which I agree is nicer. But that's not what pocketfft and fftpack do...
Anyway, yes, happy to add a comment! Since this code is for an irfft, by the assumption that the input values to the forward fft were real, the zero frequency element, which is just the sum of all the input elements, must be real too. For even number of points, by symmetry, the Nyquist frequency is also real (and of course just by counting the number of degrees of freedom one reaches the same conclusion).
1350f6f to
aa6e9b8
Compare
[skip ci]
|
OK, I added a bit of documentation (for all routines). @ngoldbaum - looks good like this? |
|
It looks like you checked in an old submodule state, that's causing the CI failure. |
|
LGTM now, feel free to self-merge once the CI is green! |
For n=1, the previously led to segfaults. Regression tests added.
aa6e9b8 to
f0322c8
Compare
|
OK, thanks! Indeed, somehow the submodules got messed up. Will merge once CI passes! |
Fixes #25661
Fixes #25679