BigFloat chebyshev transform (#118)

* BigFloat Cheb transform

* v0.10.1

* add tests, comment out broken 2nd kind transform code

Author: dlfivefifty

Project.toml

@@ -1,6 +1,6 @@
 name = "FastTransforms"
 uuid = "057dd010-8810-581a-b7be-e3fc3b93f78c"
-version = "0.10.0"
+version = "0.10.1"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"

src/chebyshevtransform.jl

@@ -75,7 +75,7 @@ end
 *(P::ChebyshevTransformPlan{T,k,false}, x::AbstractVector{T}) where {T,k} = 
     ChebyshevTransformPlan{T,k,true}(P)*copy(x)
 
-chebyshevtransform!(x::AbstractVector{T}, kind=Val(1)) where T<:fftwNumber =
+chebyshevtransform!(x::AbstractVector{T}, kind=Val(1)) where T =
     plan_chebyshevtransform!(x, kind)*x
 
 
@@ -167,7 +167,7 @@ end
 *(P::IChebyshevTransformPlan{T,k,false},x::AbstractVector{T}) where {T,k} = 
     IChebyshevTransformPlan{T,k,true}(P)*copy(x)
 
-ichebyshevtransform!(x::AbstractVector{T}, kind=Val(1)) where {T<:fftwNumber} =
+ichebyshevtransform!(x::AbstractVector{T}, kind=Val(1)) where T =
     plan_ichebyshevtransform!(x, kind)*x
 
 ichebyshevtransform(x, kind=Val(1)) = ichebyshevtransform!(copy(x), kind)
@@ -433,3 +433,85 @@ chebyshevpoints(n::Integer, kind=Val(1)) = chebyshevpoints(Float64, n, kind)
 #     x = P.plan*x
 #     rmul!(x,half(T))
 # end
+
+
+###
+# BigFloat
+# Use `Nothing` and fall back too FFT
+###
+
+plan_chebyshevtransform(x::AbstractVector{T}, ::Val{kind}) where {T,kind} =
+    ChebyshevTransformPlan{T,kind,false,Nothing}()
+plan_ichebyshevtransform(x::AbstractVector{T}, ::Val{kind}) where {T,kind} =
+    IChebyshevTransformPlan{T,kind,false,Nothing}()
+
+plan_chebyshevtransform!(x::AbstractVector{T}, ::Val{kind}) where {T,kind} =
+    ChebyshevTransformPlan{T,kind,true,Nothing}()
+plan_ichebyshevtransform!(x::AbstractVector{T}, ::Val{kind}) where {T,kind} =
+    IChebyshevTransformPlan{T,kind,true,Nothing}()
+
+#following Chebfun's @Chebtech1/vals2coeffs.m and @Chebtech2/vals2coeffs.m
+function *(P::ChebyshevTransformPlan{T,1,false,Nothing}, x::AbstractVector{T}) where T
+    n = length(x)
+    if n == 1
+        x
+    else
+        w = [2exp(im*convert(T,π)*k/2n) for k=0:n-1]
+        ret = w.*ifft([x;reverse(x)])[1:n]
+        ret = T<:Real ? real(ret) : ret
+        ret[1] /= 2
+        ret
+    end
+end
+
+
+# function *(P::ChebyshevTransformPlan{T,2,false,Nothing}, x::AbstractVector{T}) where T
+#     n = length(x)
+#     if n == 1
+#         x
+#     else
+#         ret = ifft([x;x[end:-1:2]])[1:n]
+#         ret = T<:Real ? real(ret) : ret
+#         ret[2:n-1] *= 2
+#         ret
+#     end
+# end
+
+
+*(P::ChebyshevTransformPlan{T,1,true,Nothing}, x::AbstractVector{T}) where T =
+    copyto!(x, ChebyshevTransformPlan{T,1,false,Nothing}() * x)
+# *(P::ChebyshevTransformPlan{T,2,true,Nothing}, x::AbstractVector{T}) where T =
+#     copyto!(x, ChebyshevTransformPlan{T,2,false,Nothing}() * x)
+
+
+#following Chebfun's @Chebtech1/vals2coeffs.m and @Chebtech2/vals2coeffs.m
+function *(P::IChebyshevTransformPlan{T,1,false,Nothing}, x::AbstractVector{T}) where T
+    n = length(x)
+    if n == 1
+        x
+    else
+        w = [exp(-im*convert(T,π)*k/2n)/2 for k=0:2n-1]
+        w[1] *= 2;w[n+1] *= 0;w[n+2:end] *= -1
+        ret = fft(w.*[x;one(T);x[end:-1:2]])
+        ret = T<:Real ? real(ret) : ret
+        ret[1:n]
+    end
+end
+
+# function *(P::IChebyshevTransformPlan{T,2,true,Nothing}, x::AbstractVector{T}) where T
+#     n = length(x)
+#     if n == 1
+#         x
+#     else
+#         x[1] *= 2; x[end] *= 2
+#         chebyshevtransform!(x, Val(2))
+#         x[1] *= 2; x[end] *= 2
+#         lmul!(convert(T,n-1)/2, x)
+#         x
+#     end
+# end
+
+*(P::IChebyshevTransformPlan{T,1,true,Nothing}, x::AbstractVector{T}) where T =
+    copyto!(x, IChebyshevTransformPlan{T,1,false,Nothing}() * x)
+# *(P::IChebyshevTransformPlan{T,2,false,Nothing}, x::AbstractVector{T}) where T =
+#     IChebyshevTransformPlan{T,2,true,Nothing}() * copy(x)

src/libfasttransforms.jl

@@ -595,12 +595,12 @@ function lmul!(p::FTPlan{Complex{Float64}, 2, SPINSPHEREANALYSIS}, x::Matrix{Com
     return x
 end
 
-*(p::FTPlan{T}, x::Array{T}) where T = lmul!(p, deepcopy(x))
-*(p::AdjointFTPlan{T}, x::Array{T}) where T = lmul!(p, deepcopy(x))
-*(p::TransposeFTPlan{T}, x::Array{T}) where T = lmul!(p, deepcopy(x))
-\(p::FTPlan{T}, x::Array{T}) where T = ldiv!(p, deepcopy(x))
-\(p::AdjointFTPlan{T}, x::Array{T}) where T = ldiv!(p, deepcopy(x))
-\(p::TransposeFTPlan{T}, x::Array{T}) where T = ldiv!(p, deepcopy(x))
+*(p::FTPlan{T}, x::AbstractArray{T}) where T = lmul!(p, Array(x))
+*(p::AdjointFTPlan{T}, x::AbstractArray{T}) where T = lmul!(p, Array(x))
+*(p::TransposeFTPlan{T}, x::AbstractArray{T}) where T = lmul!(p, Array(x))
+\(p::FTPlan{T}, x::AbstractArray{T}) where T = ldiv!(p, Array(x))
+\(p::AdjointFTPlan{T}, x::AbstractArray{T}) where T = ldiv!(p, Array(x))
+\(p::TransposeFTPlan{T}, x::AbstractArray{T}) where T = ldiv!(p, Array(x))
 
 *(p::FTPlan{T, 1}, x::UniformScaling{S}) where {T, S} = lmul!(p, Matrix{promote_type(T, S)}(x, p.n, p.n))
 *(p::AdjointFTPlan{T, FTPlan{T, 1, K}}, x::UniformScaling{S}) where {T, S, K} = lmul!(p, Matrix{promote_type(T, S)}(x, p.parent.n, p.parent.n))

test/chebyshevtests.jl

@@ -215,4 +215,13 @@ using FastTransforms, Test
         @test ichebyshevutransform([1,2,3]) == ichebyshevutransform([1.,2,3])
         @test ichebyshevutransform([1,2,3], Val(2)) == ichebyshevutransform([1.,2,3], Val(2))
     end
+
+    @testset "BigFloat" begin
+        x = BigFloat[1,2,3]
+        @test ichebyshevtransform(chebyshevtransform(x)) ≈ x
+        @test plan_chebyshevtransform(x)x ≈ chebyshevtransform(x)
+        @test plan_ichebyshevtransform(x)x ≈ ichebyshevtransform(x)
+        @test plan_chebyshevtransform!(x)copy(x) ≈ chebyshevtransform(x)
+        @test plan_ichebyshevtransform!(x)copy(x) ≈ ichebyshevtransform(x)
+    end
 end