Merge branch 'main' into patch-4

Author: putianyi889

.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 <solver@mac.com>"]
-version = "0.13.5"
+version = "0.13.7"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -22,6 +22,8 @@ LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
 LazyBandedMatrices = "d7e5e226-e90b-4449-9968-0f923699bf6f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
+RecurrenceRelationships = "807425ed-42ea-44d6-a357-6771516d7b2c"
+RecurrenceRelationshipArrays = "b889d2dc-af3c-4820-88a8-238fa91d3518"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
@@ -33,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"
@@ -42,6 +44,8 @@ IntervalSets = "0.7"
 LazyArrays = "2.2"
 LazyBandedMatrices = "0.10"
 QuasiArrays = "0.11"
+RecurrenceRelationships = "0.1"
+RecurrenceRelationshipArrays = "0.1"
 SpecialFunctions = "1.0, 2"
 julia = "1.10"
 

docs/src/index.md

@@ -180,9 +180,6 @@ ClassicalOrthogonalPolynomials.normalizationconstant
 ClassicalOrthogonalPolynomials.OrthogonalPolynomialRatio
 ```
 ```@docs
-ClassicalOrthogonalPolynomials.Clenshaw
-```
-```@docs
 ClassicalOrthogonalPolynomials.singularities
 ```
 ```@docs

src/ClassicalOrthogonalPolynomials.jl

@@ -5,7 +5,7 @@ using IntervalSets: UnitRange
 using ContinuumArrays, QuasiArrays, LazyArrays, FillArrays, BandedMatrices, BlockArrays,
     IntervalSets, DomainSets, ArrayLayouts, SpecialFunctions,
     InfiniteLinearAlgebra, InfiniteArrays, LinearAlgebra, FastGaussQuadrature, FastTransforms, FFTW,
-    LazyBandedMatrices, HypergeometricFunctions
+    LazyBandedMatrices, HypergeometricFunctions, RecurrenceRelationships
 
 import Base: @_inline_meta, axes, getindex, unsafe_getindex, convert, prod, *, /, \, +, -,
                 IndexStyle, IndexLinear, ==, OneTo, tail, similar, copyto!, copy, setindex,
@@ -18,7 +18,7 @@ import LazyArrays: MemoryLayout, Applied, ApplyStyle, flatten, _flatten, adjoint
                 _mul_arguments, CachedVector, CachedMatrix, LazyVector, LazyMatrix, axpy!, AbstractLazyLayout, BroadcastLayout,
                 AbstractCachedVector, AbstractCachedMatrix, paddeddata, cache_filldata!,
                 simplifiable, PaddedArray, converteltype, simplify
-import ArrayLayouts: MatMulVecAdd, materialize!, _fill_lmul!, sublayout, sub_materialize, lmul!, ldiv!, ldiv, transposelayout, triangulardata,
+import ArrayLayouts: MatMulVecAdd, materialize!, sublayout, sub_materialize, lmul!, ldiv!, ldiv, transposelayout, triangulardata,
                         subdiagonaldata, diagonaldata, supdiagonaldata, mul, rowsupport, colsupport
 import LazyBandedMatrices: SymTridiagonal, Bidiagonal, Tridiagonal, unitblocks, BlockRange1, AbstractLazyBandedLayout
 import LinearAlgebra: pinv, factorize, qr, adjoint, transpose, dot, mul!, reflectorApply!
@@ -39,9 +39,10 @@ import ContinuumArrays: Basis, Weight, basis_axes, @simplify, Identity, Abstract
     AffineQuasiVector, AffineMap, AbstractWeightLayout, AbstractWeightedBasisLayout, WeightedBasisLayout, WeightedBasisLayouts, demap, AbstractBasisLayout, BasisLayout,
     checkpoints, weight, unweighted, MappedBasisLayouts, sum_layout, invmap, plan_ldiv, layout_broadcasted, MappedBasisLayout, SubBasisLayout, broadcastbasis_layout,
     plan_grid_transform, plan_transform, MAX_PLOT_POINTS, MulPlan, grammatrix, AdjointBasisLayout, grammatrix_layout, plan_transform_layout, _cumsum
