Assume NoSingularities by default in expand (#208)

* Assume NoSingularities by default

* Update ClassicalOrthogonalPolynomials.jl

* Remove overloads

* Update ci.yml

Author: dlfivefifty

.github/workflows/ci.yml

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

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.13.5"
+version = "0.13.6"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -33,7 +33,7 @@ ContinuumArrays = "0.18.3"
 DomainSets = "0.6, 0.7"
 FFTW = "1.1"
 FastGaussQuadrature = "1"
-FastTransforms = "0.15.15, 0.16"
+FastTransforms = "0.16.6"
 FillArrays = "1"
 HypergeometricFunctions = "0.3.4"
 InfiniteArrays = " 0.14"

src/ClassicalOrthogonalPolynomials.jl

@@ -211,10 +211,13 @@ 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)))
-
 struct NoSingularities end
 
+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)))
+
+
 singularities(::Number) = NoSingularities()
 singularities(r::Base.RefValue) = r[] # pass through
 

src/classical/jacobi.jl

@@ -54,10 +54,15 @@ 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
 
+singularitiesbroadcast(::typeof(*), V::Union{NoSingularities,SubQuasiArray}...) = singularitiesbroadcast(*, map(_parent,V)...)[_parentindices(V...)...]
+
 
 _parent(::NoSingularities) = NoSingularities()
 _parent(a) = parent(a)
@@ -66,7 +71,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/clenshaw.jl

@@ -1,13 +1,4 @@
 
-##
-# For Chebyshev T. Note the shift in indexing is fine due to the AbstractFill
-##
-Base.@propagate_inbounds _forwardrecurrence_next(n, A::Vcat{<:Any,1,<:Tuple{<:Number,<:AbstractFill}}, B::Zeros, C::Ones, x, p0, p1) = 
-    _forwardrecurrence_next(n, A.args[2], B, C, x, p0, p1)
-
-Base.@propagate_inbounds _clenshaw_next(n, A::Vcat{<:Any,1,<:Tuple{<:Number,<:AbstractFill}}, B::Zeros, C::Ones, x, c, bn1, bn2) = 
-    _clenshaw_next(n, A.args[2], B, C, x, c, bn1, bn2)
-
 # Assume 1 normalization
 _p0(A) = one(eltype(A))
 

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