some rearrangements

close the longest open pull request? #89, fixing #67

also, stop the type piracy of DSP's conv for certain bitstype floating-point types

Author: MikaelSlevinsky

src/FastTransforms.jl

@@ -1,24 +1,23 @@
 module FastTransforms
 
-using FastGaussQuadrature, LinearAlgebra
-using Reexport, SpecialFunctions, ToeplitzMatrices, FillArrays, ArrayLayouts
+using ArrayLayouts, FastGaussQuadrature, FillArrays, LinearAlgebra,
+      Reexport, SpecialFunctions, ToeplitzMatrices
 
 import DSP
 
 @reexport using AbstractFFTs
 @reexport using FFTW
 
-import Base: unsafe_convert, eltype, ndims, adjoint, transpose, show, *, \,
-             inv, length, size, view, getindex, convert
+import Base: convert, unsafe_convert, eltype, ndims, adjoint, transpose, show,
+             *, \, inv, length, size, view, getindex
 
 import Base.GMP: Limb
 
 import AbstractFFTs: Plan, ScaledPlan,
-                     fft, ifft, bfft, fft!, ifft!, bfft!,
-                     plan_fft, plan_ifft, plan_bfft, plan_fft!, plan_ifft!, plan_bfft!,
-                     rfft, irfft, brfft, plan_rfft, plan_irfft, plan_brfft,
-                     fftshift, ifftshift,
-                     rfft_output_size, brfft_output_size,
+                     fft, ifft, bfft, fft!, ifft!, bfft!, rfft, irfft, brfft,
+                     plan_fft, plan_ifft, plan_bfft, plan_fft!, plan_ifft!,
+                     plan_bfft!, plan_rfft, plan_irfft, plan_brfft,
+                     fftshift, ifftshift, rfft_output_size, brfft_output_size,
                      plan_inv, normalization
 
 import DSP: conv
@@ -53,28 +52,37 @@ include("clenshaw.jl")
 
 include("libfasttransforms.jl")
 
-export plan_nufft, plan_nufft1, plan_nufft2, plan_nufft3, plan_inufft1, plan_inufft2
 export nufft, nufft1, nufft2, nufft3, inufft1, inufft2
 
+export plan_nufft, plan_nufft1, plan_nufft2, plan_nufft3,
+       plan_inufft1, plan_inufft2
+
 include("nufft.jl")
 include("inufft.jl")
 
 export paduatransform, ipaduatransform, paduatransform!, ipaduatransform!,
-       paduapoints, plan_paduatransform!, plan_ipaduatransform!
+       paduapoints
+
+export plan_paduatransform!, plan_ipaduatransform!
 
 include("PaduaTransform.jl")
 
-export plan_chebyshevtransform, plan_ichebyshevtransform, plan_chebyshevtransform!, plan_ichebyshevtransform!,
-            chebyshevtransform, ichebyshevtransform, chebyshevpoints,
-            plan_chebyshevutransform, plan_ichebyshevutransform, plan_chebyshevutransform!, plan_ichebyshevutransform!,
-            chebyshevutransform, ichebyshevutransform,
-            chebyshevtransform!, ichebyshevtransform!, chebyshevutransform!, ichebyshevutransform!
+export chebyshevtransform, ichebyshevtransform,
+       chebyshevtransform!, ichebyshevtransform!,
+       chebyshevutransform, ichebyshevutransform,
+       chebyshevutransform!, ichebyshevutransform!, chebyshevpoints
+
+export plan_chebyshevtransform, plan_ichebyshevtransform,
+       plan_chebyshevtransform!, plan_ichebyshevtransform!,
+       plan_chebyshevutransform, plan_ichebyshevutransform,
+       plan_chebyshevutransform!, plan_ichebyshevutransform!
 
 include("chebyshevtransform.jl")
 
