Fix the size one case and provide informative error message in the

size zero case.

Currently, the size one case errors and the size zero case hangs.
I've also added more tests. They revealed a bug in the
pow4 fft which I have also fixed.

Author: andreasnoack

src/algos.jl

@@ -73,7 +73,7 @@ function fft_dft!(out::AbstractVector{T}, in::AbstractVector{T}, N::Int, start_o
         tmp += in[start_in + j*stride_in]
     end
     out[start_out] = tmp
-    
+
     wk = wkn = w = convert(T, cispi(direction_sign(d)*2/N))
     @inbounds for d in 1:N-1
         tmp = in[start_in]
@@ -98,7 +98,7 @@ function fft_dft!(out::AbstractVector{Complex{T}}, in::AbstractVector{T}, N::Int
     end
     out[start_out] = convert(Complex{T}, tmpBegin)
     iseven(N) && (out[start_out + stride_out*halfN] = convert(Complex{T}, tmpHalf))
-    
+
     @inbounds for d in 1:halfN
         tmp = in[start_in]
         @inbounds for k in 1:N-1
@@ -200,7 +200,7 @@ function fft_pow4!(out::AbstractVector{T}, in::AbstractVector{U}, N::Int, start_
 
     w1 = convert(T, cispi(direction_sign(d)*2/N))
     wj = one(T)
-    
+
     w1 = convert(T, cispi(ds*2/N))
     w2 = convert(T, cispi(ds*4/N))
     w3 = convert(T, cispi(ds*6/N))
@@ -212,7 +212,7 @@ function fft_pow4!(out::AbstractVector{T}, in::AbstractVector{U}, N::Int, start_
         @muladd k2 = start_out + (k+2*m)*stride_out
         @muladd k3 = start_out + (k+3*m)*stride_out
         y_k0, y_k1, y_k2, y_k3 = out[k0], out[k1], out[k2], out[k3]
-        @muladd out[k0] = (y_k0 + y_k2*wk2) + (y_k1*wk1 + y_k3*wk2)
+        @muladd out[k0] = (y_k0 + y_k2*wk2) + (y_k1*wk1 + y_k3*wk3)
         @muladd out[k1] = (y_k0 - y_k2*wk2) + (y_k1*wk1 - y_k3*wk3) * plusi
         @muladd out[k2] = (y_k0 + y_k2*wk2) - (y_k1*wk1 + y_k3*wk3)
         @muladd out[k3] = (y_k0 - y_k2*wk2) + (y_k1*wk1 - y_k3*wk3) * minusi

src/callgraph.jl

@@ -55,18 +55,6 @@ end
 # Get the node in the graph at index i
 Base.getindex(g::CallGraph{T}, i::Int) where {T} = g.nodes[i]
 
-"""
-$(TYPEDSIGNATURES)
-Check if `N` is a power of `base`
-
-"""
-function _ispow(N, base)
-    while N % base == 0
-        N = N/base
-    end
-    return N == 1
-end
-
 """
 $(TYPEDSIGNATURES)
 Check if `N` is a power of 2 or 4
@@ -93,6 +81,9 @@ Recursively instantiate a set of `CallGraphNode`s
 
 """
 function CallGraphNode!(nodes::Vector{CallGraphNode}, N::Int, workspace::Vector{Vector{T}}, s_in::Int, s_out::Int)::Int where {T}
+    if N == 0
+        throw(DimensionMismatch("array has to be non-empty"))
+    end
     if iseven(N)
         pow = _ispow24(N)
         if !isnothing(pow)
@@ -102,15 +93,15 @@ function CallGraphNode!(nodes::Vector{CallGraphNode}, N::Int, workspace::Vector{
         end
     end
     if N % 3 == 0
-        if _ispow(N, 3)
+        if nextpow(3, N) == N
             push!(workspace, T[])
             push!(nodes, CallGraphNode(0, 0, Pow3FFT(), N, s_in, s_out))
             return 1
         end
     end
-    if isprime(N)
+    if N == 1 || isprime(N)
         push!(workspace, T[])
-        push!(nodes, CallGraphNode(0,0, DFT(),N, s_in, s_out))
+        push!(nodes, CallGraphNode(0, 0, DFT(), N, s_in, s_out))
         return 1
     end
     Ns = [first(x) for x in collect(factor(N)) for _ in 1:last(x)]

test/onedim/complex_backward.jl

@@ -8,4 +8,13 @@ test_nums = [8, 11, 15, 16, 27, 100]
         y_ref[1] = N
         @test y ≈ y_ref atol=1e-12
     end
-end
+end
+
+@testset verbose = true "against naive implementation. Size: $n" for n in 1:64
+    x = complex.(randn(n), randn(n))
+    @test naive_1d_fourier_transform(x, FFTA.FFT_BACKWARD) ≈ bfft(x)
+end
+
+@testset "error messages" begin
+    @test_throws DimensionMismatch bfft(zeros(0))
+end

test/onedim/complex_forward.jl

@@ -8,4 +8,13 @@ test_nums = [8, 11, 15, 16, 27, 100]
         y_ref[1] = N
         @test y ≈ y_ref atol=1e-12
     end
-end
+end
+
+@testset verbose = true "against naive implementation. Size: $n" for n in 1:64
+    x = complex.(randn(n), randn(n))
+    @test naive_1d_fourier_transform(x, FFTA.FFT_FORWARD) ≈ fft(x)
+end
+
+@testset "error messages" begin
+    @test_throws DimensionMismatch fft(zeros(0))
+end

test/onedim/real_backward.jl

@@ -11,4 +11,20 @@ test_nums = [8, 11, 15, 16, 27, 100]
         end
         @test y_ref ≈ y atol=1e-12
     end
-end
+end
+
+@testset verbose = true "against naive implementation. Size: $n" for n in 1:64
+    x = complex.(randn(n ÷ 2 + 1), randn(n ÷ 2 + 1))
+    x[begin] = real(x[begin])
+    if iseven(n)
+        x[end] = real(x[end])
+        xe = [x; conj.(reverse(x[begin + 1:end - 1]))]
+    else
+        xe = [x; conj.(reverse(x[begin + 1:end]))]
+    end
+    @test naive_1d_fourier_transform(xe, FFTA.FFT_BACKWARD) ≈ brfft(x, n)
+end
+
+@testset "error messages" begin
+    @test_throws DimensionMismatch brfft(zeros(ComplexF64, 0), 0)
+end

test/onedim/real_forward.jl

@@ -1,11 +1,20 @@
 using FFTA, Test
 test_nums = [8, 11, 15, 16, 27, 100]
-@testset verbose = true " forward" begin 
+@testset verbose = true " forward" begin
     for N in test_nums
         x = ones(Float64, N)
         y = rfft(x)
         y_ref = 0*y
         y_ref[1] = N
         @test y ≈ y_ref atol=1e-12
     end
+end
+
+@testset verbose = true "against naive implementation. Size: $n" for n in 1:64
+    x = randn(n)
+    @test naive_1d_fourier_transform(x, FFTA.FFT_FORWARD)[1:(n ÷ 2 + 1)] ≈ rfft(x)
+end
+
+@testset "error messages" begin
+    @test_throws DimensionMismatch rfft(zeros(0))
 end

test/runtests.jl

@@ -1,4 +1,4 @@
-using Test, Random
+using Test, Random, FFTA
 
 function padnum(m,x)
     digs = floor(Int, log10(m))
@@ -8,6 +8,42 @@ function padnum(m,x)
     String(v)
 end
 
+function naive_1d_fourier_transform(x::Vector, d::FFTA.Direction)
+    n = length(x)
+    y = zeros(Complex{Float64}, n)
+
+    for u in 0:(n - 1)
+        s = 0.0 + 0.0im
+        for v in 0:(n - 1)
+            a = FFTA.direction_sign(d) * 2π * u * v / n
+            s += x[v + 1] * exp(im * a)
+        end
+        y[u + 1] = s
+    end
+
+    return y
+end
+
+function naive_2d_fourier_transform(X::Matrix, d::FFTA.Direction)
+    rows, cols = size(X)
+    Y = zeros(Complex{Float64}, rows, cols)
+
+    for u in 0:(rows - 1)
+        for v in 0:(cols - 1)
+            s = 0.0 + 0.0im
+            for x in 0:(rows - 1)
+                for y in 0:(cols - 1)
+                    a = FFTA.direction_sign(d) * 2π * (u * x / rows + v * y / cols)
+                    s += X[x + 1, y + 1] * exp(im * a)
+                end
+            end
+            Y[u + 1, v + 1] = s
+        end
+    end
+
+    return Y
+end
+
 Random.seed!(1)
 @testset verbose = true "FFTA" begin
     @testset verbose = true "1D" begin

test/twodim/complex_backward.jl

@@ -8,4 +8,16 @@ test_nums = [8, 11, 15, 16, 27, 100]
         y_ref[1] = length(x)
         @test y ≈ y_ref
     end
-end
+end
+
+@testset verbose = true "against naive implementation" for n in 1:64
+    @testset "size: ($m, $n)" for m in n:(n + 1)
+        X = complex.(randn(m, n), randn(m, n))
+        Y = similar(X)
+        @test naive_2d_fourier_transform(X, FFTA.FFT_BACKWARD) ≈ bfft(X)
+    end
+end
+
+@testset "error messages" begin
+    @test_throws DimensionMismatch bfft(zeros(ComplexF64, 0, 0))
+end

test/twodim/complex_forward.jl

@@ -8,4 +8,16 @@ test_nums = [8, 11, 15, 16, 27, 100]
         y_ref[1] = length(x)
         @test y ≈ y_ref
     end
-end
+end
+
+@testset verbose = true "against naive implementation" for n in 1:64
+    @testset "size: ($m, $n)" for m in n:(n + 1)
+        X = complex.(randn(m, n), randn(m, n))
+        Y = similar(X)
+        @test naive_2d_fourier_transform(X, FFTA.FFT_FORWARD) ≈ fft(X)
+    end
+end
+
+@testset "error messages" begin
+    @test_throws DimensionMismatch fft(zeros(0, 0))
+end