From 7565674042b8313f92b7f2a2bd3e0c142a36decb Mon Sep 17 00:00:00 2001
From: Sheehan Olver <solver@mac.com>
Date: Tue, 22 Aug 2023 16:11:28 +0100
Subject: [PATCH] LazyBandedMatrices v0.9 (#28)

* LazyBandedMatrices v0.9

* ContinuumArrays v0.15, improvements to stieltjes

* add tests
---
 Project.toml                                |  6 ++---
 src/stieltjes.jl                            | 30 ++++++++++++---------
 test/runtests.jl                            |  2 +-
 test/{test_hilbert.jl => test_stieltjes.jl} |  7 +++--
 4 files changed, 27 insertions(+), 18 deletions(-)
 rename test/{test_hilbert.jl => test_stieltjes.jl} (96%)

diff --git a/Project.toml b/Project.toml
index fcd3786..11ceafa 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "SingularIntegrals"
 uuid = "d7440221-8b5e-42fc-909c-0567823f424a"
 authors = ["Sheehan Olver <solver@mac.com>"]
-version = "0.2"
+version = "0.2.1"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -22,13 +22,13 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 ArrayLayouts = "1.0.11"
 BandedMatrices = "0.17.17"
 ClassicalOrthogonalPolynomials = "0.11"
-ContinuumArrays = "0.14"
+ContinuumArrays = "0.14, 0.15"
 FastTransforms = "0.15"
 FillArrays = "1"
 HypergeometricFunctions = "0.3.4"
 InfiniteArrays = "0.12.11, 0.13"
 LazyArrays = "1"
-LazyBandedMatrices = "0.8.5"
+LazyBandedMatrices = "0.8.5, 0.9"
 QuasiArrays = "0.11"
 SpecialFunctions = "1, 2"
 julia = "1.9"
diff --git a/src/stieltjes.jl b/src/stieltjes.jl
index ad46495..9004e9a 100644
--- a/src/stieltjes.jl
+++ b/src/stieltjes.jl
@@ -87,17 +87,20 @@ end
 computes inv.(y - x') * P understood in a principle value sense.
 """
 stieltjes(P, y...) = stieltjes_layout(MemoryLayout(P), P, y...)
-function stieltjes_layout(lay, P, y)
+stieltjes_layout(lay, P, y) = stieltjes_size(size(P), P, y)
+function stieltjes_size(sz, P, y)
     axes(P,1) == y && return stieltjes(P)
     error("Not implemented")
 end
 
+stieltjes_size(::Tuple{Any}, P, zs::AbstractVector) = [stieltjes(P, z) for z in zs]
+
 function stieltjes_layout(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiVecOrMat, y...)
     a = arguments(LAY, V)
     *(stieltjes(a[1], y...), tail(a)...)
 end
 
-stieltjes_layout(::ExpansionLayout, A, dims...) = stieltjes_layout(ApplyLayout{typeof(*)}(), A, dims...)
+stieltjes_layout(::ExpansionLayout, A, y...) = stieltjes_layout(ApplyLayout{typeof(*)}(), A, y...)
 
 
 """
@@ -179,8 +182,7 @@ end
 
 @simplify function *(S::StieltjesPoints, w::Weight)
     zs = S.args[1].args[1] # vector of points to eval at
-    x = axes(w,1)
-    [inv.(z .- x') * w for z in zs]
+    stieltjes(w, zs)
 end
 
 function stieltjes(wP::Weighted, z::Number)
@@ -196,6 +198,17 @@ function stieltjes(wP::Weighted, z::Number)
     transpose(RecurrenceArray(z, (A,B,C), [r1,r2]))
 end
 
+function stieltjes(wP::Weighted, z::AbstractVector)
+    T = promote_type(eltype(z), eltype(wP))
+    P = wP.P
+    A,B,C = recurrencecoefficients(P)
+    w = orthogonalityweight(P)
+    data = Matrix{T}(undef, 2, length(z))
+    data[1,:] .= stieltjes(w, z) .* _p0(P)
+    data[2,:] .= (A[1] .* z .+ B[1]) .* data[1,:] .- (A[1]sum(w)*_p0(P))
+    transpose(RecurrenceArray(z, (A,B,C), data))
+end
+
 sqrtx2(z::Number) = sqrt(z-1)*sqrt(z+1)
 sqrtx2(x::Real) = sign(x)*sqrt(x^2-1)
 
@@ -204,14 +217,7 @@ stieltjes(P::Legendre, z...) = stieltjes(Weighted(P), z...)
 
 @simplify function *(S::StieltjesPoints, wP::Weighted)
     z = S.args[1].args[1] # vector of points to eval at
-    T = promote_type(eltype(S), eltype(wP))
-    P = wP.P
-    A,B,C = recurrencecoefficients(P)
-    w = orthogonalityweight(P)
-    data = Matrix{T}(undef, 2, length(z))
-    data[1,:] .= (S * w) .* _p0(P)
-    data[2,:] .= (A[1] .* z .+ B[1]) .* data[1,:] .- (A[1]sum(w)*_p0(P))
-    transpose(RecurrenceArray(z, (A,B,C), data))
+    stieltjes(wP, z)
 end
 
 @simplify function *(S::StieltjesPoints, P::Legendre)
diff --git a/test/runtests.jl b/test/runtests.jl
index bf992bc..4f5163a 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -25,7 +25,7 @@ include("test_recurrence.jl")
     @test sum(w) == 1
 end
 
-include("test_hilbert.jl")
+include("test_stieltjes.jl")
 include("test_logkernel.jl")
 include("test_power.jl")
 include("test_piecewise.jl")
\ No newline at end of file
diff --git a/test/test_hilbert.jl b/test/test_stieltjes.jl
similarity index 96%
rename from test/test_hilbert.jl
rename to test/test_stieltjes.jl
index 69f9930..6117307 100644
--- a/test/test_hilbert.jl
+++ b/test/test_stieltjes.jl
@@ -141,11 +141,14 @@ end
         @test inv.(z .- x') * f ≈ [2.8826861116458593, 1.6307809018753612]
     end
 
-    @testset "StieltjesPoints Legendre" begin
+    @testset "StieltjesPoints" begin
         z = [2.,3.]
         P = Legendre()
         x = axes(P,1)
-        @test (inv.(z .- x') * P)[:,1:5] ≈ [(inv.(2 .- x')*P)[1:5]'; (inv.(3 .- x')*P)[1:5]']
+        S = inv.(z .- x')
+        @test (S * P)[:,1:5] ≈ [(inv.(2 .- x')*P)[1:5]'; (inv.(3 .- x')*P)[1:5]']
+
+        @test S * JacobiWeight(0.1,0.2) == stieltjes.(Ref(JacobiWeight(0.1,0.2)), z)
     end
 end