export jac2jac

Author: dlfivefifty

src/FastTransforms.jl

@@ -41,7 +41,7 @@ import AbstractFFTs: Plan
 import Compat: range, transpose, adjoint, axes
 
 export cjt, icjt, jjt, plan_cjt, plan_icjt
-export leg2cheb, cheb2leg, leg2chebu, ultra2ultra, jac2jac
+export leg2cheb, cheb2leg, leg2chebu, ultra2ultra, jac2jac, plan_jac2jac
 export normleg2cheb, cheb2normleg, normleg12cheb2, cheb22normleg1
 export plan_leg2cheb, plan_cheb2leg
 export plan_normleg2cheb, plan_cheb2normleg
@@ -99,9 +99,11 @@ include("toeplitzhankel.jl")
 
 #leg2cheb(x...)=th_leg2cheb(x...)
 #cheb2leg(x...)=th_cheb2leg(x...)
-leg2chebu(x...)=th_leg2chebu(x...)
-ultra2ultra(x...)=th_ultra2ultra(x...)
-jac2jac(x...)=th_jac2jac(x...)
+leg2chebu(x...) = th_leg2chebu(x...)
+ultra2ultra(x...) = th_ultra2ultra(x...)
+jac2jac(x...) = th_jac2jac(x...)
+plan_jac2jac(x...) = th_jac2jacplan(x...)
+plan_ultra2ultra(x...) = th_ultra2ultraplan(x...)
 
 include("hierarchical.jl")
 include("SphericalHarmonics/SphericalHarmonics.jl")

src/toeplitzhankel.jl

@@ -145,15 +145,25 @@ function ultra2ultraTH(::Type{S},n,λ₁,λ₂) where S
 end
 
 function jac2jacTH(::Type{S},n,α,β,γ,δ) where S
-    @assert β == δ
-    @assert abs(α-γ) < 1
-    @assert α+β > -1
-    jk = zero(S):n-one(S)
-    DL = (2jk .+ γ .+ β .+ one(S)).*Λ.(jk,γ+β+one(S),β+one(S))
-    T = TriangularToeplitz(Λ.(jk,α-γ,one(S)),:U)
-    H = Hankel(Λ.(jk,α+β+one(S),γ+β+two(S)),Λ.(jk.+n.-one(S),α+β+one(S),γ+β+two(S)))
-    DR = Λ.(jk,β+one(S),α+β+one(S))./gamma(α-γ)
-    T,H,DL,DR
+    if β == δ
+        @assert abs(α-γ) < 1
+        @assert α+β > -1
+        jk = zero(S):n-one(S)
+        DL = (2jk .+ γ .+ β .+ one(S)).*Λ.(jk,γ+β+one(S),β+one(S))
+        T = TriangularToeplitz(Λ.(jk,α-γ,one(S)),:U)
+        H = Hankel(Λ.(jk,α+β+one(S),γ+β+two(S)),Λ.(jk.+n.-one(S),α+β+one(S),γ+β+two(S)))
+        DR = Λ.(jk,β+one(S),α+β+one(S))./gamma(α-γ)
+        T,H,DL,DR
+    elseif α == γ
+        T,H,DL,DR = jac2jacTH(S,n,β,α,δ,γ)
+        ve = T.ve
+        @inbounds for k = 2:2:length(ve)
+            ve[k] *= -1
+        end
+        TriangularToeplitz(ve, :U),H,DL,DR
+    else
+        throw(ArgumentError("Cannot create Toeplitz dot Hankel, use a sequence of plans."))
+    end
 end
 
 struct ChebyshevToLegendrePlanTH{TH}