Fix #38 (#40)

* logkernel defaults to complexlogkernel

* Update Project.toml

* better broadcasting with real

* Weighted(ChebyshevU) now uses RecurrenceArray

* Tests pass

* Weighted Chebyshev T

* split out recurrence from IC code

* Fix recurrence forward mode bug

.github/workflows/ci.yml

@@ -20,7 +20,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.9'
+          - '1.10'
         os:
           - ubuntu-latest
           - macOS-latest

Project.toml

@@ -1,7 +1,7 @@
 name = "SingularIntegrals"
 uuid = "d7440221-8b5e-42fc-909c-0567823f424a"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.2.6"
+version = "0.3"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -20,18 +20,18 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
 ArrayLayouts = "1.4"
-BandedMatrices = "0.17.17, 1"
-ClassicalOrthogonalPolynomials = "0.11, 0.12"
-ContinuumArrays = "0.16, 0.17"
-FastTransforms = "0.15"
+BandedMatrices = "1"
+ClassicalOrthogonalPolynomials = "0.12"
+ContinuumArrays = "0.17"
+FastTransforms = "0.15, 0.16"
 FillArrays = "1"
 HypergeometricFunctions = "0.3.4"
 InfiniteArrays = "0.13"
-LazyArrays = "1.8"
+LazyArrays = "1.10"
 LazyBandedMatrices = "0.9"
 QuasiArrays = "0.11"
 SpecialFunctions = "1, 2"
-julia = "1.9"
+julia = "1.10"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

src/SingularIntegrals.jl

@@ -2,9 +2,9 @@ module SingularIntegrals
 using ClassicalOrthogonalPolynomials, ContinuumArrays, QuasiArrays, LazyArrays, LazyBandedMatrices, FillArrays, BandedMatrices, LinearAlgebra, SpecialFunctions, HypergeometricFunctions, InfiniteArrays
 using ContinuumArrays: @simplify, Weight, AbstractAffineQuasiVector, inbounds_getindex, broadcastbasis, MappedBasisLayouts, MemoryLayout, MappedWeightLayout, AbstractWeightLayout, ExpansionLayout, demap, basismap, AbstractBasisLayout, SubBasisLayout
 using QuasiArrays: AbstractQuasiMatrix, BroadcastQuasiMatrix, LazyQuasiArrayStyle, AbstractQuasiVecOrMat
-import ClassicalOrthogonalPolynomials: AbstractJacobiWeight, WeightedBasis, jacobimatrix, orthogonalityweight, recurrencecoefficients, _p0, Clenshaw, chop, initiateforwardrecurrence, MappedOPLayouts, unweighted
+import ClassicalOrthogonalPolynomials: AbstractJacobiWeight, WeightedBasis, jacobimatrix, orthogonalityweight, recurrencecoefficients, _p0, Clenshaw, chop, initiateforwardrecurrence, MappedOPLayouts, unweighted, WeightedOPLayout, MappedOPLayout
 using LazyBandedMatrices: Tridiagonal, SymTridiagonal, subdiagonaldata, supdiagonaldata, diagonaldata, ApplyLayout
-import LazyArrays: AbstractCachedMatrix, AbstractCachedArray, paddeddata, arguments, resizedata!, cache_filldata!, zero!, cacheddata
+import LazyArrays: AbstractCachedMatrix, AbstractCachedArray, paddeddata, arguments, resizedata!, cache_filldata!, zero!, cacheddata, LazyArrayStyle
 import Base: *, +, -, /, \, Slice, axes, getindex, sum, ==, oneto, size, broadcasted, copy, tail, view
 import LinearAlgebra: dot
 using BandedMatrices: _BandedMatrix
