Implement weighted conversion between b=-1 polynomials (#118)

* Implement weighted conversion for b=-1

* Fix integer case, add more tests

Author: DanielVandH

src/SemiclassicalOrthogonalPolynomials.jl

@@ -461,10 +461,26 @@ function \(w_A::WeightedSemiclassicalJacobi{T}, w_B::WeightedSemiclassicalJacobi
     Δb = B.b-A.b
     Δc = B.c-A.c
 
-    if (wA.a == A.a) && (wA.b == A.b) && (wA.c == A.c) && (wB.a == B.a) && (wB.b == B.b) && (wB.c == B.c) && isinteger(A.a) && isinteger(A.b) && isinteger(A.c) && isinteger(B.a) && isinteger(B.b) && isinteger(B.c)
-        # k = (A \ SemiclassicalJacobiWeight(A.t,Δa,Δb,Δc))[1]
-        k = sumquotient(SemiclassicalJacobiWeight(B.t,B.a,B.b,B.c),SemiclassicalJacobiWeight(A.t,A.a,A.b,A.c))
-        return (ApplyArray(*,Diagonal(Fill(k,∞)),(B \ A)))'
+    # k = (A \ SemiclassicalJacobiWeight(A.t,Δa,Δb,Δc))[1]
+    if isone(-wA.b) && isone(-wB.b)
+        @assert A.a + 1 == B.a && A.c + 1 == B.c
+        Q = SemiclassicalJacobi(B.t, B.a, one(B.b), B.c, B)
+        P = SemiclassicalJacobi(A.t, A.a, one(A.b), A.c, A)
+        wP = Weighted(P) 
+        wQ = Weighted(Q)
+        R22 = wP \ wQ
+        r11 = A.t - 1 
+        qw0 = SemiclassicalJacobiWeight(Q.t, Q.a, zero(Q.b), Q.c)
+        pw0 = SemiclassicalJacobiWeight(P.t, P.a, zero(P.b), P.c)
+        r21 = wP[:, 1:2] \ (qw0 .- r11 .* pw0)
+        d0 = Vcat(r11, R22[band(0)])
+        d1 = Vcat(r21[begin], R22[band(-1)])
+        d2 = Vcat(r21[begin+1], R22[band(-2)])
+        data = Hcat(d0, d1, d2)
+        return _BandedMatrix(data', 1:∞, 2, 0)
+    elseif (wA.a == A.a) && (wA.b == A.b) && (wA.c == A.c) && (wB.a == B.a) && (wB.b == B.b) && (wB.c == B.c) && isinteger(A.a) && isinteger(A.b) && isinteger(A.c) && isinteger(B.a) && isinteger(B.b) && isinteger(B.c)
+            k = sumquotient(SemiclassicalJacobiWeight(B.t,B.a,B.b,B.c),SemiclassicalJacobiWeight(A.t,A.a,A.b,A.c))
+            return (ApplyArray(*,Diagonal(Fill(k,∞)),(B \ A)))'
     elseif wA.a == wB.a && wA.b == wB.b && wA.c == wB.c # fallback to Christoffel–Darboux
         A \ B
     elseif wA.a+1 == wB.a && wA.b == wB.b && wA.c == wB.c

test/test_neg1b.jl

@@ -234,3 +234,18 @@ end
     @test_throws ArgumentError expand(HalfWeighted{:b}(P), g)
 end
 
+@testset "Weighted conversion between b=-1" begin
+    for (t, a, b, c) in ((2.0, 1 / 2, -1.0, 1 / 2), (2.5, 3 / 2, -1.0, 1 / 2), (2.5, 1.0, -1.0, 2.0))
+        Q = SemiclassicalJacobi(t, a, b, c)
+        P = SemiclassicalJacobi(t, a - 1, b, c - 1)
+        L = Weighted(P) \ Weighted(Q)
+        wP = SemiclassicalJacobiWeight(t, a - 1, b, c - 1)
+        wQ = SemiclassicalJacobiWeight(t, a, b, c)
+        lhs = wQ .* Q
+        rhs = wP .* (P * L)
+        x = LinRange(eps(), 1 - eps(), 250)
+        lhs_vals = lhs[x, 1:20]
+        rhs_vals = rhs[x, 1:20]
+        @test lhs_vals ≈ rhs_vals
+    end
+end