From 50f81dc0bbc165d8af731986c0e5d0126452e3dd Mon Sep 17 00:00:00 2001
From: Sheehan Olver <solver@mac.com>
Date: Wed, 19 Apr 2023 21:16:51 +0100
Subject: [PATCH] Laguerre Hilbert

---
 src/stieltjes.jl     | 21 +++++++++++++++++++--
 test/test_hilbert.jl | 11 +++++++++++
 2 files changed, 30 insertions(+), 2 deletions(-)

diff --git a/src/stieltjes.jl b/src/stieltjes.jl
index e0c5b36..992e695 100644
--- a/src/stieltjes.jl
+++ b/src/stieltjes.jl
@@ -80,6 +80,24 @@ end
     log.(x .+ one(T)) .- log.(one(T) .- x)
 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}, wT::Weighted{<:Any,<:ChebyshevT})
     T = promote_type(eltype(H), eltype(wT))
     ChebyshevU{T}() * _BandedMatrix(Fill(-convert(T,π),1,∞), ℵ₀, -1, 1)
@@ -91,8 +109,7 @@ end
 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]
diff --git a/test/test_hilbert.jl b/test/test_hilbert.jl
index 4064070..4dc9a6f 100644
--- a/test/test_hilbert.jl
+++ b/test/test_hilbert.jl
@@ -273,6 +273,17 @@ end
         @test inv.(t .- x')*f ≈ -2.797995066227555
         @test log.(abs.(t .- x'))*f ≈ -5.9907385495482821485
     end
+
+    @testset "Laguerre" begin
+        x = axes(Laguerre(), 1)
+        @test (inv.(x .- x') * LaguerreWeight())[0.1] ≈ exp(-0.1)*expinti(0.1)
+        @test (inv.(x .- x') * LaguerreWeight(-0.5))[0.1] ≈ 3.3177694149902   # Mathematica
+        @test (inv.(x .- x') * LaguerreWeight(-0.5))[2.1] ≈ 1.08294830124810  # Mathematica
+        @test (inv.(x .- x') * LaguerreWeight(0.5))[0.1]  ≈ -1.44067690942561 # Mathematica
+        @test (inv.(x .- x') * LaguerreWeight(0.5))[2.1]  ≈ 0.501737581714500 # Mathematica
+        @test (inv.(x .- x') * LaguerreWeight(3.5))[0.1]  ≈ -3.46658795669242 # Mathematica
+        @test (inv.(x .- x') * LaguerreWeight(3.5))[2.1]  ≈ -5.37663478253357 # Mathematica
+    end
 end
 
 @testset "three-interval" begin