Merge branch 'main' into dl/moveoutrecurrence

Author: dlfivefifty

Project.toml

@@ -35,7 +35,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

@@ -212,10 +212,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/chebyshev.jl

@@ -159,7 +159,7 @@ weightedgrammatrix(::ChebyshevU{V}) where V = Diagonal(Fill(convert(V,π)/2,∞)
 @simplify *(A::QuasiAdjoint{<:Any,<:Weighted{<:Any,<:ChebyshevU}}, B::ChebyshevU) = weightedgrammatrix(ChebyshevU{promote_type(eltype(A),eltype(B))}())
 
 function grammatrix(A::ChebyshevT{T}) where T
-    f = (k,j) -> isodd(j-k) ? zero(T) : -(((1 + (-1)^(j + k))*(-1 + j^2 + k^2))/(j^4 + (-1 + k^2)^2 - 2j^2*(1 + k^2)))
+    f = (k,j) -> isodd(j-k) ? zero(T) : -(T(2)*(-1 + j^2 + k^2))/((-1 + j - k)*(-1 + j + k)*(1 + j - k)*(1 + j + k))
     BroadcastMatrix{T}(f, 0:∞, (0:∞)')
 end
 
@@ -170,13 +170,13 @@ end
 
 @simplify function *(A::QuasiAdjoint{<:Any,<:ChebyshevT}, B::ChebyshevU)
     T = promote_type(eltype(A), eltype(B))
-    f = (k,j) -> isodd(j-k) ? zero(T) : ((one(T) + (-1)^(j + k))*(1 + j))/((1 + j - k)*(1 + j + k))
+    f = (k,j) -> isodd(j-k) ? zero(T) : (T(2)*(1 + j))/((1 + j - k)*(1 + j + k))
     BroadcastMatrix{T}(f, 0:∞, (0:∞)')
 end
 
 
 function grammatrix(A::Weighted{T,<:ChebyshevU}) where T
-    f = (k,j) ->  isodd(j-k) ? zero(T) : -((2*(one(T) + (-1)^(j + k))*(1 + j)*(1 + k))/((-1 + j - k)*(1 + j - k)*(1 + j + k)*(3 + j + k)))
+    f = (k,j) ->  isodd(j-k) ? zero(T) : -((T(4)*(1 + j)*(1 + k))/((-1 + j - k)*(1 + j - k)*(1 + j + k)*(3 + j + k)))
     BroadcastMatrix{T}(f, 0:∞, (0:∞)')
 end
 

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,6 +1,4 @@
 
-
-
 # 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