Fix jacobi integral where one order needs shifting (#261)

* fix jacobi integral where one order needs shifting

* version bump to v0.6.34

Author: jishnub

Project.toml

@@ -1,6 +1,6 @@
 name = "ApproxFunOrthogonalPolynomials"
 uuid = "b70543e2-c0d9-56b8-a290-0d4d6d4de211"
-version = "0.6.33"
+version = "0.6.34"
 
 [deps]
 ApproxFunBase = "fbd15aa5-315a-5a7d-a8a4-24992e37be05"

src/Spaces/Jacobi/Jacobi.jl

@@ -88,10 +88,18 @@ reverseorientation(f::Fun{<:Jacobi}) =
 # p_{n+1} = (A_n x + B_n)p_n - C_n p_{n-1}
 #####
 @inline function jacobirecA(::Type{T},α,β,k)::T where T
-    k==0&&((α+β==0)||(α+β==-1)) ? (α+β)/2+one(T) : (2k+α+β+one(T))*(2k+α+β+2one(T))/(2*(k+one(T))*(k+α+β+one(T)))
+    if k==0 && ((α+β==0)||(α+β==-1))
+        (α+β)/2+one(T)
+    else
+        (2k+α+β+one(T))*(2k+α+β+2one(T))/(2*(k+one(T))*(k+α+β+one(T)))
+    end
 end
 @inline function jacobirecB(::Type{T},α,β,k)::T where T
-    k==0&&((α+β==0)||(α+β==-1)) ? (α-β)*one(T)/2 : (α-β)*(α+β)*(2k+α+β+one(T))/(2*(k+one(T))*(k+α+β+one(T))*(2one(T)*k+α+β))
+    if k==0 && ((α+β==0)||(α+β==-1))
+        (α-β)*one(T)/2
+    else
+        (α-β)*(α+β)*(2k+α+β+one(T))/(2*(k+one(T))*(k+α+β+one(T))*(2one(T)*k+α+β))
+    end
 end
 @inline function jacobirecC(::Type{T},α,β,k)::T where T
     (one(T)*k+α)*(one(T)*k+β)*(2k+α+β+2one(T))/((k+one(T))*(k+α+β+one(T))*(2one(T)*k+α+β))

src/Spaces/Jacobi/JacobiOperators.jl

@@ -52,9 +52,11 @@ end
     elseif J.a > 0 && J.b > 0   # we have a simple definition
         ConcreteIntegral(J,1)
     else   # convert and then integrate
-        a_max = maximum(J.a:1:(1 + (J.a > 1)))
-        b_max = maximum(J.b:1:(1 + (J.b > 1)))
-        sp=Jacobi(b_max,a_max,domain(J))
+        abmin = min(J.a, J.b)
+        nsteps = length(abmin:0)
+        a = J.a + nsteps
+        b = J.b + nsteps
+        sp=Jacobi(b,a,domain(J))
         C=_conversion_shiftordersbyone(J,sp)
         Qconc=ConcreteIntegral(sp,1)
         IntegralWrapper(TimesOperator(Qconc,C),1,J,rangespace(Qconc))

test/JacobiTest.jl

@@ -743,7 +743,8 @@ include("testutils.jl")
         @inferred (() -> Integral(NormalizedLegendre()))()
         @inferred (() -> Integral(NormalizedJacobi(1,1)))()
         @inferred (() -> Integral(NormalizedJacobi(NormalizedUltraspherical(1))))()
-        for sp in (Legendre(), Jacobi(1,1), Jacobi(Ultraspherical(1)))
+        for sp in (Legendre(), Jacobi(1,1), Jacobi(Ultraspherical(1)),
+                    Jacobi(0,2), Jacobi(2,0), Jacobi(-0.5,0))
             @testset for _sp in (sp, NormalizedPolynomialSpace(sp))
                 Ij = Integral(_sp, 1)
                 @test !isdiag(Ij)