-import FastTransforms: Λ, forwardrecurrence, forwardrecurrence!, _forwardrecurrence!, clenshaw, clenshaw!,
-                        _forwardrecurrence_next, _clenshaw_next, check_clenshaw_recurrences, ChebyshevGrid, chebyshevpoints, Plan, ScaledPlan, th_cheb2leg
-
+import FastTransforms: Λ, ChebyshevGrid, chebyshevpoints, Plan, ScaledPlan, th_cheb2leg, pochhammer
+import RecurrenceRelationships: forwardrecurrence, forwardrecurrence!, clenshaw, clenshaw!,
+                        check_clenshaw_recurrences
+import RecurrenceRelationshipArrays: initiateforwardrecurrence, Clenshaw
 import FastGaussQuadrature: jacobimoment
 
 import BlockArrays: blockedrange, _BlockedUnitRange, unblock, _BlockArray, block, blockindex, BlockSlice, blockvec
@@ -211,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

@@ -194,17 +194,16 @@ end
 ##########
 
 # Ultraspherical(1)\(D*Chebyshev())
-function diff(S::ChebyshevT{T}; dims=1) where T
+function diff(::ChebyshevT{T}; dims=1) where T
     D = _BandedMatrix((zero(T):∞)', ℵ₀, -1,1)
     ApplyQuasiMatrix(*, ChebyshevU{T}(), D)
 end
 
-function diff(W::Weighted{T,<:ChebyshevU}; dims=1) where T
+function diff(::Weighted{T,<:ChebyshevU}; dims=1) where T
     D =  _BandedMatrix((-one(T):-one(T):(-∞))', ℵ₀, 1,-1)
     ApplyQuasiMatrix(*, Weighted(ChebyshevT{T}()), D)
 end
 
-
 #####
 # Conversion
 #####

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
@@ -384,6 +388,11 @@ end
 # Jacobi(a+1,b+1)\(D*Jacobi(a,b))
 diff(S::Jacobi; dims=1) = ApplyQuasiMatrix(*, Jacobi(S.a+1,S.b+1), _BandedMatrix((((1:∞) .+ (S.a + S.b))/2)', ℵ₀, -1,1))
 
+function diff(S::Jacobi{T}, m::Integer; dims=1) where T
+    D = _BandedMatrix((pochhammer.((S.a + S.b+1):∞, m)/convert(T, 2)^m)', ℵ₀, -m, m)
+    ApplyQuasiMatrix(*, Jacobi(S.a+m,S.b+m), D)
+end
+
 
 #L_6^t
 function diff(WS::HalfWeighted{:a,<:Any,<:Jacobi}; dims=1)

src/classical/ultraspherical.jl

@@ -119,9 +119,12 @@ end
 ##########
 
 # Ultraspherical(1)\(D*Chebyshev())
-diff(S::ChebyshevU; dims=1) = diff(Ultraspherical(S))
+diff(S::ChebyshevU, m...; dims=1) = diff(Ultraspherical(S), m...; dims)
+diff(S::Legendre, m...; dims=1) = diff(Ultraspherical(S), m...; dims)
+
 
 # Ultraspherical(1/2)\(D*Legendre())
+# Special cased as its a Ones
 function diff(S::Legendre{T}; dims=1) where T
     A = _BandedMatrix(Ones{T}(1,∞), ℵ₀, -1,1)
     ApplyQuasiMatrix(*, Ultraspherical{T}(convert(T,3)/2), A)
@@ -134,6 +137,20 @@ function diff(S::Ultraspherical{T}; dims=1) where T
     ApplyQuasiMatrix(*, Ultraspherical{T}(S.λ+1), A)
 end
 
+# higher order 
+
+function diff(::ChebyshevT{T}, m::Integer; dims=1) where T
+    μ = pochhammer(one(T),m-1)*convert(T,2)^(m-1)
+    D = _BandedMatrix((μ * (0:∞))', ℵ₀, -m, m)
+    ApplyQuasiMatrix(*, Ultraspherical{T}(m), D)
+end
+
+function diff(C::Ultraspherical{T}, m::Integer; dims=1) where T
+    μ = pochhammer(convert(T,C.λ),m)*convert(T,2)^m
+    D = _BandedMatrix(Fill(μ,1,∞), ℵ₀, -m, m)
+    ApplyQuasiMatrix(*, Ultraspherical{T}(C.λ+m), D)
+end
+
 # Ultraspherical(λ-1)\ (D*wUltraspherical(λ))
 function diff(WS::Weighted{T,<:Ultraspherical}; dims=1) where T
     S = WS.P