Legendre log kernel (#14)

* start Legendre log kernel

* add ComplexLogKernelPoint

* Update test_hilbert.jl

* Update test_logkernel.jl

* fix tests

Project.toml

@@ -19,7 +19,7 @@ QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
-ArrayLayouts = "1"
+ArrayLayouts = "1.0.5"
 BandedMatrices = "0.17.17"
 ClassicalOrthogonalPolynomials = "0.8, 0.9"
 ContinuumArrays = "0.12.6"

src/logkernel.jl

@@ -1,3 +1,4 @@
+const ComplexLogKernelPoint{T,C,W<:Number,V,D} = BroadcastQuasiMatrix{T,typeof(log),Tuple{ConvKernel{C,W,V,D}}}
 const LogKernelPoint{T<:Real,C,W<:Number,V,D} = BroadcastQuasiMatrix{T,typeof(log),Tuple{BroadcastQuasiMatrix{T,typeof(abs),Tuple{ConvKernel{C,W,V,D}}}}}
 const LogKernel{T,D1,D2} = BroadcastQuasiMatrix{T,typeof(log),Tuple{BroadcastQuasiMatrix{T,typeof(abs),Tuple{ConvKernel{T,Inclusion{T,D1},T,D2}}}}}
 
@@ -47,6 +48,26 @@ end
 
 end
 
+@simplify function *(L::ComplexLogKernelPoint, P::Legendre)
+    T = promote_type(eltype(L), eltype(P))
+    z, xc = parent(L).args[1].args
+    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]))
+end
+
+@simplify function *(L::LogKernelPoint, P::Legendre)
+    T = promote_type(eltype(L), eltype(P))
+    z, x = parent(L).args[1].args[1].args
+    real.(log.(complex(z) .- x) * P)
+end
+
+
+###
+# Maps
+###
+
 
 @simplify function *(L::LogKernelPoint, wT::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Any}})
     P = parent(wT)

test/test_hilbert.jl

@@ -1,4 +1,5 @@
 using SingularIntegrals, ClassicalOrthogonalPolynomials, Test
+using LazyBandedMatrices: blockcolsupport, Block, BlockHcat, blockbandwidths, paddeddata
 
 @testset "Hilbert" begin
     @testset "weights" begin

test/test_logkernel.jl

@@ -1,3 +1,25 @@
+using SingularIntegrals, ClassicalOrthogonalPolynomials, FillArrays, Test
+import SingularIntegrals: RecurrenceArray
+
+@testset "ComplexLogKernelPoint" begin
+
+    P = Legendre()
+    Pc = Legendre{ComplexF64}()
+    x = axes(P,1)
+
+
+    L = (z,k) -> sum(Pc / Pc \ (log.(z .- x).*P[:,k+1]))
+
+    for z in (5, 1+2im, -1+2im, 1-2im, -3+0.0im, -3-0.0im)
+        @test (log.(z .- x') * P)[1:5] ≈ L.(z, 0:4)
+    end
+
+    for z in (-5,-1,0,0.1)
+        @test_throws DomainError log.(z .- x') * P
+    end
+
+end
+
 @testset "LogKernelPoint" begin
     @testset "Complex point" begin
         wU = Weighted(ChebyshevU())
@@ -35,6 +57,44 @@
         f = wU / wU \ @.(sqrt(2-x)sqrt(x-1)exp(x))
         @test L*f ≈ 2.2374312398976586 # MAthematica
     end
+
+    @testset "Legendre" begin
+        z = 5
+        P = Legendre()
+        x = axes(P,1)
+
+        L = (z,k) -> sum(P / P \ (log.(abs.(z .- x)).*P[:,k+1]))
+
+        for z in (5, 1+2im, -1+2im, 1-2im, -3+0.0im, -3-0.0im, -3)
+            @test (log.(abs.(z .- x')) * P)[1:5] ≈ L.(z, 0:4)
+        end
+
+        @testset "derivation" begin
+            W = Weighted(Jacobi(1,1))
+
+            z = 5
+
+            @test L(z,0) ≈ 2log(z+1) + inv.(z .- x') * (P / P \ (x .- 1)) ≈ (1 + z)log(1 + z) - (z-1)log(z-1) - 2
+            @test L(z,1) ≈ (inv.(z .- x') * W)[1]/(-2)
+
+            @test z * L(z,0) ≈ 2z*log(z+1) + z*(inv.(z .- x') * (P / P \ (x .- 1))) ≈
+                        2z*log(z+1) + (inv.(z .- x') * (P / P \ (x.^2 .- 1))) + (inv.(z .- x') * (P / P \ (1 .- x))) - 2 ≈
+                        2L(z,1) - L(z,0) + 2(z+1)*log(z+1) - 2
+            @test z * L(z,1) ≈ z * (inv.(z .- x') * W)[1]/(-2) ≈ -2/3 + (inv.(z .- x') * (x .* W))[1]/(-2) ≈
+                        -2/3 + (inv.(z .- x') * W)[2]/(-4) ≈ -2/3 + L(z,2)
+
+            for k = 2:5
+                @test z * L(z,k) ≈ (k-1)/(2k+1)*L(z,k-1)+ (k+2)/(2k+1)*L(z,k+1)
+                @test (2k+1)/(k+2)*z * L(z,k) ≈ (k-1)/(k+2)*L(z,k-1)+ L(z,k+1)
+            end
+
+            r0 = (1 + z)log(1 + z) - (z-1)log(z-1) - 2
+            r1 = (z+1)*r0/2 + 1 - (z+1)log(z+1)
+            r2 = z*r1 + 2*one(z)/3
+            r = RecurrenceArray(z, ((1:2:∞)./(2:∞), Zeros(∞), (-1:∞)./(2:∞)), [r0,r1,r2])
+            @test r[1:10] ≈ L.(z,0:9)
+        end
+    end
 end