Skip to content

Commit 56e762b

Browse files
authored
WIP: close #261 (#270)
* close #261 * add OffsetArrays support * clean up constructors; deprecate old patterns * add test coverage
1 parent 507683d commit 56e762b

File tree

8 files changed

+110
-70
lines changed

8 files changed

+110
-70
lines changed

src/polynomials/ChebyshevT.jl

Lines changed: 9 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -26,7 +26,7 @@ ChebyshevT(1.0⋅T_0(x))
2626
struct ChebyshevT{T <: Number} <: AbstractPolynomial{T}
2727
coeffs::Vector{T}
2828
var::Symbol
29-
function ChebyshevT{T}(coeffs::AbstractVector{T}, var::Symbol) where {T <: Number}
29+
function ChebyshevT{T}(coeffs::AbstractVector{T}, var::SymbolLike=:x) where {T <: Number}
3030
length(coeffs) == 0 && return new{T}(zeros(T, 1), var)
3131
last_nz = findlast(!iszero, coeffs)
3232
last = max(1, last_nz === nothing ? 0 : last_nz)
@@ -36,6 +36,13 @@ end
3636

3737
@register ChebyshevT
3838

39+
function ChebyshevT{T}(coeffs::OffsetArray{T,1, Array{T, 1}}, var::SymbolLike=:x) where {T <: Number}
40+
cs = zeros(T, 1 + lastindex(coeffs))
41+
cs[1 .+ (firstindex(coeffs):lastindex(coeffs))] = coeffs.parent
42+
ChebyshevT{T}(cs, var)
43+
end
44+
45+
3946
function Base.convert(P::Type{<:Polynomial}, ch::ChebyshevT)
4047
if length(ch) < 3
4148
return P(ch.coeffs, ch.var)
@@ -58,6 +65,7 @@ function Base.convert(C::Type{<:ChebyshevT}, p::Polynomial)
5865
return res
5966
end
6067

