Description
While coding an iterative FFT function, I came upon a strange bug that depends on the arguments to the function in a complicated way. I've managed to reduce the problematic code to the following:
using FFTW
unitize(x::Complex) = x==0 ? 0 : x/abs(x)
function problematic(ys::Tuple,phis::Tuple,nit)
n = length(ys)
x = sum(ifft(ys[i]) for i=1:n)/n
for j=1:nit
thetas = ((unitize.(fft(x.*phi)) for phi in phis)...,)
x = sum(ifft(ys[k].*thetas[k]) for k=1:n)/n
end
return x
end
x = [exp(-j^2) for j=-10:.1:10]
quad = exp.(im*[-10:.1:10;].^2)
y0 = abs.(fft(x))
y1 = abs.(fft(x.*quad))
x1 = problematic((y0,y1),(ones(length(x)),quad),4); # This works
x1 = problematic((y0,y1),(ones(length(x)),quad),5); # But this does not work
The various parameters and the auxiliary function unitize are, unfortunately, necessary for this demonstration, as slight tweaks to the function eliminates the error. The final parameter to problematic determines the number of iterations, and changing it slightly affects whether or not an error is thrown. With slightly different inputs y0, y1, quad, etc. the exact number of iterations at which the error shows up will change also. Under some conditions, I also find that the number of iterations required for an error differs between the REPL and a Jupyter notebook.
I also find it strange that the error mentions plan_bfft, as I do not invoke that at all. Is the compiler automatically planning an FFT when it sees many FFTs in a loop?
I'm using Julia v1.7.2 and FFTW v1.4.6.