Support expansion of matrices and complex-valued functions (#151)

* Support expansion of matrices and complex-valued functions

* Update Project.toml

* Update test_legendre.jl

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.11.3"
+version = "0.11.4"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -29,7 +29,7 @@ ArrayLayouts = "1.0.1"
 BandedMatrices = "0.17.17"
 BlockArrays = "0.16.9"
 BlockBandedMatrices = "0.12"
-ContinuumArrays = "0.15"
+ContinuumArrays = "0.15.2"
 DomainSets = "0.6"
 FFTW = "1.1"
 FastGaussQuadrature = "0.5"

src/ClassicalOrthogonalPolynomials.jl

@@ -161,7 +161,7 @@ gives the singularity structure of an expansion, e.g.,
 `JacobiWeight`.
 """
 singularities(::AbstractWeightLayout, w) = w
-singularities(lay::BroadcastLayout, a::AbstractQuasiVector) = singularitiesbroadcast(call(a), map(singularities, arguments(lay, a))...)
+singularities(lay::BroadcastLayout, a) = singularitiesbroadcast(call(a), map(singularities, arguments(lay, a))...)
 singularities(::WeightedBasisLayouts, a) = singularities(BroadcastLayout{typeof(*)}(), a)
 singularities(::WeightedOPLayout, a) = singularities(weight(a))
 singularities(w) = singularities(MemoryLayout(w), w)
@@ -170,7 +170,7 @@ singularities(::ExpansionLayout, f) = singularities(basis(f))
 singularitiesview(w, ::Inclusion) = w # for now just assume it doesn't change
 singularities(S::SubQuasiArray) = singularitiesview(singularities(parent(S)), parentindices(S)[1])
 
-basis_axes(::Inclusion{<:Any,<:AbstractInterval}, v) = basis_singularities(singularities(v))
+basis_axes(::Inclusion{<:Any,<:AbstractInterval}, v) = convert(AbstractQuasiMatrix{eltype(v)}, basis_singularities(singularities(v)))
 
 struct NoSingularities end
 

src/classical/legendre.jl

@@ -143,7 +143,7 @@ jacobimatrix(::Legendre{T}) where T =  Tridiagonal((one(T):∞)./(1:2:∞), Zero
 
 # These return vectors A[k], B[k], C[k] are from DLMF. Cause of MikaelSlevinsky we need an extra entry in C ... for now.
 function recurrencecoefficients(::Legendre{T}) where T
-    n = zero(T):∞
+    n = zero(real(T)):∞
     ((2n .+ 1) ./ (n .+ 1), Zeros{T}(∞), n ./ (n .+ 1))
 end
 

test/test_legendre.jl

@@ -189,4 +189,14 @@ import QuasiArrays: MulQuasiArray
         P̃ = Normalized(Legendre())
         @test P̃'P̃ ≡ grammatrix(P̃) ≡ Eye(∞)
     end
+
+    @testset "log sum/diff" begin
+    	P = Legendre()
+        z = 2 + im
+        f = expand(P, exp)
+        x = axes(f,1)
+        @test sum(log.(z .- x) .* f) ≈ sum(expand(Legendre{ComplexF64}(), x -> log(z - x)*exp(x)))
+        @test diff(log.(z .- x) .* P)[0.1,1:5] ≈ log.(z .- 0.1) .* diff(P)[0.1,1:5] - P[0.1,1:5] ./ (z .- 0.1)
+        @test diff(log.(z .- x) .* f)[0.1] ≈ log(z - 0.1) * exp(0.1) - exp(0.1) / (z - 0.1)
+    end
 end