68+
6169
domain(::Type{<:ChebyshevT}) = Interval(-1, 1)
6270
variable(P::Type{<:ChebyshevT}, var::SymbolLike=:x ) = P([0,1], var)
6371
"""

src/polynomials/ImmutablePolynomial.jl

Lines changed: 15 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -81,12 +81,18 @@ function ImmutablePolynomial{T}(coeffs::Tuple, var::SymbolLike=:x) where {T}
8181
ImmutablePolynomial{T}(T.(coeffs), var)
8282
end
8383

84-
# entry point from abstract.jl; note Vector -- not AbstractVector -- type
85-
function ImmutablePolynomial{T}(coeffs::Vector{T}, var::SymbolLike=:x) where {T}
84+
# entry point from abstract.jl; note T <: Number
85+
function ImmutablePolynomial{T}(coeffs::AbstractVector{T}, var::SymbolLike=:x) where {T <: Number}
8686
M = length(coeffs)
8787
ImmutablePolynomial{T}(NTuple{M,T}(tuple(coeffs...)), var)
8888
end
8989

90+
function ImmutablePolynomial{T}(coeffs::OffsetArray{T,1, Array{T, 1}}, var::SymbolLike=:x) where {T <: Number}
91+
cs = zeros(T, 1 + lastindex(coeffs))
92+
cs[1 .+ (firstindex(coeffs):lastindex(coeffs))] = coeffs.parent
93+
ImmutablePolynomial{T}(cs, var)
94+
end
95+
9096
## --
9197

9298
function ImmutablePolynomial(coeffs::Tuple, var::SymbolLike=:x)
@@ -95,21 +101,17 @@ function ImmutablePolynomial(coeffs::Tuple, var::SymbolLike=:x)
95101
ImmutablePolynomial{T}(cs, var)
96102
end
97103

98-
# # need to catch map case wihch return Variables, not Values
99-
# function ImmutablePolynomial(coeffs::Cs, var::SymbolLike=:x) where {
100-
# M, T, Cs <: Union{Values{M,T}, Variables{M,T}}}
101-
# ImmutablePolynomial{T}(Values(coeffs.v), var)
102-
# end
103-
104-
105-
106-
107104
# Convenience; pass tuple to Polynomial
108105
# Not documented, not sure this is a good idea as P(...)::P is not true...
109-
Polynomial(coeffs::NTuple{N,T}, var::SymbolLike = :x) where{N,T} =
106+
# Deprecated
107+
function Polynomial(coeffs::NTuple{N,T}, var::SymbolLike = :x) where{N,T}
108+
Base.depwarn("Use of `Polynomial(NTuple, var)` is deprecated. Use the `ImmutablePolynomial` constructor",
109+
:Polynomial)
110110
ImmutablePolynomial(coeffs, var)
111-
111+
end
112112
function Polynomial{T}(coeffs::NTuple{N,S}, var::SymbolLike = :x) where{N,T,S}
113+
Base.depwarn("Use of `Polynomial(NTuple, var)` is deprecated. Use the `ImmutablePolynomial` constructor",
114+
:Polynomial)
113115
ImmutablePolynomial{N,T}(T.(coeffs), var)
114116
end
115117

@@ -208,7 +210,6 @@ function Base.:+(p1::ImmutablePolynomial{T,N}, p2::ImmutablePolynomial{S,M}) whe
208210

209211
p1.var != p2.var && error("Polynomials must have same variable")
210212

211-
212213
if N == M
213214
cs = NTuple{N,R}(p1[i] + p2[i] for i in 0:N-1)
214215
ImmutablePolynomial{R}(cs, p1.var)
@@ -220,8 +221,6 @@ function Base.:+(p1::ImmutablePolynomial{T,N}, p2::ImmutablePolynomial{S,M}) whe
220221
ImmutablePolynomial{R,N}(cs, p1.var)
221222
end
222223

223-
224-
225224
end
226225

227226

@@ -291,9 +290,6 @@ function Base.:+(p::ImmutablePolynomial{T,N}, c::S) where {T, N, S<:Number}
291290
N == 0 && return ImmutablePolynomial{R,1}((c,), p.var)
292291
N == 1 && return ImmutablePolynomial((p[0]+c,), p.var)
293292

294-
# cs = NTuple{N,R}(iszero(i) ? p[i]+c : p[i] for i in 0:N-1)
295-
# return ImmutablePolynomial{R,N}(cs, p.var)
296-
297293
q = ImmutablePolynomial{R,1}((c,), p.var)
298294
return p + q
299295

src/polynomials/LaurentPolynomial.jl

Lines changed: 22 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,7 @@ export LaurentPolynomial
55
66
A [Laurent](https://en.wikipedia.org/wiki/Laurent_polynomial) polynomial is of the form `a_{m}x^m + ... + a_{n}x^n` where `m,n` are integers (not necessarily positive) with ` m <= n`.
77
8-
The `coeffs` specify `a_{m}, a_{m-1}, ..., a_{n}`. The range specified is of the form `m` (or `m:n`), if left empty, `m` is taken to be `0` (i.e., the coefficients refer to the standard basis). Alternatively, the coefficients can be specified using an `OffsetVector` from the `OffsetArrays` package.
8+
The `coeffs` specify `a_{m}, a_{m-1}, ..., a_{n}`. The range specified is of the form `m`, if left empty, `m` is taken to be `0` (i.e., the coefficients refer to the standard basis). Alternatively, the coefficients can be specified using an `OffsetVector` from the `OffsetArrays` package.
99
1010
Laurent polynomials and standard basis polynomials promote to Laurent polynomials. Laurent polynomials may be converted to a standard basis polynomial when `m >= 0`
1111
.
@@ -98,50 +98,42 @@ function LaurentPolynomial{T}(coeffs::AbstractVector{S}, m::Int, var::SymbolLike
9898
LaurentPolynomial{T}(T.(coeffs), m, var)
9999
end
100100

101-
function LaurentPolynomial{T}(coeffs::AbstractVector{S}, var::SymbolLike=:x) where {
102-
T <: Number, S <: Number}
103-
LaurentPolynomial{T}(T.(coeffs), 0, var)
104-
end
105-
106-
function LaurentPolynomial(coeffs::AbstractVector{T}, m::Int, var::SymbolLike=:x) where {T <: Number}
107-
LaurentPolynomial{T}(coeffs, m, Symbol(var))
101+
function LaurentPolynomial{T}(coeffs::AbstractVector{T}, var::SymbolLike=:x) where {
102+
T <: Number}
103+
LaurentPolynomial{T}(coeffs, 0, var)
108104
end
109105

110-
function LaurentPolynomial(coeffs::AbstractVector{T}, var::SymbolLike=:x) where {T <: Number}
111-
LaurentPolynomial{T}(coeffs, 0, Symbol(var))
106+
function LaurentPolynomial{T}(coeffs::OffsetArray{T, 1, Array{T,1}}, var::SymbolLike=:x) where {
107+
T<:Number}
108+
m,n = axes(coeffs, 1)
109+
LaurentPolynomial{T}(T.(coeffs.parent), m, Symbol(var))
112110
end
113111

112+
function LaurentPolynomial(coeffs::AbstractVector{T}, m::Int, var::SymbolLike=:x) where {T <: Number}
113+
LaurentPolynomial{T}(coeffs, m, Symbol(var))
114+
end
114115

115116

116-
# Add interface for OffsetArray
117-
function LaurentPolynomial{T}(coeffs::OffsetArray{S, 1, Array{S,1}}, var::SymbolLike=:x) where {
118-
T<:Number, S<:Number}
119-
m,n = axes(coeffs, 1)
120-
LaurentPolynomial{T}(T.(coeffs.parent), m:n, Symbol(var))
121-
end
122-
function LaurentPolynomial(coeffs::OffsetArray{S, 1, Array{S,1}}, var::SymbolLike=:x) where {S}
123-
LaurentPolynomial{S}(coeffs, var)
124-
end
125117

126118

127119
## Alternate with range specified
120+
## Deprecate
128121
function LaurentPolynomial{T}(coeffs::AbstractVector{S},
129122
rng::UnitRange{Int},
130123
var::Symbol=:x) where {T <: Number, S <: Number}
124+
Base.depwarn("Using a range to indicate the offset is deprecated. Use just the lower value",
125+
:LaurentPolynomial)
126+
error("")
131127
LaurentPolynomial{T}(T.(coeffs), first(rng), var)
132128
end
133129

134130
function LaurentPolynomial(coeffs::AbstractVector{T}, rng::UnitRange, var::SymbolLike=:x) where {T <: Number}
131+
Base.depwarn("Using a range to indicate the offset is deprecated. Use just the lower value",
132+
:LaurentPolynomial)
135133
LaurentPolynomial{T}(coeffs, rng, Symbol(var))
136134
end
137135

138136

139-
## Alternate interface for Polynomial
140-
Polynomial(coeffs::OffsetArray{T,1,Array{T,1}}, var::SymbolLike=:x) where {T <: Number} =
141-
LaurentPolynomial{T}(coeffs, var)
142-
143-
Polynomial{T}(coeffs::OffsetArray{S,1,Array{S,1}}, var::SymbolLike=:x) where {T <: Number, S <: Number} =
144-
LaurentPolynomial{T}(coeffs, var)
145137

146138
##
147139
## conversion
@@ -191,10 +183,10 @@ Base.:(==)(p1::LaurentPolynomial, p2::LaurentPolynomial) =
191183
Base.hash(p::LaurentPolynomial, h::UInt) = hash(p.var, hash(degreerange(p), hash(coeffs(p), h)))
192184

193185
isconstant(p::LaurentPolynomial) = iszero(lastindex(p)) && iszero(firstindex(p))
194-
basis(P::Type{<:LaurentPolynomial{T}}, n::Int, var::SymbolLike=:x) where{T} = LaurentPolynomial(ones(T,1), n:n, var)
195-
basis(P::Type{LaurentPolynomial}, n::Int, var::SymbolLike=:x) = LaurentPolynomial(ones(Float64, 1), n:n, var)
186+
basis(P::Type{<:LaurentPolynomial{T}}, n::Int, var::SymbolLike=:x) where{T} = LaurentPolynomial(ones(T,1), n, var)
187+
basis(P::Type{LaurentPolynomial}, n::Int, var::SymbolLike=:x) = LaurentPolynomial(ones(Float64, 1), n, var)
196188

197-
Base.zero(::Type{LaurentPolynomial{T}}, var=Symbollike=:x) where {T} = LaurentPolynomial{T}(zeros(T,1), 0:0, Symbol(var))
189+
Base.zero(::Type{LaurentPolynomial{T}}, var=Symbollike=:x) where {T} = LaurentPolynomial{T}(zeros(T,1), 0, Symbol(var))
198190
Base.zero(::Type{LaurentPolynomial}, var=Symbollike=:x) = zero(LaurentPolynomial{Float64}, var)
199191
Base.zero(p::P, var=Symbollike=:x) where {P <: LaurentPolynomial} = zero(P, var)
200192

@@ -545,7 +537,7 @@ function derivative(p::P, order::Integer = 1) where {T, P<:LaurentPolynomial{T}}
545537
order < 0 && error("Order of derivative must be non-negative")
546538
order == 0 && return p
547539

548-
hasnan(p) && return (P)(T[NaN], 0:0, p.var)
540+
hasnan(p) && return (P)(T[NaN], 0, p.var)
549541

550542
m,n = (extrema degreerange)(p)
551543
m = m - order
@@ -604,7 +596,7 @@ function Base.gcd(p::LaurentPolynomial{T}, q::LaurentPolynomial{T}, args...; kwa
604596
mp, Mp = (extrema degreerange)(p)
605597
mq, Mq = (extrema degreerange)(q)
606598
if mp < 0 || mq < 0
607-
throw(ArgumentError("GCD is not defined when there are `x⁻¹` terms"))
599+
throw(ArgumentError("GCD is not defined when there are `x⁻` terms"))
608600
end
609601

610602
degree(p) == 0 && return iszero(p) ? q : one(q)

src/polynomials/Polynomial.jl

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -44,6 +44,22 @@ end
4444

4545
@register Polynomial
4646

47+
48+
function Polynomial{T}(coeffs::OffsetArray{T,1,Array{T,1}}, var::SymbolLike=:x) where {T <: Number}
49+
m = firstindex(coeffs)
50+
if m < 0
51+
## depwarn
52+
Base.depwarn("Use the `LaurentPolynomial` type for offset vectors with negative first index",
53+
:Polynomial)
54+
LaurentPolynomial{T}(coeffs, var)
55+
else
56+
cs = zeros(T, 1 + lastindex(coeffs))
57+
cs[1 .+ (firstindex(coeffs):lastindex(coeffs))] = coeffs.parent
58+
Polynomial{T}(cs, var)
59+
end
60+
end
61+
62+
4763
"""
4864
(p::Polynomial)(x)
4965

