ChebyshevGrid in points (#216)

* ChebyshevGrid in points

* scaled Chebyshev grid

* Add tests

Project.toml

@@ -1,6 +1,6 @@
 name = "ApproxFunOrthogonalPolynomials"
 uuid = "b70543e2-c0d9-56b8-a290-0d4d6d4de211"
-version = "0.6.17"
+version = "0.6.18"
 
 [deps]
 ApproxFunBase = "fbd15aa5-315a-5a7d-a8a4-24992e37be05"

src/ApproxFunOrthogonalPolynomials.jl

@@ -82,7 +82,30 @@ using StaticArrays: SVector
 import LinearAlgebra: isdiag
 
 points(d::IntervalOrSegmentDomain{T},n::Integer) where {T} =
-    fromcanonical.(Ref(d), chebyshevpoints(float(real(eltype(T))), n))  # eltype to handle point
+    _maybefromcanonical(d, chebyshevpoints(float(real(eltype(T))), n))  # eltype to handle point
+_maybefromcanonical(d, pts) = fromcanonical.(Ref(d), pts)
+_maybefromcanonical(::ChebyshevInterval, pts::FastTransforms.ChebyshevGrid) = pts
+function _maybefromcanonical(d::IntervalOrSegment{<:Union{Number, SVector}}, pts::FastTransforms.ChebyshevGrid)
+    shift = mean(d)
+    scale = complexlength(d) / 2
+    ShiftedChebyshevGrid(pts, shift, scale)
+end
+
+struct ShiftedChebyshevGrid{T, S, C<:FastTransforms.ChebyshevGrid} <: AbstractVector{T}
+    grid :: C
+    shift :: S
+    scale :: S
+end
+function ShiftedChebyshevGrid(grid::G, shift::S, scale::S) where {G,S}
+    T = typeof(zero(eltype(G)) * zero(S))
+    ShiftedChebyshevGrid{T,S,G}(grid, shift, scale)
+end
+Base.size(S::ShiftedChebyshevGrid) = size(S.grid)
+Base.@propagate_inbounds Base.getindex(S::ShiftedChebyshevGrid, i::Int) = S.shift + S.grid[i] * S.scale
+function Base.showarg(io::IO, S::ShiftedChebyshevGrid, toplevel::Bool)
+    print(io, "ShiftedChebyshevGrid{", eltype(S), "}")
+end
+
 bary(v::AbstractVector{Float64},d::IntervalOrSegmentDomain,x::Float64) = bary(v,tocanonical(d,x))
 
 strictconvert(T::Type, x) = convert(T, x)::T

src/roots.jl

@@ -217,10 +217,18 @@ function rootsunit_coeffs(c::Vector{T}, htol::Float64,clplan::ClenshawPlan{S,T})
         # Find the roots of the polynomial on each piece and then concatenate. Recurse if necessary.
 
         # Evaluate the polynomial at Chebyshev grids on both intervals:
-        #(clenshaw! overwrites points, which only makes sence if c is real)
+        #(clenshaw! overwrites points, which only makes sense if c is real)
 
-        v1 = isa(c,Vector{Float64}) ? clenshaw!( c, points(-1..splitPoint,n),clplan) : clenshaw( c, points(-1..splitPoint,n),clplan)
-        v2 = isa(c,Vector{Float64}) ? clenshaw!( c, points(splitPoint..1 ,n),clplan) : clenshaw( c, points(splitPoint..1 ,n),clplan)
+        v1 = if isa(c, Vector{Float64})
+            clenshaw!(c, convert(Vector, points(-1..splitPoint,n)), clplan)
+        else
+            clenshaw(c, points(-1..splitPoint,n),clplan)
+        end
+        v2 = if isa(c, Vector{Float64})
+            clenshaw!(c, convert(Vector, points(splitPoint..1 ,n)), clplan)
+        else
+            clenshaw(c, points(splitPoint..1 ,n),clplan)
+        end
 
         # Recurse (and map roots back to original interval):
         p = plan_chebyshevtransform( v1 )

test/ChebyshevTest.jl

@@ -39,6 +39,12 @@ using ApproxFunOrthogonalPolynomials: forwardrecurrence
         @test y == 2Fun()
     end
 
+    @testset "points" begin
+        p = points(Chebyshev(-1..1), 4)
+        @test points(Chebyshev(), 4) ≈ p
+        @test contains(summary(p), "ShiftedChebyshevGrid{Float64}")
+    end
+
     @testset "inference in Space(::Interval)" begin
         S = @inferred Space(0..1)
         f = Fun(S)