Support norms and fix empty local minima

Project.toml

@@ -1,6 +1,6 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
-version = "0.15.13"
+version = "0.15.14"
 authors = ["Sheehan Olver <solver@mac.com>"]
 
 [deps]
@@ -51,7 +51,7 @@ IntervalSets = "0.7"
 LazyArrays = "2.8"
 LazyBandedMatrices = "0.11"
 MutableArithmetics = "1"
-QuasiArrays = "0.13"
+QuasiArrays = "0.13.2"
 RecurrenceRelationshipArrays = "0.1.2"
 RecurrenceRelationships = "0.2"
 SpecialFunctions = "1.0, 2"

src/roots.jl

@@ -10,6 +10,7 @@
 
 function colleaguematrix(P, c)
     cₙ = paddeddata(c)
+    isempty(cₙ) && return Matrix{eltype(P)}(undef, 0, 0)
     n = findlast(!iszero, cₙ)-1
     J = jacobimatrix(P)'
     C = Matrix(J[1:n,1:n])
@@ -29,15 +30,27 @@ end
 ####
 function minimum_layout(::ExpansionLayout{<:AbstractOPLayout}, f::AbstractQuasiVector, dims)
     r = findall(iszero, diff(f))
-    min(first(f), minimum(f[r]), last(f))
+    if isempty(r)
+        min(first(f), last(f))
+    else
+        min(first(f), minimum(f[r]), last(f))
+    end
 end
 
 function maximum_layout(::ExpansionLayout{<:AbstractOPLayout}, f::AbstractQuasiVector, dims)
     r = findall(iszero, diff(f))
-    max(first(f), maximum(f[r]), last(f))
+    if isempty(r)
+        max(first(f), last(f))
+    else
+        max(first(f), maximum(f[r]), last(f))
+    end
 end
 
 function extrema_layout(::ExpansionLayout{<:AbstractOPLayout}, f::AbstractQuasiVector, dims...)
     r = findall(iszero, diff(f))
-    extrema([first(f); f[r]; last(f)])
+    if isempty(r)
+        extrema([first(f); last(f)])
+    else    
+        extrema([first(f); f[r]; last(f)])
+    end
 end

test/test_legendre.jl

@@ -237,4 +237,16 @@ import QuasiArrays: MulQuasiArray
         f = P * (P \ exp.(x))
         @test P \ cumsum(f) ≈ P \ (exp.(x) .- exp(-1))
     end
+
+    @testset "norm" begin
+        f = expand(Legendre(), x -> exp(im*x))
+        @test norm(f) isa Float64
+        @test norm(f,1) isa Float64
+        @test norm(f,3) isa Float64
+        @test norm(f,Inf) isa Float64
+        @test norm(f) ≈ norm(expand(legendre(-1..1), x -> exp(im*x))) ≈ sqrt(2)
+        @test norm(f,1) ≈ 2
+        @test norm(f,Inf) ≈ 1
+        @test norm(f,3) ≈ cbrt(2) 
+    end
 end

test/test_roots.jl

@@ -37,4 +37,9 @@ end
     @test minimum(f) ≈ -2.682833127491678
     @test maximum(f) ≈ 2.6401248792053362
     @test extrema(f) == (minimum(f), maximum(f))
+
+    f = expand(ChebyshevT(), exp)
+    @test minimum(f) ≈ 1/ℯ
+    @test maximum(f) ≈ ℯ
+    @test extrema(f) == (minimum(f), maximum(f))
 end