Fix highest coefficient normalization in Chebyshev U transforms

Project.toml

@@ -1,6 +1,6 @@
 name = "FastTransforms"
 uuid = "057dd010-8810-581a-b7be-e3fc3b93f78c"
-version = "0.17.3"
+version = "0.17.4"
 
 
 [deps]

src/chebyshevtransform.jl

@@ -408,7 +408,11 @@ end
 function *(P::ChebyshevUTransformPlan{T,1,K,true,N}, x::AbstractArray{T,N}) where {T,K,N}
     length(x) ≤ 1 && return x
     _chebu1_prescale!(P.plan.region, x)
-    P.plan * x
+    x = P.plan * x
+    for d in P.plan.region
+        ldiv_dim_end!(2, d, x)
+    end
+    x
 end
 
 function mul!(y::AbstractArray{T}, P::ChebyshevUTransformPlan{T,1,K,false}, x::AbstractArray{T}) where {T,K}
@@ -418,7 +422,7 @@ function mul!(y::AbstractArray{T}, P::ChebyshevUTransformPlan{T,1,K,false}, x::A
     _plan_mul!(y, P.plan, x)
     _chebu1_postscale!(P.plan.region, x)
     for d in P.plan.region
-        size(y,d) == 1 && ldiv!(2, y) # fix doubling
+        ldiv_dim_end!(2, d, y)
     end
     y
 end
@@ -527,18 +531,24 @@ end
 
 function *(P::IChebyshevUTransformPlan{T,1,K,true}, x::AbstractArray{T}) where {T<:fftwNumber,K}
     length(x) ≤ 1 && return x
+    for d in P.plan.region
+        lmul_dim_end!(2, d, x)
+    end
     x = P.plan * x
     _ichebu1_postscale!(P.plan.region, x)
 end
 
 function mul!(y::AbstractArray{T}, P::IChebyshevUTransformPlan{T,1,K,false}, x::AbstractArray{T}) where {T<:fftwNumber,K}
     size(y) == size(x) || throw(DimensionMismatch("output must match dimension"))
     isempty(x) && return y
+    for d in P.plan.region
+        lmul_dim_end!(2, d, x)
+    end
     _plan_mul!(y, P.plan, x)
-    _ichebu1_postscale!(P.plan.region, y)
     for d in P.plan.region
-        size(y,d) == 1 && lmul!(2, y) # fix doubling
+        ldiv_dim_end!(2, d, x)
     end
+    _ichebu1_postscale!(P.plan.region, y)
     y
 end
 

test/chebyshevtests.jl

@@ -213,6 +213,16 @@ using FastTransforms, Test
             @test ichebyshevutransform!(v) === v
         end
     end
+    @testset "Chebyshev U highest coefficient" begin
+        for T in (Float32, Float64, ComplexF32, ComplexF64)
+            f = T[2,-1,2]
+            g = T[0,0,1]
+            @test chebyshevutransform(f) ≈ g
+            @test chebyshevutransform!(copy(f)) ≈ g
+            @test ichebyshevutransform(g) ≈ f
+            @test ichebyshevutransform!(copy(g)) ≈ f
+        end
+    end
     @testset "Chebyshev second kind points <-> second kind coefficients" begin
         for T in (Float32, Float64, ComplexF32, ComplexF64)
             n = 20