-export plan_clenshawcurtis, clenshawcurtisnodes, clenshawcurtisweights
-export plan_fejer1, fejernodes1, fejerweights1,
-       plan_fejer2, fejernodes2, fejerweights2
+export clenshawcurtisnodes, clenshawcurtisweights, fejernodes1, fejerweights1,
+       fejernodes2, fejerweights2
+
+export plan_clenshawcurtis, plan_fejer1, plan_fejer2
 
 include("clenshawcurtis.jl")
 include("fejer.jl")
@@ -97,10 +105,6 @@ export sphones, sphzeros, sphrand, sphrandn, sphevaluate,
        tetones, tetzeros, tetrand, tetrandn,
        spinsphones, spinsphzeros, spinsphrand, spinsphrandn
 
-lgamma(x) = logabsgamma(x)[1]
-
 include("specialfunctions.jl")
 
-
-
 end # module

src/fftBigFloat.jl

@@ -51,7 +51,7 @@ function generic_fft(x::Vector{T}) where T<:AbstractFloats
     ks = range(zero(real(T)),stop=n-one(real(T)),length=n)
     Wks = exp.((-im).*convert(T,π).*ks.^2 ./ n)
     xq, wq = x.*Wks, conj([exp(-im*convert(T,π)*n);reverse(Wks);Wks[2:end]])
-    return Wks.*conv(xq,wq)[n+1:2n]
+    return Wks.*_conv!(xq,wq)[n+1:2n]
 end
 
 generic_bfft(x::StridedArray{T, N}, region) where {T <: AbstractFloats, N} = conj!(generic_fft(conj(x), region))
@@ -69,16 +69,23 @@ function generic_irfft(v::Vector{T}, n::Integer, region) where T<:ComplexFloats
 end
 generic_brfft(v::StridedArray, n::Integer, region) = generic_irfft(v, n, region)*n
 
-function conv(u::StridedVector{T}, v::StridedVector{T}) where T<:AbstractFloats
-    nu,nv = length(u),length(v)
+function _conv!(u::StridedVector{T}, v::StridedVector{T}) where T<:AbstractFloats
+    nu = length(u)
+    nv = length(v)
     n = nu + nv - 1
-    np2 = nextpow(2,n)
-    append!(u,zeros(T,np2-nu)),append!(v,zeros(T,np2-nv))
+    np2 = nextpow(2, n)
+    append!(u, zeros(T, np2-nu))
+    append!(v, zeros(T, np2-nv))
     y = generic_ifft_pow2(generic_fft_pow2(u).*generic_fft_pow2(v))
     #TODO This would not handle Dual/ComplexDual numbers correctly
     y = T<:Real ? real(y[1:n]) : y[1:n]
 end
 
+conv(u::AbstractArray{T, N}, v::AbstractArray{T, N}) where {T<:AbstractFloat, N} = _conv!(deepcopy(u), deepcopy(v))
+conv(u::AbstractArray{T, N}, v::AbstractArray{Complex{T}, N}) where {T<:AbstractFloat, N} = _conv!(complex(deepcopy(u)), deepcopy(v))
+conv(u::AbstractArray{Complex{T}, N}, v::AbstractArray{T, N}) where {T<:AbstractFloat, N} = _conv!(deepcopy(u), complex(deepcopy(v)))
+conv(u::AbstractArray{Complex{T}, N}, v::AbstractArray{Complex{T}, N}) where {T<:AbstractFloat, N} = _conv!(deepcopy(u), deepcopy(v))
+
 # This is a Cooley-Tukey FFT algorithm inspired by many widely available algorithms including:
 # c_radix2.c in the GNU Scientific Library and four1 in the Numerical Recipes in C.
 # However, the trigonometric recurrence is improved for greater efficiency.
@@ -120,16 +127,16 @@ function generic_fft_pow2!(x::Vector{T}) where T<:AbstractFloat
 end
 
 function generic_fft_pow2(x::Vector{Complex{T}}) where T<:AbstractFloat
-    y = interlace(real(x),imag(x))
+    y = interlace(real(x), imag(x))
     generic_fft_pow2!(y)
