Wrong results for certain floating point types and FFT sizes

[mfherbst]

When testing this package against FFTW.jl I noticed certain inconsistencies depending on the floating point type and the array size.

To reproduce

import FFTW
import FourierTransforms

# (1) Size 8
data = randn(Float64, 8) .+ 0im
diff = maximum(abs, FFTW.fft(data) - FourierTransforms.fft(Complex{BigFloat}.(data)))
@assert diff < 1e-14  # Works fine

# (2) Size 9, standard floating point types
data = randn(Float64, 9) .+ 0im
diff = maximum(abs, FFTW.fft(data) - FourierTransforms.fft(Complex{Float32}.(data)))
@assert diff < 1e-5  # Works fine

# (3) Size 9, non-standard floating point types
diff = maximum(abs, FFTW.fft(data) - FourierTransforms.fft(Complex{BigFloat}.(data)))
@assert diff < 1e-14  # diff is about 8

Nine is the first array size which behaves as such. Arrays of sizes 1 to 8 are fine. Using other non-standard floating point types (e.g. DoubleFloats) seem to behave like in case (3) with roughly agreeing errors for what I checked, so it's not just BigFloat. Other problematic array sizes in the range 1 to 30 are 15, 18, 21, 25, 27, 30 ... I'd say that points to an issue related to the kernels for primes 3 and 5, but I have not checked so far.

[mfherbst]

Looking at the code in more detail, I see that there are different code paths for Complex{Float32} and Complex{Float64} and more general floating point types, so that's that mystery solved. I also tried to track down a little further and I think the error occurs only if NontwiddleBluesteinStep / TwiddleBluesteinStep are involved.

Update: I narrowed it down to TwiddleBluesteinStep, that's why 9 is the first failing number.

[mfherbst]

Potentially related: For multidimensional FFTs there are more sizes, that give wrong results. For example an array of size (4, 32) (i.e. with pure 2^n sizes) gives the wrong results for non-Float32 and non-Float64.