src/polynomials/SparsePolynomial.jl

Lines changed: 31 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -42,41 +42,53 @@ julia> p(1)
4242
struct SparsePolynomial{T <: Number} <: StandardBasisPolynomial{T}
4343
coeffs::Dict{Int, T}
4444
var::Symbol
45-
function SparsePolynomial{T}(coeffs::Dict{Int, T}, var::Symbol) where {T <: Number}
46-
45+
function SparsePolynomial{T}(coeffs::Dict{Int, T}, var::SymbolLike) where {T <: Number}
4746
for (k,v) in coeffs
4847
iszero(v) && pop!(coeffs, k)
4948
end
50-
5149
new{T}(coeffs, var)
52-
5350
end
54-
function SparsePolynomial{T}(coeffs::AbstractVector{T}, var::Symbol) where {T <: Number}
51+
end
52+
53+
@register SparsePolynomial
5554

56-
D = Dict{Int,T}()
57-
for (i,val) in enumerate(coeffs)
58-
if !iszero(val)
59-
D[i-1] = val
60-
end
55+
function SparsePolynomial{T}(coeffs::AbstractVector{T}, var::SymbolLike=:x) where {T <: Number}
56+
D = Dict{Int,T}()
57+
for (i,val) in enumerate(coeffs)
58+
if !iszero(val)
59+
D[i-1] = val
6160
end
62-
63-
return new{T}(D, var)
64-
6561
end
62+
return SparsePolynomial{T}(D, var)
6663
end
6764