-    return complex.(y[1:2:end],y[2:2:end])
+    return complex.(y[1:2:end], y[2:2:end])
 end
-generic_fft_pow2(x::Vector{T}) where {T<:AbstractFloat} = generic_fft_pow2(complex(x))
+generic_fft_pow2(x::Vector{T}) where T<:AbstractFloat = generic_fft_pow2(complex(x))
 
 function generic_ifft_pow2(x::Vector{Complex{T}}) where T<:AbstractFloat
-    y = interlace(real(x),-imag(x))
+    y = interlace(real(x), -imag(x))
     generic_fft_pow2!(y)
-    return complex.(y[1:2:end],-y[2:2:end])/length(x)
+    return ldiv!(length(x), conj!(complex.(y[1:2:end], y[2:2:end])))
 end
 
 function generic_dct(x::StridedVector{T}, region::Integer) where T<:AbstractFloats

src/specialfunctions.jl

@@ -56,6 +56,8 @@ function pochhammer(x::Number,n::UnitRange{T}) where T<:Real
     ret
 end
 
+lgamma(x) = logabsgamma(x)[1]
+
 ogamma(x::Number) = (isinteger(x) && x<0) ? zero(float(x)) : inv(gamma(x))
 
 """

test/fftBigFloattests.jl

@@ -1,4 +1,4 @@
-using FastTransforms, FFTW, Test
+using DSP, FFTW, FastTransforms, LinearAlgebra, Test
 
 @testset "BigFloat FFT and DCT" begin
 
@@ -9,6 +9,12 @@ using FastTransforms, FFTW, Test
     c = collect(range(-big(1.0),stop=1.0,length=201))
     @test norm(ifft(fft(c))-c) < 200norm(c)eps(BigFloat)
 
+    s = big(1) ./ (1:10)
+    s64 = Float64.(s)
+    @test Float64.(conv(s, s)) ≈ conv(s64, s64)
+    @test s == big(1) ./ (1:10) #67, ensure conv doesn't overwrite input
+    @test all(s64 .=== Float64.(big(1) ./ (1:10)))
+
     p = plan_dct(c)
     @test norm(FastTransforms.generic_dct(c) - p*c) == 0
 

test/gaunttests.jl

@@ -1,4 +1,4 @@
-using FastTransforms, Test
+using FastTransforms, LinearAlgebra, Test
 
 import FastTransforms: δ
 

test/nuffttests.jl

@@ -1,4 +1,4 @@
-using FastTransforms, Test, FFTW
+using FFTW, FastTransforms, LinearAlgebra, Test
 
 FFTW.set_num_threads(ceil(Int, Sys.CPU_THREADS/2))
 
@@ -75,7 +75,7 @@ FFTW.set_num_threads(ceil(Int, Sys.CPU_THREADS/2))
         fftc = fft(c)
         if Sys.WORD_SIZE == 64
             @test_skip norm(nufft1(c, ω, ϵ) - fftc) == 0 # skip because fftw3 seems to change this
-            @test_skip norm(nufft2(c, x, ϵ) - fftc) == 0 # skip because fftw3 seems to change this
+            @test norm(nufft2(c, x, ϵ) - fftc) == 0
             @test_skip norm(nufft3(c, x, ω, ϵ) - fftc) == 0 # skip because fftw3 seems to change this
         end
         err_bnd = 500*eps(Float64)*norm(c)

test/quadraturetests.jl

@@ -1,4 +1,4 @@
-using FastTransforms, Test
+using FastTransforms, LinearAlgebra, Test
 
 import FastTransforms: chebyshevmoments1, chebyshevmoments2,
                        chebyshevjacobimoments1, chebyshevjacobimoments2,

test/specialfunctionstests.jl

@@ -1,4 +1,4 @@
-using FastTransforms, Test
+using FastTransforms, LinearAlgebra, Test
 
 import FastTransforms: pochhammer, sqrtpi, SpecialFunctions.gamma
 import FastTransforms: Cnλ, Λ, lambertw, Cnαβ, Anαβ