Support Jacobi(-1,0)\Ultraspherical(-1/2) and Jacobi(0,-1) (#245)

* Support Jacobi(-1,0)\Ultraspherical(-1/2)

* Fix slash

* Add Jacobi(0,-1)

* Clarify

Author: DanielVandH

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.15.3"
+version = "0.15.4"
 
 
 [deps]

src/classical/ultraspherical.jl

@@ -199,8 +199,19 @@ function \(A::Ultraspherical, B::Jacobi)
     Diagonal(Ã[1,:]./A[1,:]) * (Ã\B)
 end
 function \(A::Jacobi, B::Ultraspherical)
-    B̃ = Jacobi(B)
-    (A\B̃)*Diagonal(B[1,:]./B̃[1,:])
+    if B == Ultraspherical(-1/2) && (A == Jacobi(-1, 0) || A == Jacobi(0, -1))
+        # In this case, Jacobi(-1, -1) is (currently) undefined, so the conversion via B̃ = Jacobi(B) leads to NaNs 
+        # from evaluating in B̃[1, :]
+        T = promote_type(eltype(A), eltype(B))
+        n = -2one(T) ./ (2 .* (2:∞) .- one(T))
+        sgn = A == Jacobi(-1, 0) ? one(T) : -one(T)
+        dv = Vcat(one(T), -2one(T), n)
+        ev = Vcat(-sgn, sgn .* n)
+        LazyBandedMatrices.Bidiagonal(dv, ev, :U)
+    else
+        B̃ = Jacobi(B)
+        (A\B̃)*Diagonal(B[1,:]./B̃[1,:])
+    end 
 end
 
 function \(wA::Weighted{<:Any,<:Ultraspherical}, wB::Weighted{<:Any,<:Jacobi})

test/test_ultraspherical.jl

@@ -207,4 +207,15 @@ using ClassicalOrthogonalPolynomials: grammatrix
     @testset "ChebyshevInterval constructior" begin
         @test ultraspherical(2,ChebyshevInterval()) ≡ Ultraspherical(2)
     end
+end
+
+@testset "Jacobi(-1, 0) \\ Ultraspherical(-1/2) and Jacobi(0, -1) \\ Ultraspherical(-1/2)" begin 
+    for A in (Jacobi(-1,0),Jacobi(0,-1))
+        B = Ultraspherical(-1/2)
+        R = A \ B 
+        AR = A * R 
+        lhs = AR[0.3918, 1:100]
+        rhs = B[0.3918, 1:100]
+        @test lhs ≈ rhs
+    end
 end