68-
@register SparsePolynomial
65+
function SparsePolynomial{T}(coeffs::OffsetArray{T,1, Array{T, 1}}, var::SymbolLike=:x) where {T <: Number}
66+
firstindex(coeffs) >= 0 || throw(ArgumentError("Use the `LaurentPolynomial` type for offset arrays with negative first index"))
67+
D = Dict{Int, T}()
68+
for i in eachindex(coeffs)
69+
D[i] = coeffs[i]
70+
end
71+
SparsePolynomial{T}(D, var)
72+
end
6973

7074
function SparsePolynomial(coeffs::Dict{Int, T}, var::SymbolLike=:x) where {T <: Number}
7175
SparsePolynomial{T}(coeffs, Symbol(var))
7276
end
73-
function SparsePolynomial(coeffs::AbstractVector{T}, var::SymbolLike=:x) where {T <: Number}
74-
SparsePolynomial{T}(coeffs, Symbol(var))
75-
end
77+
7678

7779
# Interface through `Polynomial`. As with ImmutablePolynomial, this may not be good idea...
78-
Polynomial{T}(coeffs::Dict{Int,T}, var::SymbolLike = :x) where {T} = SparsePolynomial{T}(coeffs, var)
79-
Polynomial(coeffs::Dict{Int,T}, var::SymbolLike = :x) where {T} = SparsePolynomial{T}(coeffs, var)
80+
# Deprecated
81+
function Polynomial{T}(coeffs::Dict{Int,T}, var::SymbolLike = :x) where {T}
82+
Base.depwarn("Use of `Polynomial(Dict, var)` is deprecated. Use the `SparsePolynomial` constructor",
83+
:Polynomial)
84+
SparsePolynomial{T}(coeffs, var)
85+
end
86+
87+
function Polynomial(coeffs::Dict{Int,T}, var::SymbolLike = :x) where {T}
88+
Base.depwarn("Use of `Polynomial(Dict, var)` is deprecated. Use the `SparsePolynomial` constructor",
89+
:Polynomial)
90+
SparsePolynomial{T}(coeffs, var)
91+
end
8092

