simplify jacobi conversion (#189)

* simplify jacobi conversion

* fix Ultraspherical conversion

Author: jishnub

src/Spaces/Jacobi/JacobiOperators.jl

@@ -149,16 +149,18 @@ function Conversion(L::Jacobi,M::Jacobi)
         elseif (isapprox(M.b,L.b+1) && isapprox(M.a,L.a)) ||
             (isapprox(M.b,L.b) && isapprox(M.a,L.a+1))
             ConcreteConversion(L,M)
-        elseif M.b > L.b+1
-            ConversionWrapper(
-                TimesOperator(
-                    ConcreteConversion(Jacobi(M.b-1,M.a,dm),M),
-                    Conversion(L,Jacobi(M.b-static(1),M.a,dm))))
-        else  #if M.a >= L.a+1
-            ConversionWrapper(
-                TimesOperator(
-                    ConcreteConversion(Jacobi(M.b,M.a-1,dm),M),
-                    Conversion(L,Jacobi(M.b,M.a-static(1),dm))))
+        elseif M.b >= L.b && M.a >= L.a
+            # We split this into steps where a and b are changed by 1:
+            # Define the intermediate space J = Jacobi(M.b, L.a, dm)
+            # Conversion(L, M) == Conversion(J, M) * Conversion(L, J)
+            # Conversion(L, J) = Conversion(Jacobi(L.b, L.a, dm), Jacobi(M.b, L.a, dm))
+            # Conversion(J, M) = Conversion(Jacobi(M.b, L.a, dm), Jacobi(M.b, M.a, dm))
+            CLJ = [ConcreteConversion(Jacobi(b-1,L.a,dm), Jacobi(b, L.a, dm)) for b in M.b:-1:L.b+1]
+            CJM = [ConcreteConversion(Jacobi(M.b,a-1,dm), Jacobi(M.b, a, dm)) for a in M.a:-1:L.a+1]
+            C = [CJM; CLJ]
+            ConversionWrapper(TimesOperator(C))
+        else
+            error("Implement for $L → $M")
         end
     elseif L.a ≈ L.b ≈ 0 && M.a ≈ M.b ≈ 0.5
         Conversion(L,Ultraspherical(L),Ultraspherical(M),M)

src/Spaces/Ultraspherical/UltrasphericalOperators.jl

@@ -156,13 +156,22 @@ function Conversion(A::Ultraspherical,B::Ultraspherical)
     a=order(A); b=order(B)
     if b==a
         ConversionWrapper(Operator(I,A))
-    elseif -1 ≤ b-a ≤ 1 && (a,b) ≠ (2,1)
-        ConcreteConversion(A,B)
-    elseif b-a > 1
-        d=domain(A)
-        US=Ultraspherical(b-static(1),d)
-        v = [ConcreteConversion(Ultraspherical(i-1,d), Ultraspherical(i,d)) for i in b:-1:a+1]
-        ConversionWrapper(TimesOperator(v))
+    elseif isapproxinteger(b-a) || isapproxinteger_addhalf(b-a)
+        if -1 ≤ b-a ≤ 1 && (a,b) ≠ (2,1)
+            ConcreteConversion(A,B)
+        elseif b-a > 1
+            d=domain(A)
+            US=Ultraspherical(b-static(1),d)
+            r = b:-1:a+1
+            v = [ConcreteConversion(Ultraspherical(i-1,d), Ultraspherical(i,d)) for i in r]
+            if !(last(r) ≈ a+1)
+                vlast = ConcreteConversion(A, Ultraspherical(last(r)-1, d))
+                v = [v; vlast]
+            end
+            ConversionWrapper(TimesOperator(v))
+        else
+            throw(ArgumentError("Cannot convert from $A to $B"))
+        end
     else
         throw(ArgumentError("Cannot convert from $A to $B"))
     end

test/UltrasphericalTest.jl

@@ -69,10 +69,12 @@ using Static
         @test g ≈ Fun(x->x^2, Ultraspherical(0.5))
 
         f = Fun(x->x^2, Ultraspherical(0.5)) # Legendre
-        CLU = Conversion(Ultraspherical(0.5), Ultraspherical(2.5))
-        @test !isdiag(CLU)
-        g = CLU * f
-        @test g ≈ Fun(x->x^2, Ultraspherical(2.5))
+        for n in (2.5, 3)
+            CLU = Conversion(Ultraspherical(0.5), Ultraspherical(n))
+            @test !isdiag(CLU)
+            g = CLU * f
+            @test g ≈ Fun(x->x^2, Ultraspherical(n))
+        end
     end
 
     @testset "Normalized space" begin