Use Jacobi transforms from FastTransforms

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.11.9"
+version = "0.11.10"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -33,7 +33,7 @@ ContinuumArrays = "0.16"
 DomainSets = "0.6"
 FFTW = "1.1"
 FastGaussQuadrature = "0.5, 1"
-FastTransforms = "0.15.2"
+FastTransforms = "0.15.10"
 FillArrays = "1"
 HypergeometricFunctions = "0.3.4"
 InfiniteArrays = "0.12.11, 0.13"

src/classical/jacobi.jl

@@ -69,6 +69,16 @@ singularitiesbroadcast(::typeof(*), V::Union{NoSingularities,SubQuasiArray}...)
 
 abstract type AbstractJacobi{T} <: OrthogonalPolynomial{T} end
 
+struct JacobiTransformPlan{T, CHEB2JAC, DCT} <: Plan{T}
+    cheb2jac::CHEB2JAC
+    chebtransform::DCT
+end
+
+JacobiTransformPlan(c2l, ct) = JacobiTransformPlan{promote_type(eltype(c2l),eltype(ct)),typeof(c2l),typeof(ct)}(c2l, ct)
+
+*(P::JacobiTransformPlan, x::AbstractArray) = P.cheb2jac*(P.chebtransform*x)
+
+
 include("legendre.jl")
 
 singularitiesbroadcast(::typeof(*), ::LegendreWeight, b::AbstractJacobiWeight) = b
@@ -94,6 +104,13 @@ Jacobi(P::Legendre{T}) where T = Jacobi(zero(T), zero(T))
 
 basis_singularities(w::JacobiWeight) = Weighted(Jacobi(w.a, w.b))
 
+function plan_grid_transform(P::Jacobi{T}, szs::NTuple{N,Int}, dims=1:N) where {T,N}
+    arr = Array{T}(undef, szs...)
+    x = grid(P, size(arr,1))
+    x, JacobiTransformPlan(FastTransforms.plan_th_cheb2jac!(arr, P.a, P.b, dims), plan_chebyshevtransform(arr, dims))
+end
+
+
 """
      jacobip(n, a, b, z)
 
@@ -179,7 +196,7 @@ ldiv(P::Jacobi{V}, f::Inclusion{T}) where {T,V} = _op_ldiv(P, f)
 ldiv(P::Jacobi{V}, f::AbstractQuasiFill{T,1}) where {T,V} = _op_ldiv(P, f)
 function transform_ldiv(P::Jacobi{V}, f::AbstractQuasiArray) where V
     T = ChebyshevT{V}()
-    pad(cheb2jac(paddeddata(T \ f), P.a, P.b), axes(P,2), tail(axes(f))...)
+    pad(FastTransforms.th_cheb2jac(paddeddata(T \ f), P.a, P.b, 1), axes(P,2), tail(axes(f))...)
 end
 
 

src/classical/legendre.jl

@@ -98,19 +98,10 @@ function transform_ldiv(::Legendre{V}, f::Union{AbstractQuasiVector,AbstractQuas
     pad(th_cheb2leg(paddeddata(dat)), axes(dat)...)
 end
 
-struct LegendreTransformPlan{T, CHEB2LEG, DCT} <: Plan{T}
-    cheb2leg::CHEB2LEG
-    chebtransform::DCT
-end
-
-LegendreTransformPlan(c2l, ct) = LegendreTransformPlan{promote_type(eltype(c2l),eltype(ct)),typeof(c2l),typeof(ct)}(c2l, ct)
-
-*(P::LegendreTransformPlan, x::AbstractArray) = P.cheb2leg*(P.chebtransform*x)
-
 function plan_grid_transform(P::Legendre{T}, szs::NTuple{N,Int}, dims=1:N) where {T,N}
     arr = Array{T}(undef, szs...)
     x = grid(P, size(arr,1))
-    x, LegendreTransformPlan(FastTransforms.plan_th_cheb2leg!(arr, dims), plan_chebyshevtransform(arr, dims))
+    x, JacobiTransformPlan(FastTransforms.plan_th_cheb2leg!(arr, dims), plan_chebyshevtransform(arr, dims))
 end
 
 

test/test_jacobi.jl

@@ -197,6 +197,7 @@ import ClassicalOrthogonalPolynomials: recurrencecoefficients, basis, MulQuasiMa
         P = Jacobi(1/2,0.)
         x = axes(P,1)
         @test (P * (P \ exp.(x)))[0.1] ≈ exp(0.1)
+        @test P[:,1:20] \ exp.(x) ≈ (P \ exp.(x))[1:20]
 
         @test P[0.1,:]' * (P \ [exp.(x) cos.(x)]) ≈ [exp(0.1) cos(0.1)]
 
@@ -219,8 +220,8 @@ import ClassicalOrthogonalPolynomials: recurrencecoefficients, basis, MulQuasiMa
             x = axes(P,1)
             u = P * (P \ exp.(x))
             @test u[0.1] ≈ exp(0.1)
-            # U = P * (P \ [exp.(x) cos.(x)]) # not good at multiscale
-            @test_broken U[0.1,:] ≈ [exp(0.1),cos(0.1)]
+            U = P * (P \ [exp.(x) cos.(x)])
+            @test U[0.1,:] ≈ [exp(0.1),cos(0.1)]
         end
 
         @testset "special cases" begin