Add chebyshev u transform

Author: dlfivefifty

examples/chebyshev.jl

@@ -3,6 +3,8 @@
 # explaining precisely the normalization and points
 #############
 
+using FastTransforms
+
 # first kind points -> first kind polynomials
 n = 20
 p_1 = [sinpi((n-2k-1)/2n) for k=0:n-1]
@@ -12,6 +14,9 @@ f̌ = chebyshevtransform(f; kind=1)
 f̃ = x -> [cos(k*acos(x)) for k=0:n-1]' * f̌
 f̃(0.1) ≈ exp(0.1)
 
+# first kind polynomials -> first kind points
+ichebyshevtransform(f̌; kind=1) ≈ exp.(p_1)
+
 # second kind points -> first kind polynomials
 p_2 = [cospi(k/(n-1)) for k=0:n-1]
 f = exp.(p_2)
@@ -20,8 +25,28 @@ f̌ = chebyshevtransform(f; kind=2)
 f̃ = x -> [cos(k*acos(x)) for k=0:n-1]' * f̌
 f̃(0.1) ≈ exp(0.1)
 
-# first kind polynomials -> first kind points
-ichebyshevtransform(f̌; kind=1) ≈ exp.(p_1)
-
 # first kind polynomials -> second kind points
 ichebyshevtransform(f̌; kind=2) ≈ exp.(p_2)
+
+
+# first kind points -> second kind polynomials
+n = 20
+p_1 = [sinpi((n-2k-1)/2n) for k=0:n-1]
+f = exp.(p_1)
+f̌ = chebyshevutransform(f; kind=1)
+f̃ = x -> [sin((k+1)*acos(x))/sin(acos(x)) for k=0:n-1]' * f̌
+f̃(0.1) ≈ exp(0.1)
+
+# second kind polynomials -> first kind points
+ichebyshevutransform(f̌; kind=1) ≈ exp.(p_1)
+
+
+# second kind points -> second kind polynomials
+p_2 = [cospi(k/(n-1)) for k=1:n-2]
+f = exp.(p_2)
+f̌ = chebyshevutransform(f; kind=2)
+f̃ = x -> [sin((k+1)*acos(x))/sin(acos(x)) for k=0:n-3]' * f̌
+f̃(0.1) ≈ exp(0.1)
+
+# second kind polynomials -> second kind points
+ichebyshevutransform(f̌; kind=2) ≈ exp.(p_2)

src/FastTransforms.jl

@@ -34,7 +34,7 @@ import HierarchicalMatrices: mul!, At_mul_B!, Ac_mul_B!
 import HierarchicalMatrices: ThreadSafeVector, threadsafezeros
 import LowRankApprox: ColPerm
 import AbstractFFTs: Plan
-import Compat: range, transpose, adjoint
+import Compat: range, transpose, adjoint, axes
 
 export cjt, icjt, jjt, plan_cjt, plan_icjt
 export leg2cheb, cheb2leg, leg2chebu, ultra2ultra, jac2jac

src/chebyshevtransform.jl

@@ -187,7 +187,7 @@ ChebyshevUTransformPlan{k,inp}(plan) where {k,inp} =
 
 function plan_chebyshevutransform!(x::AbstractVector{T}; kind::Integer=1) where T<:fftwNumber
     if kind == 1
-        plan = plan_r2r!(x, REDFT10)
+        plan = plan_r2r!(x, RODFT10)
         ChebyshevUTransformPlan{1,true}(plan)
     elseif kind == 2
         if length(x) ≤ 1
@@ -205,13 +205,10 @@ end
 
 function *(P::ChebyshevUTransformPlan{T,1,true},x::AbstractVector{T}) where T
     n = length(x)
-    if n == 1
-        x
-    else
-        x = P.plan*x
-        x[1]/=2
-        lmul!(inv(convert(T,n)), x)
+    for k=1:n # sqrt(1-x_j^2) weight
+        x[k] *= sinpi(one(T)/(2n) + (k-one(T))/n)/n
     end