8193
# conversion
8294
function Base.convert(P::Type{<:Polynomial}, q::SparsePolynomial)

test/ChebyshevT.jl

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -42,6 +42,11 @@ end
4242
p0 = ChebyshevT([0])
4343
@test iszero(p0)
4444
@test degree(p0) == -1
45+
46+
as = ones(3:4) # offsetvector
47+
bs = [0,0,0,1,1]
48+
@test ChebyshevT(as) == ChebyshevT(bs)
49+
@test ChebyshevT{Float64}(as) == ChebyshevT{Float64}(bs)
4550
end
4651

4752
@testset "Roots $i" for i in 1:5

test/StandardBasis.jl

Lines changed: 11 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,5 @@
11
using LinearAlgebra
2+
using OffsetArrays
23

34
## Test standard basis polynomials with (nearly) the same tests
45

@@ -42,7 +43,9 @@ isimmutable(::Type{<:ImmutablePolynomial}) = true
4243
@test eltype(p) == eltype(coeff)
4344
@test all([-200, -0.3, 1, 48.2] .∈ domain(p))
4445
end
46+
4547
end
48+
4649

4750
@testset "Mapdomain" begin
4851
for P in Ps
@@ -106,6 +109,13 @@ end
106109
@test degree(Polynomials.basis(P,5)) == 5
107110
@test Polynomials.isconstant(P(1))
108111
@test !Polynomials.isconstant(variable(P))
112+
113+
# OffsetVector
114+
as = ones(3:4) # offsetvector
115+
bs = [0,0,0,1,1]
116+
@test P(as) == P(bs)
117+
@test P{Float64}(as) == P{Float64}(bs)
118+
109119
end
110120
end
111121

@@ -662,7 +672,7 @@ end
662672

663673
## Issue #225 and different meanings for "conjugate"
664674
P = LaurentPolynomial
665-
p = P(rand(Complex{Float64}, 4), -1:2)
675+
p = P(rand(Complex{Float64}, 4), -1)
666676
z = rand(Complex{Float64})
667677
s = imag(z)*im
668678
@test conj(p)(z) (conj p conj)(z)

test/runtests.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@ using SpecialFunctions
66
using RecipesBase: apply_recipe
77

88
import SparseArrays: sparse, nnz
9+
using OffsetArrays
910

1011
@testset "Standard basis" begin include("StandardBasis.jl") end
1112
@testset "ChebyshevT" begin include("ChebyshevT.jl") end

0 commit comments

Comments
 (0)