Laguerre Derivatives

src/classical/laguerre.jl

@@ -9,7 +9,8 @@ end
 
 LaguerreWeight{T}() where T = LaguerreWeight{T}(zero(T))
 LaguerreWeight() = LaguerreWeight{Float64}()
-axes(::LaguerreWeight{T}) where T = (Inclusion(ℝ),)
+# axes(::LaguerreWeight{T}) where T = (Inclusion(ℝ),)
+axes(::LaguerreWeight{T}) where T = (Inclusion(HalfLine{T}()),)
 function getindex(w::LaguerreWeight, x::Number)
     x ∈ axes(w,1) || throw(BoundsError())
     x^w.α * exp(-x)

src/stieltjes.jl

@@ -94,9 +94,25 @@ end
     ChebyshevT{T}() * _BandedMatrix(Fill(convert(T,π),1,∞), ℵ₀, 1, -1)
 end
 
+@simplify function *(H::Hilbert{<:Any,<:HalfLine,<:HalfLine}, w::LaguerreWeight)
+    T = promote_type(eltype(H), eltype(w))
+    α = convert(T, w.α)
+    z = axes(w,1)
+
+    α == 0 && return real(exp.(log.(expinti.(z) .+ zero(T)*im) .-z))
+
+    β = -1 < α < 0 ? α : α - ceil(α)
+    # S = (2* π * im) * (-z)^β*exp(-z)*gamma(-β,-z)/(gamma(-β)*(exp(im*π*β)-exp(-im*π*β)))
+    S = (2 * convert(T,π) * im) .* exp.(β*log.(-z .+ zero(T)*im) .- z + loggamma.(-β, -z .+ zero(T)*im)) ./ (gamma(-β)*(exp(im*convert(T,π)*β)-exp(-im*convert(T,π)*β)))
+    if α > 0
+        for j in 1:Int(ceil(α))
+            S = -gamma(β+j) .+ z .* S
+        end
+    end
+    real(S)
+end
 
-
-@simplify function *(H::Hilbert{<:Any,<:ChebyshevInterval,<:ChebyshevInterval}, wP::Weighted{<:Any,<:OrthogonalPolynomial})
+@simplify function *(H::Hilbert{<:Any,<:Any,<:Any}, wP::Weighted{<:Any,<:OrthogonalPolynomial})
     P = wP.P
     w = orthogonalityweight(P)
     A = recurrencecoefficients(P)[1]

test/test_stieltjes.jl

@@ -141,6 +141,15 @@ end
             @test H * ChebyshevTWeight() ≡ QuasiZeros{Float64}((x,))
             @test H * ChebyshevUWeight() == π*x
             @test (H * LegendreWeight())[0.1] ≈ log((0.1+1)/(1-0.1))
+
+            x = axes(Laguerre(), 1)
+            @test (inv.(x .- x') * LaguerreWeight())[0.1] ≈ exp(-0.1)*expinti(0.1)
+            @test isapprox((inv.(x .- x') * LaguerreWeight(-0.5))[0.1],  3.317769, rtol=1e-4) # Mathematica
+            @test isapprox((inv.(x .- x') * LaguerreWeight(-0.5))[2.1],  1.0829, rtol=1e-4) # Mathematica
+            @test isapprox((inv.(x .- x') * LaguerreWeight(0.5))[0.1],  -1.44067, rtol=1e-4) # Mathematica
+            @test isapprox((inv.(x .- x') * LaguerreWeight(0.5))[2.1],  0.5017, rtol=1e-3) # Mathematica
+            @test isapprox((inv.(x .- x') * LaguerreWeight(3.5))[0.1],  -3.46658, rtol=1e-4) # Mathematica
+            @test isapprox((inv.(x .- x') * LaguerreWeight(3.5))[2.1],  -5.37663, rtol=1e-4) # Mathematica
         end
 
         @test (Ultraspherical(1) \ (H*wT))[1:10,1:10] == diagm(1 => fill(-π,9))