+    P.plan * x
 end
 
 function *(P::ChebyshevUTransformPlan{T,2,true},x::AbstractVector{T}) where T
@@ -242,28 +239,42 @@ function plan_ichebyshevutransform!(x::AbstractVector{T};kind::Integer=1) where
         if length(x) == 0
             error("Cannot create a length 0 inverse chebyshevu transform")
         end
-        plan = plan_r2r!(x, REDFT01)
+        plan = plan_r2r!(x, RODFT01)
         IChebyshevUTransformPlan{T,1,true,typeof(plan)}(plan)
     elseif kind == 2
         if length(x) ≤ 1
             error("Cannot create a length $(length(x)) inverse chebyshevu transform")
         end
-        plan = plan_chebyshevutransform!(x;kind=2)
+        plan = plan_r2r!(x, RODFT00)
         IChebyshevUTransformPlan{T,2,true,typeof(plan)}(plan)
     end
 end
 
-function plan_ichebyshevutransform(x::AbstractVector{T};kind::Integer=1) where T<:fftwNumber
-    plan = plan_ichebyshevutransform!(similar(Vector{T},axes(x));kind=kind)
+function plan_ichebyshevutransform(x::AbstractVector{T}; kind::Integer=1) where T<:fftwNumber
+    plan = plan_ichebyshevutransform!(similar(Vector{T},axes(x)); kind=kind)
     IChebyshevUTransformPlan{T,kind,false,typeof(plan)}(plan)
 end
 
-function *(P::IChebyshevUTransformPlan{T,1,true},x::AbstractVector{T}) where T<:fftwNumber
-    error("Not implemented")
+function *(P::IChebyshevUTransformPlan{T,1,true}, x::AbstractVector{T}) where T<:fftwNumber
+    n = length(x)
+    x = P.plan * x
+    for k=1:n # sqrt(1-x_j^2) weight
+        x[k] /= 2sinpi(one(T)/(2n) + (k-one(T))/n)
+    end
+    x
 end
 
-function *(P::IChebyshevUTransformPlan{T,2,true},x::AbstractVector{T}) where T<:fftwNumber
-    error("Not implemented")
+
+
+function *(P::IChebyshevUTransformPlan{T,2,true}, x::AbstractVector{T}) where T<:fftwNumber
+    n = length(x)
+    c = one(T)/ (n+1)
+    lmul!((n+1)/(2n+2*one(T)), x)
+    x = P.plan * x
+    for k=1:n # sqrt(1-x_j^2) weight
+        x[k] /= sinpi(k*c)
+    end
+    x
 end
 
 ichebyshevutransform!(x::AbstractVector{T};kind::Integer=1) where {T<:fftwNumber} =

test/chebyshevtests.jl

@@ -26,6 +26,19 @@ using FastTransforms, Compat, Compat.Test
         end
     end
 
+    @testset "Chebyshev first kind points <-> second kind coefficients" begin
+        for T in (Float32, Float64, ComplexF32, ComplexF64)
+            n = 20
+            p_1 = chebyshevpoints(T, n, kind=1)
+            f = exp.(p_1)
+            f̌ = chebyshevutransform(f; kind=1)
+
+            f̃ = x -> [sin((k+1)*acos(x))/sin(acos(x)) for k=0:n-1]' * f̌
+            @test f̃(0.1) ≈ exp(T(0.1))
+            @test ichebyshevutransform(f̌; kind=1) ≈ exp.(p_1)
+        end
+    end
+
     @testset "Chebyshev second kind points <-> second kind coefficients" begin
         for T in (Float32, Float64, ComplexF32, ComplexF64)
             n = 20
@@ -35,7 +48,7 @@ using FastTransforms, Compat, Compat.Test
 
             f̃ = x -> [sin((k+1)*acos(x))/sin(acos(x)) for k=0:n-3]' * f̌
             @test f̃(0.1) ≈ exp(T(0.1))
-            # @test ichebyshevtransform(f̌; kind=2) ≈ exp.(p_2)
+            @test ichebyshevutransform(f̌; kind=2) ≈ exp.(p_2)
         end
     end
 end