Assume NoSingularities by default

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <solver@mac.com>"]
-version = "0.13.5"
+version = "0.13.6"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"

src/ClassicalOrthogonalPolynomials.jl

@@ -211,7 +211,10 @@ singularitiesview(w, ::Inclusion) = w # for now just assume it doesn't change
 singularitiesview(w, ind) = view(w, ind)
 singularities(S::SubQuasiArray) = singularitiesview(singularities(parent(S)), parentindices(S)[1])
 
-basis_axes(::Inclusion{<:Any,<:AbstractInterval}, v) = convert(AbstractQuasiMatrix{eltype(v)}, basis_singularities(singularities(v)))
+
+basis_singularities(ax, ::NoSingularities) = basis(ax)
+basis_singularities(ax, sing) = basis_singularities(sing)
+basis_axes(ax::Inclusion{<:Any,<:AbstractInterval}, v) = convert(AbstractQuasiMatrix{eltype(v)}, basis_singularities(ax, singularities(v)))
 
 struct NoSingularities end
 

src/classical/jacobi.jl

@@ -54,6 +54,9 @@ hasboundedendpoints(w::AbstractJacobiWeight) = w.a ≥ 0 && w.b ≥ 0
 singularities(a::AbstractAffineQuasiVector) = singularities(a.x)
 
 
+## default is to just assume no singularities
+singularitiesbroadcast(_...) = NoSingularities()
+
 for op in (:+, :*)
     @eval singularitiesbroadcast(::typeof($op), A, B, C, D...) = singularitiesbroadcast(*, singularitiesbroadcast(*, A, B), C, D...)
 end
@@ -66,7 +69,6 @@ _parentindices(a, b...) = parentindices(a)
 # for singularitiesbroadcast(literal_pow), ^, ...)
 singularitiesbroadcast(F::Function, G::Function, V::SubQuasiArray, K) = singularitiesbroadcast(F, G, parent(V), K)[parentindices(V)...]
 singularitiesbroadcast(F, V::Union{NoSingularities,SubQuasiArray}...) = singularitiesbroadcast(F, map(_parent,V)...)[_parentindices(V...)...]
-singularitiesbroadcast(::typeof(*), V::Union{NoSingularities,SubQuasiArray}...) = singularitiesbroadcast(*, map(_parent,V)...)[_parentindices(V...)...]
 
 
 abstract type AbstractJacobi{T} <: OrthogonalPolynomial{T} end

src/classical/legendre.jl

@@ -37,6 +37,7 @@ for op in (:+, :-, :*)
 end
 singularitiesbroadcast(::typeof(^), L::LegendreWeight, ::NoSingularities) = L
 singularitiesbroadcast(::typeof(/), ::NoSingularities, L::LegendreWeight) = L # can't find roots
+singularitiesbroadcast(::typeof(*), V::Union{NoSingularities,SubQuasiArray,LegendreWeight}...) = singularitiesbroadcast(*, map(_parent,V)...)[_parentindices(V...)...]
 
 singularities(::AbstractJacobi{T}) where T = LegendreWeight{T}()
 singularities(::Inclusion{T,<:ChebyshevInterval}) where T = LegendreWeight{T}()

test/test_ultraspherical.jl

@@ -183,4 +183,12 @@ using ClassicalOrthogonalPolynomials: grammatrix
         L = Weighted(Jacobi(1,0)) \ Weighted(Ultraspherical(3/2))
         @test L[1:10,1:10] == (Weighted(Jacobi(1,0)) \ Weighted(Jacobi(1,1)))[1:10,1:10] * (Jacobi(1,1) \ Ultraspherical(3/2))[1:10,1:10]
     end
+
+    @testset "expand w/ no singularities" begin
+        z = 2
+        P = Legendre()
+        C = Ultraspherical(-1/2)
+        x = axes(P,1)
+        @test sum(@.(ultrasphericalc(1, -1/2, x)/(z-x))) ≈ sum(C[:,2] ./ (z .- x))
+    end
 end