1+ # ###
2+ # Associated
3+ # ###
4+
5+
6+ """
7+ AssociatedWeighted(P)
8+
9+ We normalise so that `orthogonalityweight(::Associated)` is a probability measure.
10+ """
11+ struct AssociatedWeight{T,OPs<: AbstractQuasiMatrix{T} } <: Weight{T}
12+ P:: OPs
13+ end
14+ axes (w:: AssociatedWeight ) = (axes (w. P,1 ),)
15+
16+ sum (:: AssociatedWeight{T} ) where T = one (T)
17+
18+ """
19+ Associated(P)
20+
21+ constructs the associated orthogonal polynomials for P, which have the Jacobi matrix
22+
23+ jacobimatrix(P)[2:end,2:end]
24+
25+ and constant first term. Or alternatively
26+
27+ w = orthogonalityweight(P)
28+ A = recurrencecoefficients(P)[1]
29+ Associated(P) == (w/(sum(w)*A[1]))'*((P[:,2:end]' - P[:,2:end]) ./ (x' - x))
30+
31+ where `x = axes(P,1)`.
32+ """
33+
34+ struct Associated{T, OPs<: AbstractQuasiMatrix{T} } <: OrthogonalPolynomial{T}
35+ P:: OPs
36+ end
37+
38+ associated (P) = Associated (P)
39+
40+ axes (Q:: Associated ) = axes (Q. P)
41+ == (A:: Associated , B:: Associated ) = A. P == B. P
42+
43+ orthogonalityweight (Q:: Associated ) = AssociatedWeight (Q. P)
44+
45+ function associated_jacobimatrix (X:: Tridiagonal )
46+ c,a,b = subdiagonaldata (X),diagonaldata (X),supdiagonaldata (X)
47+ Tridiagonal (c[2 : end ], a[2 : end ], b[2 : end ])
48+ end
49+
50+ function associated_jacobimatrix (X:: SymTridiagonal )
51+ a,b = diagonaldata (X),supdiagonaldata (X)
52+ SymTridiagonal (a[2 : end ], b[2 : end ])
53+ end
54+ jacobimatrix (a:: Associated ) = associated_jacobimatrix (jacobimatrix (a. P))
55+
56+ associated (:: ChebyshevT{T} ) where T = ChebyshevU {T} ()
57+ associated (:: ChebyshevU{T} ) where T = ChebyshevU {T} ()
58+
59+
160const StieltjesPoint{T,V,D} = BroadcastQuasiMatrix{T,typeof (inv),Tuple{BroadcastQuasiMatrix{T,typeof (- ),Tuple{T,QuasiAdjoint{V,Inclusion{V,D}}}}}}
261const ConvKernel{T,D} = BroadcastQuasiMatrix{T,typeof (- ),Tuple{D,QuasiAdjoint{T,D}}}
362const Hilbert{T,D} = BroadcastQuasiMatrix{T,typeof (inv),Tuple{ConvKernel{T,Inclusion{T,D}}}}
463const LogKernel{T,D} = BroadcastQuasiMatrix{T,typeof (log),Tuple{BroadcastQuasiMatrix{T,typeof (abs),Tuple{ConvKernel{T,Inclusion{T,D}}}}}}
564const PowKernel{T,D,F<: Real } = BroadcastQuasiMatrix{T,typeof (^ ),Tuple{BroadcastQuasiMatrix{T,typeof (abs),Tuple{ConvKernel{T,Inclusion{T,D}}}},F}}
665
766
8- @simplify function * (S:: StieltjesPoint{<:Any,<:Any,<:ChebyshevInterval} , wT:: WeightedBasis{<:Any,<:ChebyshevTWeight,<:ChebyshevT} )
9- w,T = wT. args
10- J = jacobimatrix (T)
11- z, x = parent (S). args[1 ]. args
12- transpose ((J' - z* I) \ [- π; zeros (∞)])
67+ @simplify function * (H:: Hilbert , w:: ChebyshevTWeight )
68+ T = promote_type (eltype (H), eltype (w))
69+ zeros (T, axes (w,1 ))
70+ end
71+
72+ @simplify function * (H:: Hilbert , w:: ChebyshevUWeight )
73+ T = promote_type (eltype (H), eltype (w))
74+ fill (convert (T,π), axes (w,1 ))
1375end
1476
15- @simplify function * (H:: Hilbert{<:Any,<:ChebyshevInterval} , wT:: WeightedBasis{<:Any,<:ChebyshevTWeight,<:ChebyshevT} )
77+ @simplify function * (H:: Hilbert , w:: LegendreWeight )
78+ T = promote_type (eltype (H), eltype (w))
79+ x = axes (w,1 )
80+ log .(x .+ 1 ) .- log .(1 .- x)
81+ end
82+
83+ @simplify function * (H:: Hilbert , wT:: Weighted{<:Any,<:ChebyshevT} )
1684 T = promote_type (eltype (H), eltype (wT))
17- ApplyQuasiArray ( * , ChebyshevU {T} (), _BandedMatrix (Fill (- convert (T,π),1 ,∞), ℵ₀, - 1 , 1 ) )
85+ ChebyshevU {T} () * _BandedMatrix (Fill (- convert (T,π),1 ,∞), ℵ₀, - 1 , 1 )
1886end
1987
20- @simplify function * (H:: Hilbert{<:Any,<:ChebyshevInterval} , wU:: WeightedBasis {<:Any,<:ChebyshevUWeight ,<:ChebyshevU}
88+ @simplify function * (H:: Hilbert , wU:: Weighted {<:Any,<:ChebyshevU}
2189 T = promote_type (eltype (H), eltype (wU))
22- ApplyQuasiArray (* , ChebyshevT {T} (), _BandedMatrix (Fill (convert (T,π),1 ,∞), ℵ₀, 1 , - 1 ))
90+ ChebyshevT {T} () * _BandedMatrix (Fill (convert (T,π),1 ,∞), ℵ₀, 1 , - 1 )
91+ end
92+
93+
94+ @simplify function * (H:: Hilbert , wP:: Weighted{<:Any,<:OrthogonalPolynomial} )
95+ P = wP. P
96+ w = orthogonalityweight (P)
97+ A = recurrencecoefficients (P)[1 ]
98+ (- A[1 ]* sum (w))* [zero (axes (P,1 )) associated (P)] + (H* w) .* P
2399end
24100
101+ @simplify * (H:: Hilbert , P:: Legendre ) = H * Weighted (P)
102+
25103# ##
26104# LogKernel
27105# ##
28106
29- @simplify function * (L:: LogKernel{<:Any,<:ChebyshevInterval} , wT:: WeightedBasis {<:Any,<:ChebyshevTWeight ,<:ChebyshevT}
107+ @simplify function * (L:: LogKernel , wT:: Weighted {<:Any,<:ChebyshevT}
30108 T = promote_type (eltype (L), eltype (wT))
31- ApplyQuasiArray ( * , ChebyshevT {T} (), Diagonal (Vcat (- π* log (2 * one (T)),- convert (T,π)./ (1 : ∞) )))
109+ ChebyshevT {T} () * Diagonal (Vcat (- π* log (2 * one (T)),- convert (T,π)./ (1 : ∞)))
32110end
33111
34112
113+
35114# ##
36115# PowKernel
37116# ##
38117
39- @simplify function * (K:: PowKernel{<:Any,<:ChebyshevInterval} , wT:: WeightedBasis {<:Any,<:JacobiWeight ,<:Jacobi}
118+ @simplify function * (K:: PowKernel , wT:: Weighted {<:Any,<:Jacobi}
40119 T = promote_type (eltype (K), eltype (wT))
41120 cnv,α = K. args
42121 x,y = K. args[1 ]. args[1 ]. args
43122 @assert x' == y
44123 β = (- α- 1 )/ 2
45- ApplyQuasiArray (* , ChebyshevT {T} (), Diagonal (1 : ∞))
124+ error (" Not implemented"
125+ # ChebyshevT{T}() * Diagonal(1:∞)
46126end
47127
48128
49129# ###
50130# StieltjesPoint
51131# ###
52132
53- @simplify function * (S:: StieltjesPoint , wT:: SubQuasiArray{<:Any,2,<:WeightedBasis,<:Tuple{<:AbstractAffineQuasiVector,<:Any}} )
133+ stieltjesmoment_jacobi_normalization (n:: Int ,α:: Real ,β:: Real ) = 2 ^ (α+ β)* gamma (n+ α+ 1 )* gamma (n+ β+ 1 )/ gamma (2 n+ α+ β+ 2 )
134+
135+ @simplify function * (S:: StieltjesPoint , w:: AbstractJacobiWeight )
136+ α,β = w. a,w. b
137+ z,_ = parent (S). args[1 ]. args
138+ (x = 2 / (1 - z);stieltjesmoment_jacobi_normalization (0 ,α,β)* HypergeometricFunctions. mxa_₂F₁ (1 ,α+ 1 ,α+ β+ 2 ,x))
139+ end
140+
141+ @simplify function * (S:: StieltjesPoint , wP:: Weighted )
142+ P = wP. P
143+ w = orthogonalityweight (P)
144+ X = jacobimatrix (P)
145+ z, x = parent (S). args[1 ]. args
146+ transpose ((X' - z* I) \ [- sum (w)* _p0 (P); zeros (∞)])
147+ end
148+
149+
150+ @simplify function * (S:: StieltjesPoint , wT:: SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Any}} )
54151 P = parent (wT)
55152 z, x = parent (S). args[1 ]. args
56153 z̃ = inbounds_getindex (parentindices (wT)[1 ], z)
57154 x̃ = axes (P,1 )
58155 (inv .(z̃ .- x̃' ) * P)[:,parentindices (wT)[2 ]]
59156end
60157
61- @simplify function * (H:: Hilbert , wT:: SubQuasiArray{<:Any,2,<:WeightedBasis ,<:Tuple{<:AbstractAffineQuasiVector,<:Any}} )
158+ @simplify function * (H:: Hilbert , wT:: SubQuasiArray{<:Any,2,<:Any ,<:Tuple{<:AbstractAffineQuasiVector,<:Any}} )
62159 P = parent (wT)
63160 x = axes (P,1 )
64161 apply (* , inv .(x .- x' ), P)[parentindices (wT)... ]
65162end
66163
67164
68- @simplify function * (L:: LogKernel , wT:: SubQuasiArray{<:Any,2,<:WeightedBasis ,<:Tuple{<:AbstractAffineQuasiVector,<:Slice}} )
165+ @simplify function * (L:: LogKernel , wT:: SubQuasiArray{<:Any,2,<:Any ,<:Tuple{<:AbstractAffineQuasiVector,<:Slice}} )
69166 V = promote_type (eltype (L), eltype (wT))
70- P = parent (wT)
167+ wP = parent (wT)
71168 kr, jr = parentindices (wT)
72- @assert P isa WeightedBasis{<: Any ,<: ChebyshevWeight ,<: Chebyshev }
73- x = axes (P,1 )
74- w,T = P. args
75- D = T \ apply (* , log .(abs .(x .- x' )), P)
169+ x = axes (wP,1 )
170+ w,T = arguments (ApplyLayout {typeof(*)} (), wP)
171+ @assert w isa ChebyshevTWeight
172+ @assert T isa ChebyshevT
173+ D = T \ (log .(abs .(x .- x' )) * wP)
76174 c = inv (2 * kr. A)
77175 T[kr,:] * Diagonal (Vcat (2 * convert (V,π)* c* log (c), 2 c* D. diag. args[2 ]))
78- end
176+ end
177+
178+
179+
180+ # ## generic fallback
181+ for Op in (:Hilbert , :StieltjesPoint , :LogKernel , :PowKernel )
182+ @eval @simplify function * (H:: $Op , wP:: WeightedBasis{<:Any,<:Weight,<:Any} )
183+ w,P = wP. args
184+ Q = OrthogonalPolynomial (w)
185+ (H * Weighted (Q)) * (Q \ P)
186+ end
187+ end
0 commit comments