@@ -32,12 +32,15 @@ for Op in (:Stieltjes, :StieltjesPoint, :LogKernelPoint, :PowerKernelPoint, :Log
 end
 
 
+for lk in (:complexlogkernel, :stieltjes)
+    lk_layout = Symbol(lk, :_layout)
+    @eval $lk_layout(::AbstractBasisLayout, P, z...) = error("not implemented")
+end
 
 # general routines
 for lk in (:logkernel, :complexlogkernel, :stieltjes)
     lk_layout = Symbol(lk, :_layout)
     @eval begin
-        $lk_layout(::AbstractBasisLayout, P, z...) = error("not implemented")
         $lk_layout(::AbstractWeightLayout, w, zs::AbstractVector) = [stieltjes(w, z) for z in zs]
         function $lk_layout(::AbstractWeightLayout, w, z::Inclusion)
             axes(w,1) == z || error("Not implemented")

src/logkernel.jl

@@ -41,7 +41,8 @@ applies the log kernel log(x-t) to the columns of a quasi matrix, i.e., `(log.(x
 complexlogkernel(P, z...) = complexlogkernel_layout(MemoryLayout(P), P, z...)
 
 logkernel(wT::Weighted{T,<:ChebyshevT}) where T = ChebyshevT{T}() * Diagonal(Vcat(-convert(T,π)*log(2*one(T)),-convert(T,π)./(1:∞)))
-function logkernel_layout(::Union{MappedBasisLayouts, MappedOPLayouts}, wT::AbstractQuasiMatrix{V}) where V
+function logkernel_layout(::Union{MappedBasisLayouts, MappedOPLayouts}, wT)
+    V = eltype(wT)
     kr = basismap(wT)
     @assert kr isa AbstractAffineQuasiVector
     W = demap(wT)
@@ -53,8 +54,7 @@ function logkernel_layout(::Union{MappedBasisLayouts, MappedOPLayouts}, wT::Abst
     unweighted(wT) * Diagonal(Vcat(-convert(V,π)*(log(2*one(V))+log(abs(A)))/A,-convert(V,π) ./ (A * (1:∞))))
 end
 
-
-function logkernel_layout(::Union{MappedBasisLayouts, MappedOPLayouts}, wT::AbstractQuasiMatrix{V}, z::Number) where V
+function logkernel_demap(wT, z)
     P = demap(wT)
     kr = basismap(wT)
     z̃ = inbounds_getindex(kr, z)
@@ -65,6 +65,9 @@ function logkernel_layout(::Union{MappedBasisLayouts, MappedOPLayouts}, wT::Abst
 end
 
 
+logkernel_layout(::Union{MappedBasisLayouts, MappedOPLayouts}, wT, z::Number) = logkernel_demap(wT, z)
+logkernel_layout(::WeightedOPLayout{MappedOPLayout}, wT, z::Real) = logkernel_demap(wT, z)
+
 
 
 
@@ -74,71 +77,109 @@ end
 ####
 
 
-function logkernel(wP::Weighted{T,<:ChebyshevU}, z::Number) where T
-    if z in axes(wP,1)
-        Tn = Vcat(convert(T,π)*log(2*one(T)), convert(T,π)*ChebyshevT{T}()[z,2:end]./oneto(∞))
-        return transpose((Tn[3:end]-Tn[1:end])/2)
-    else
-        # for U_k where k>=1
-        ξ = inv(z + sqrtx2(z))
-        ζ = (convert(T,π)*ξ.^oneto(∞))./oneto(∞)
-        ζ = (ζ[3:end]- ζ[1:end])/2
+function complexlogkernel_ics(wP::Weighted{<:Any,<:ChebyshevT}, z::Number)
+    V = promote_type(eltype(wP), typeof(z))
+    ξ = inv(z + sqrtx2(z))
+    r0 = convert(V,π)*(-log(ξ)-log(convert(V,2)))
+    r1 = -convert(V,π)*ξ
+    r2 = convert(V,π)/2 + z*r1
+    r0,r1,r2
+end
+
+function complexlogkernel_recurrence(wP::Weighted{<:Any,<:ChebyshevT})
+    # We have for n ≠ 0 L_n(z) = stieltjes(U_n, z)
+    # where U_n(z) = ∫_1^x T_n(x)/sqrt(1-x^2) dx = -(1-x^2)^(-1/2)T_{n-1}(x)/n
+    # We have the 3-term recurrence
+    # T_{n+1}(x) == 2 x T_n(x) -  T_{n-1}(x)
+    # Thus
+    # U_{n+1}(x) == -(1-x^2)^(-1/2) T_n(x)/(n+1)
+    # == -(1-x^2)^(-1/2) * ( 2x T_{n-1}(x) - T_{n-2}(x))/(n+1)
+    # == -(1-x^2)^(-1/2) * ( 2n/(n+1) x T_{n-1}(x)/n -  (n-1)/(n+1) T_{n-2}(x)/(n-1))
+    # == (2n/(n+1)  * x  U_n(x) -  (n - 1)/(n+1) * U_{n-1}(x)
+
+    R = real(eltype(wP))
+    n = zero(R):∞
+    (2n) ./ (n .+ 1), Zeros{R}(∞), (n .- 1) ./ (n .+ 1)
+end
 
-        # for U_0
-        ζ = Vcat(convert(T,π)*(ξ^2/4 - (log.(abs.(ξ)) + log(2*one(T)))/2), ζ)
-        return transpose(ζ)
-    end
 
+function complexlogkernel_ics(wP::Weighted{<:Any,<:ChebyshevU}, z::Number)
+    T = promote_type(eltype(wP), typeof(z))
+    ξ = inv(z + sqrtx2(z))
+    r0 = convert(T,π)*(ξ^2/4 - (log.(abs.(ξ)) + log(2*one(T)))/2)
+    r1 = convert(T,π)*(ξ^3/3 - ξ)/2
+    r2 = convert(T,π)*(ξ^4/4 - ξ^2/2)/2
+    r0,r1,r2
 end
 
-function complexlogkernel(P::Legendre, z::Number)
+function complexlogkernel_recurrence(wP::Weighted{<:Any,<:ChebyshevU})
+    # We have for n ≠ 0 L_n(z) = stieltjes(U_n, z)
+    # where U_n(z) = ∫_1^x sqrt(1-x^2) U_n(x) dx = -(1-x^2)^(3/2)C_{n-1}(x) * 2/(n * (n + 2))
+    # where C_n(x) = C_n^{(2)}(x)
+    # We have the 3-term recurrence
+    # C_{n+1}(x) == 2(n + 2) / (n + 1) * x C_n(x) - (n + 3) / (n + 1) C_{n-1}(x)
+    # Thus
+    # U_{n+1}(x) == -(1-x^2)^(3/2) C_n(x) * 2/((n+1) * (n + 3))
+    # == -(1-x^2)^(3/2) * 2/((n+1) * (n + 3)) ( 2(n + 1) / n  * x C_{n-1}(x) - (n + 2) / n C_{n-2}(x))
+    # == -(1-x^2)^(3/2) *  (4 / (n*(n+3))  * x C_{n-1}(x) - 2 (n + 2) /(n*(n+1)*(n+3)) C_{n-2}(x))
+    # == -(1-x^2)^(3/2) *  (2 (n-1)*(n+1) / (n*(n+3))  * x  2/((n-1)*(n+1)) C_{n-1}(x) -  (n + 2)*(n-1) /(n*(n+3)) * 2/((n-1)*(n+1)) C_{n-2}(x))
+    # == (2 (n+2)/(n+3)  * x  U_n(x) -  (n + 2)*(n-1) /(n*(n+3)) * U_{n-1}(x)
+
+    R = real(eltype(wP))
+    n = zero(R):∞
+    (2*(n .+ 2)) ./ (n .+ 3), Zeros{R}(∞), (n .+ 2) .* (n .- 1) ./ (n .* (n .+ 3))
+end
+
+function complexlogkernel_ics(P::Weighted{<:Any,<:Legendre}, z::Number)
     T = promote_type(eltype(P), typeof(z))
     r0 = (1 + z)log(1 + z) - (z-1)log(z-1) - 2one(T)
     r1 = (z+1)*r0/2 + 1 - (z+1)log(z+1)
     r2 = z*r1 + 2*one(T)/3
-    transpose(RecurrenceArray(z, ((one(real(T)):2:∞)./(2:∞), Zeros{real(T)}(∞), (-one(real(T)):∞)./(2:∞)), [r0,r1,r2]))
+    r0,r1,r2
 end
 
-function complexlogkernel(P::Legendre, zs::AbstractVector)
-    T = promote_type(eltype(P), eltype(zs))
+function complexlogkernel_recurrence(wP::Weighted{<:Any,<:Legendre})
+    # We have for n ≠ 0 L_n(z) = stieltjes(U_n, z)
+    # where U_n(z) = ∫_1^x P_n(x) dx = C_{n+1}^{(-1/2)}(x)
+    # Since these are equivalent to weihted OPs (1-x^2)C_{n-1}^(3/2)(x)
+    # we know they satisfy the same recurrence coefficients.
+    # Thus the following could also be written:
+    # A,B,C = recurrencecoefficients(Ultraspherical(-1/2))
+    # A[2:end],B[2:end],C[2:end]
+    R = real(eltype(wP))
+    ((one(R):2:∞)./(2:∞), Zeros{R}(∞), (-one(R):∞)./(2:∞))
+end
+complexlogkernel(P::Legendre, z...) = complexlogkernel(Weighted(P), z...)
+logkernel(P::Legendre, x...) = logkernel(Weighted(P), x...)
+
+
+
+complexlogkernel_layout(::WeightedOPLayout, wP, z::Number) = transpose(RecurrenceArray(z, complexlogkernel_recurrence(wP), [complexlogkernel_ics(wP,z)...]))
+
+function complexlogkernel_layout(::WeightedOPLayout, wP, zs::AbstractVector)
+    T = promote_type(eltype(wP), eltype(zs))
     m = length(zs)
     data = Matrix{T}(undef, 3, m)
     for j = 1:m
         z = zs[j]
-        r0 = (1 + z)log(1 + z) - (z-1)log(z-1) - 2one(T)
-        r1 = (z+1)*r0/2 + 1 - (z+1)log(z+1)
-        r2 = z*r1 + 2*one(T)/3
-        data[1,j] = r0; data[2,j] = r1; data[3,j] = r2;
+        data[1:3,j] .= complexlogkernel_ics(wP, z)
     end
 
-    transpose(RecurrenceArray(zs, ((one(real(T)):2:∞)./(2:∞), Zeros{real(T)}(∞), (-one(real(T)):∞)./(2:∞)), data))
+    transpose(RecurrenceArray(zs, complexlogkernel_recurrence(wP), data))
 end
 
 
-logkernel(P::Legendre, z) = real.(complexlogkernel(P, complex(z)))
 
-function logkernel(P::Legendre, x::Real)
-    T = promote_type(eltype(P), typeof(x))
-    z = complex(x)
-    r0 = (1 + z)log(1 + z) - (z-1)log(z-1) - 2one(T)
-    r1 = (z+1)*r0/2 + 1 - (z+1)log(z+1)
-    r2 = z*r1 + 2*one(T)/3
-    transpose(RecurrenceArray(x, ((one(real(T)):2:∞)./(2:∞), Zeros{real(T)}(∞), (-one(real(T)):∞)./(2:∞)), [real(r0),real(r1),real(r2)]))
-end
+logkernel_layout(::AbstractBasisLayout, P, z...) = real.(complexlogkernel(P, z...))
 
-function logkernel(P::Legendre, x::AbstractVector{<:Real})
-    T = promote_type(eltype(P), eltype(x))
-    m = length(x)
-    data = Matrix{T}(undef, 3, m)
-    for j = 1:m
-        z = complex(x[j])
-        r0 = (1 + z)log(1 + z) - (z-1)log(z-1) - 2one(T)
-        r1 = (z+1)*r0/2 + 1 - (z+1)log(z+1)
-        r2 = z*r1 + 2*one(T)/3
-        data[1,j] = real(r0); data[2,j] = real(r1); data[3,j] = real(r2);
-    end
+function logkernel_layout(::WeightedOPLayout, P, x::Real)
+    L = transpose(complexlogkernel(P, complex(x)))
+    transpose(RecurrenceArray(x, (L.A, L.B, L.C), real.(L.data)))
+end
 
-    transpose(RecurrenceArray(x, ((one(real(T)):2:∞)./(2:∞), Zeros{real(T)}(∞), (-one(real(T)):∞)./(2:∞)), data))
+function logkernel_layout(::WeightedOPLayout, P, x::AbstractVector{<:Real})
+    L = transpose(complexlogkernel(P, complex(x)))
+    transpose(RecurrenceArray(x, (L.A, L.B, L.C), real.(L.data)))
 end
 
 

src/recurrence.jl

@@ -24,6 +24,9 @@ const RecurrenceMatrix{T, Z<:AbstractVector, A<:AbstractVector, B<:AbstractVecto
 function RecurrenceArray(z::Number, (A,B,C), data::AbstractVector{T}) where T
     N = length(data)
     p0, p1 = initiateforwardrecurrence(N, A, B, C, z, one(z))
+    if iszero(p1)
+        p1 = one(p1) # avoid degeneracy in recurrence. Probably needs more thought
+    end
     RecurrenceVector{T,typeof(A),typeof(B),typeof(C)}(z, A, B, C, data, size(data), T[p0], T[p1], T[])
 end
 
@@ -33,6 +36,9 @@ function RecurrenceArray(z::AbstractVector, (A,B,C), data::AbstractMatrix{T}) wh
     p1 = Vector{T}(undef, N)
     for j = axes(z,1)
         p0[j], p1[j] = initiateforwardrecurrence(M, A, B, C, z[j], one(T))
+        if iszero(p1[j])
+            p1[j] = one(p1[j]) # avoid degeneracy in recurrence. Probably needs more thought
+        end
     end
     RecurrenceMatrix{T,typeof(z),typeof(A),typeof(B),typeof(C)}(z, A, B, C, data, size(data), p0, p1, T[])
 end
@@ -165,4 +171,9 @@ _getindex_resize_iffinite!(A, kr, jr, _) = layout_getindex(A, kr, jr)
 function view(A::RecurrenceVector, kr::AbstractVector)
     resizedata!(A, maximum(kr))
     view(A.data, kr)
-end
+end
+
+###
+# broadcasted
+###
+broadcasted(::LazyArrayStyle, op, A::Transpose{<:Any,<:RecurrenceArray}) = transpose(op.(parent(A)))

test/runtests.jl

@@ -1,6 +1,6 @@
 using SingularIntegrals, ClassicalOrthogonalPolynomials, ContinuumArrays, QuasiArrays, BandedMatrices, LinearAlgebra, Test
 using SingularIntegrals: Stieltjes, StieltjesPoint, ChebyshevInterval, associated, Associated,
-        orthogonalityweight, Weighted, *, dot, LogKernelPoint
+        orthogonalityweight, Weighted, *, dot
 using LazyArrays: MemoryLayout, PaddedLayout, colsupport, rowsupport, paddeddata
 using LazyBandedMatrices: blockcolsupport, Block, BlockHcat, blockbandwidths
 

test/test_logkernel.jl

@@ -1,5 +1,5 @@
 using SingularIntegrals, ClassicalOrthogonalPolynomials, FillArrays, Test
-using SingularIntegrals: RecurrenceArray
+using SingularIntegrals: RecurrenceArray, LogKernelPoint
 using ClassicalOrthogonalPolynomials: affine
 
 @testset "ComplexLogKernelPoint" begin
@@ -29,22 +29,36 @@ end
         wU = Weighted(ChebyshevU())
         x = axes(wU,1)
         z = 0.1+0.2im
-        L = log.(abs.(z.-x'))
+        L = log.(abs.(z .- x'))
         @test L isa LogKernelPoint{Float64,ComplexF64,ComplexF64,Float64,ChebyshevInterval{Float64}}
+        @test (L * wU)[1:5] ≈ [  -1.2919202947616695, -0.20965486677056738, 0.6799687631764493, 0.13811497572177128, -0.2289481463304956]
+
+        @test L * (wU / wU \ @.(exp(x) * sqrt(1 - x^2))) ≈ -1.4812979070884382
+
+        wT = Weighted(ChebyshevT())
+        @test L * (wT / wT \ @.(exp(x) / sqrt(1 - x^2))) ≈ -1.9619040529776954 #mathematica
     end
 
     @testset "Real point" begin
+        T = ChebyshevT()
         U = ChebyshevU()
         x = axes(U,1)
 
         t = 2.0
+        @test (log.(abs.(t .- x') )* Weighted(T))[1,1:3] ≈ [1.9597591637624774, -0.8417872144769223, -0.11277810215896047] #mathematica
         @test (log.(abs.(t .- x') )* Weighted(U))[1,1:3] ≈ [1.0362686329607178,-0.4108206734393296, -0.054364775221816465] #mathematica
 
         t = 0.5
+        @test (log.(abs.(t .- x') )* Weighted(T))[1,1:3] ≈ [-2.1775860903036017, -1.5707963267948832, 0.7853981633974272] #mathematica
         @test (log.(abs.(t .- x') )* Weighted(U))[1,1:3] ≈ [-1.4814921268505252, -1.308996938995747, 0.19634954084936207] #mathematica
 
         t = 0.5+0im
+        @test (log.(abs.(t .- x') )* Weighted(T))[1,1:3] ≈ [-2.1775860903036017, -1.5707963267948832, 0.7853981633974272] #mathematica
         @test (log.(abs.(t .- x') )* Weighted(U))[1,1:3] ≈ [-1.4814921268505252, -1.308996938995747, 0.19634954084936207] #mathematica
+
+
+        # checks for degeneracy in recurrence initial conditions
+        @test iszero(logkernel(Weighted(T), 10.3)[100])
     end
 
     @testset "mapped" begin