Concrete derivative in Jacobi space (#334)

Author: jishnub

src/Spaces/Jacobi/JacobiOperators.jl

@@ -1,31 +1,21 @@
 ## Derivative
 
-# specialize Derivative so that this is type-inferred even without constant propagation
-Derivative(J::MaybeNormalized{<:Jacobi}) = ConcreteDerivative(J,1)
-@inline function _Derivative(J::Jacobi, k::Number)
+function Derivative(J::MaybeNormalized{<:Jacobi}, k::Number)
     assert_integer(k)
-    if k==1
-        return ConcreteDerivative(J,1)
-    else
-        d = domain(J)
-        v = [ConcreteDerivative(Jacobi(J.b+i-1, J.a+i-1, d)) for i in k:-1:1]
-        DerivativeWrapper(TimesOperator(v), k, J)
-    end
+    return ConcreteDerivative(J,k)
 end
-@static if VERSION >= v"1.8"
-    Base.@constprop :aggressive Derivative(J::Jacobi, k::Number) =
-        _Derivative(J, k)
-else
-    Derivative(J::Jacobi, k::Number) = _Derivative(J, k)
-end
-
 
 rangespace(D::ConcreteDerivative{<:MaybeNormalized{<:Jacobi}}) = Jacobi(D.space.b+D.order,D.space.a+D.order,domain(D))
 bandwidths(D::ConcreteDerivative{<:MaybeNormalized{<:Jacobi}}) = -D.order,D.order
 isdiag(D::ConcreteDerivative{<:MaybeNormalized{<:Jacobi}}) = false
 
-getindex(T::ConcreteDerivative{<:Jacobi}, k::Integer, j::Integer) =
-    j==k+1 ? eltype(T)((k+1+T.space.a+T.space.b)/complexlength(domain(T))) : zero(eltype(T))
+function getindex(D::ConcreteDerivative{<:Jacobi}, k::Integer, j::Integer)
+    m = D.order
+    sp = domainspace(D)
+    a, b = sp.a, sp.b
+    x = j + a + b
+    j == k+m ? eltype(D)(pochhammer(x, m)/complexlength(domain(D))^m) : zero(eltype(D))
+end
 
 
 # Evaluation

test/ChebyshevTest.jl

@@ -302,16 +302,10 @@ include("testutils.jl")
             s1 = NormalizedChebyshev(-1..1)
             s2 = NormalizedChebyshev()
             @test s1 == s2
-            D1 = if VERSION >= v"1.8"
-                @inferred Derivative(s1)
-            else
-                Derivative(s1)
-            end
-            D2 = if VERSION >= v"1.8"
-                @inferred Derivative(s2)
-            else
-                Derivative(s2)
-            end
+
+            D1 = @inferred Derivative(s1)
+            D2 = Derivative(s2)
+
             @test D1 isa ApproxFunBase.ConcreteDerivative
             @test D2 isa ApproxFunBase.ConcreteDerivative
             @test !isdiag(D1)
@@ -323,13 +317,28 @@ include("testutils.jl")
             f2 = Fun(f, s2)
             @test f1 == f2
             @test D1 * f1 == D2 * f2
+            @test D1 * f1 ≈ Fun(x -> 6x + 5, Chebyshev())
 
             if VERSION < v"1.10-"
                 ElT = @inferred (D1 -> eltype(D1 + D1))(D1)
                 @test ElT == eltype(D1)
             else
                 @test_broken @inferred((D1 -> eltype(D1 + D1))(D1)) == eltype(D1)
             end
+
+            ST = Chebyshev
+            for d in ((), (0..1,))
+                S1 = ST(d...)
+                for S in (S1, NormalizedPolynomialSpace(S1))
+                    @test Derivative(S) == Derivative(S,1)
+                    @test Derivative(S)^2 == Derivative(S,2)
+                    f = Fun(x->x^3, S)
+                    @test Derivative(S) * f ≈ Fun(x->3x^2, S)
+                    @test Derivative(S,2) * f ≈ Fun(x->6x, S)
+                    @test Derivative(S,3) * f ≈ Fun(x->6, S)
+                    @test Derivative(S,4) * f ≈ zeros(S)
+                end
+            end
         end
 
         @testset "space promotion" begin

test/JacobiTest.jl

@@ -792,15 +792,15 @@ include("testutils.jl")
             @test Derivative() * g == Derivative() * Fun(space(g), collect(coefficients(g)))
             @static if isdefined(FillArrays, :OneElement)
                 if coefficients(g) isa OneElement
-                    @test_broken coefficients(Derivative() * g) isa OneElement
+                    @test coefficients(Derivative() * g) isa OneElement
                 end
             end
 
             g = NormalizedJacobi(1,2.5)(3)
             @test Derivative() * g == Derivative() * Fun(space(g), collect(coefficients(g)))
             @static if isdefined(FillArrays, :OneElement)
                 if coefficients(g) isa OneElement
-                    @test_broken coefficients(Derivative() * g) isa OneElement
+                    @test coefficients(Derivative() * g) isa OneElement
                 end
             end
         end

test/MiscAFBTest.jl

@@ -253,21 +253,6 @@ Base.:(==)(a::UniqueInterval, b::UniqueInterval) = (@assert a.parentinterval ==
 
         @testset "Derivative" begin
             @test Derivative() == Derivative()
-            for d in Any[(), (0..1,)]
-                for ST in Any[Chebyshev, Legendre,
-                        (x...) -> Jacobi(2,2,x...), (x...) -> Jacobi(1.5,2.5,x...)]
-                    S1 = ST(d...)
-                    @testset for S in [S1, NormalizedPolynomialSpace(S1)]
-                        @test Derivative(S) == Derivative(S,1)
-                        @test Derivative(S)^2 == Derivative(S,2)
-                        f = Fun(x->x^3, S)
-                        @test Derivative(S) * f ≈ Fun(x->3x^2, S)
-                        @test Derivative(S,2) * f ≈ Fun(x->6x, S)
-                        @test Derivative(S,3) * f ≈ Fun(x->6, S)
-                        @test Derivative(S,4) * f ≈ zeros(S)
-                    end
-                end
-            end
             @test Derivative(Chebyshev()) != Derivative(Chebyshev(), 2)
             @test Derivative(Chebyshev()) != Derivative(Legendre())
         end