Skip to content

Commit 496aabf

Browse files
authored
Generalize the polynomial constructors to accept AbstractVectors as lists of coefficients (#303)
* Generalize the StandardBasisPolynomial constrctors to accept AbstractVectors * Fix test on v1.0 * Fix doctest meta * Fix ChebyshevT
1 parent 2e7cbb4 commit 496aabf

File tree

8 files changed

+114
-109
lines changed

8 files changed

+114
-109
lines changed

.github/workflows/ci.yml

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -60,8 +60,9 @@ jobs:
6060
Pkg.instantiate()'
6161
- run: |
6262
julia --project=docs -e '
63-
using Documenter: doctest
63+
using Documenter: doctest, DocMeta
6464
using Polynomials
65+
DocMeta.setdocmeta!(Polynomials, :DocTestSetup, :(using Polynomials); recursive = true)
6566
doctest(Polynomials)'
6667
- run: julia --project=docs docs/make.jl
6768
env:

src/polynomials/ChebyshevT.jl

Lines changed: 13 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,16 @@
11
export ChebyshevT
22

33
"""
4-
ChebyshevT{<:Number}(coeffs::AbstractVector, var=:x)
4+
ChebyshevT{<:Number}(coeffs::AbstractVector, [var = :x])
55
66
Chebyshev polynomial of the first kind.
77
8-
Construct a polynomial from its coefficients `a`, lowest order first, optionally in
9-
terms of the given variable `x`. `x` can be a character, symbol, or string.
8+
Construct a polynomial from its coefficients `coeffs`, lowest order first, optionally in
9+
terms of the given variable `var`, which can be a character, symbol, or string.
10+
11+
!!! note
12+
`ChebyshevT` is not axis-aware, and it treats `coeffs` simply as a list of coefficients with the first
13+
index always corresponding to the coefficient of `T_0(x)`.
1014
1115
# Examples
1216
@@ -28,21 +32,18 @@ struct ChebyshevT{T <: Number} <: AbstractPolynomial{T}
2832
var::Symbol
2933
function ChebyshevT{T}(coeffs::AbstractVector{T}, var::SymbolLike=:x) where {T <: Number}
3034
length(coeffs) == 0 && return new{T}(zeros(T, 1), var)
31-
last_nz = findlast(!iszero, coeffs)
35+
if Base.has_offset_axes(coeffs)
36+
@warn "ignoring the axis offset of the coefficient vector"
37+
end
38+
c = OffsetArrays.no_offset_view(coeffs)
39+
last_nz = findlast(!iszero, c)
3240
last = max(1, last_nz === nothing ? 0 : last_nz)
33-
return new{T}(coeffs[1:last], var)
41+
return new{T}(c[1:last], var)
3442
end
3543
end
3644

3745
@register ChebyshevT
3846

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-
4647
function Base.convert(P::Type{<:Polynomial}, ch::ChebyshevT)
4748
if length(ch) < 3
4849
return P(ch.coeffs, ch.var)

src/polynomials/ImmutablePolynomial.jl

Lines changed: 10 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -22,8 +22,11 @@ As the coefficient size is a compile-time constant, several performance
2222
improvements are possible. For example, immutable polynomials can take advantage of
2323
faster polynomial evaluation provided by `evalpoly` from Julia 1.4.
2424
25+
!!! note
26+
`ImmutablePolynomial` is not axis-aware, and it treats `coeffs` simply as a list of coefficients with the first
27+
index always corresponding to the constant term.
2528
26-
# Examples
29+
# Examples
2730
2831
```jldoctest
2932
julia> using Polynomials
@@ -63,6 +66,9 @@ function ImmutablePolynomial{T,N}(coeffs::Tuple, var::SymbolLike=:x) where {T,N
6366
end
6467

6568
function ImmutablePolynomial{T,N}(coeffs::AbstractVector{S}, var::SymbolLike=:x) where {T <: Number, N, S}
69+
if Base.has_offset_axes(coeffs)
70+
@warn "ignoring the axis offset of the coefficient vector"
71+
end
6672
ImmutablePolynomial{T,N}(NTuple{N,T}(tuple(coeffs...)), var)
6773
end
6874

@@ -83,16 +89,13 @@ end
8389

8490
# entry point from abstract.jl; note T <: Number
8591
function ImmutablePolynomial{T}(coeffs::AbstractVector{T}, var::SymbolLike=:x) where {T <: Number}
92+
if Base.has_offset_axes(coeffs)
93+
@warn "ignoring the axes of the coefficient vector and treating it as a list"
94+
end
8695
M = length(coeffs)
8796
ImmutablePolynomial{T}(NTuple{M,T}(tuple(coeffs...)), var)
8897
end
8998

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-
9699
## --
97100

98101
function ImmutablePolynomial(coeffs::Tuple, var::SymbolLike=:x)

src/polynomials/LaurentPolynomial.jl

Lines changed: 22 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -1,18 +1,25 @@
11
export LaurentPolynomial
22

33
"""
4-
LaurentPolynomial(coeffs, range, var)
4+
LaurentPolynomial(coeffs::AbstractVector, [m::Integer = 0], [var = :x])
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`, 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}`.
9+
The argument `m` represents the lowest exponent of the variable in the series, and is taken to be zero by default.
910
10-
Laurent polynomials and standard basis polynomials promote to Laurent polynomials. Laurent polynomials may be converted to a standard basis polynomial when `m >= 0`
11+
Laurent polynomials and standard basis polynomials promote to Laurent polynomials. Laurent polynomials may be converted to a standard basis polynomial when `m >= 0`
1112
.
1213
1314
Integration will fail if there is a `x⁻¹` term in the polynomial.
1415
15-
Example:
16+
!!! note
17+
`LaurentPolynomial` is not axis-aware by default, and it treats `coeffs` simply as a
18+
list of coefficients with the first index always corresponding to the constant term.
19+
In order to use the axis of `coeffs` as the exponents of the variable `var`,
20+
set `m` to `firstindex(coeff)` in the constructor.
21+
22+
# Examples:
1623
```jldoctest laurent
1724
julia> using Polynomials
1825
@@ -74,19 +81,22 @@ struct LaurentPolynomial{T <: Number} <: StandardBasisPolynomial{T}
7481
m::Int,
7582
var::Symbol=:x) where {T <: Number}
7683

77-
84+
if Base.has_offset_axes(coeffs)
85+
@warn "ignoring the axis offset of the coefficient vector"
86+
end
87+
c = OffsetArrays.no_offset_view(coeffs) # ensure 1-based indexing
7888
# trim zeros from front and back
79-
lnz = findlast(!iszero, coeffs)
80-
fnz = findfirst(!iszero, coeffs)
81-
(lnz == nothing || length(coeffs) == 0) && return new{T}(zeros(T,1), var, Ref(0), Ref(0))
82-
coeffs = coeffs[fnz:lnz]
89+
lnz = findlast(!iszero, c)
90+
fnz = findfirst(!iszero, c)
91+
(lnz == nothing || length(c) == 0) && return new{T}(zeros(T,1), var, Ref(0), Ref(0))
92+
c = c[fnz:lnz]
93+
8394
m = m + fnz - 1
8495
n = m + (lnz-fnz)
8596

86-
(n-m+1 == length(coeffs)) || throw(ArgumentError("Lengths do not match"))
87-
88-
new{T}(coeffs, var, Ref(m), Ref(n))
97+
(n-m+1 == length(c)) || throw(ArgumentError("Lengths do not match"))
8998

99+
new{T}(c, var, Ref(m), Ref(n))
90100
end
91101
end
92102

@@ -103,19 +113,10 @@ function LaurentPolynomial{T}(coeffs::AbstractVector{T}, var::SymbolLike=:x) whe
103113
LaurentPolynomial{T}(coeffs, 0, var)
104114
end
105115

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))
110-
end
111-
112116
function LaurentPolynomial(coeffs::AbstractVector{T}, m::Int, var::SymbolLike=:x) where {T <: Number}
113117
LaurentPolynomial{T}(coeffs, m, Symbol(var))
114118
end
115119

116-
117-
118-
119120
## Alternate with range specified
120121
## Deprecate
121122
function LaurentPolynomial{T}(coeffs::AbstractVector{S},

src/polynomials/Polynomial.jl

Lines changed: 15 additions & 40 deletions
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,10 @@
11
export Polynomial
22

33
"""
4-
Polynomial{T<:Number}(coeffs::AbstractVector{T}, var=:x)
4+
Polynomial{T<:Number}(coeffs::AbstractVector{T}, [var = :x])
55
6-
Construct a polynomial from its coefficients `a`, lowest order first, optionally in
7-
terms of the given variable `x`. `x` can be a character, symbol, or string.
6+
Construct a polynomial from its coefficients `coeffs`, lowest order first, optionally in
7+
terms of the given variable `var` which may be a character, symbol, or a string.
88
99
If ``p = a_n x^n + \\ldots + a_2 x^2 + a_1 x + a_0``, we construct this through
1010
`Polynomial([a_0, a_1, ..., a_n])`.
@@ -13,13 +13,12 @@ The usual arithmetic operators are overloaded to work with polynomials as well a
1313
with combinations of polynomials and scalars. However, operations involving two
1414
polynomials of different variables causes an error except those involving a constant polynomial.
1515
16-
# Examples
17-
```@meta
18-
DocTestSetup = quote
19-
using Polynomials
20-
end
21-
```
16+
!!! note
17+
`Polynomial` is not axis-aware, and it treats `coeffs` simply as a list of coefficients with the first
18+
index always corresponding to the constant term. In order to use the axis of `coeffs` as exponents,
19+
consider using a [`LaurentPolynomial`](@ref) or possibly a [`SparsePolynomial`](@ref).
2220
21+
# Examples
2322
```jldoctest
2423
julia> Polynomial([1, 0, 3, 4])
2524
Polynomial(1 + 3*x^2 + 4*x^3)
@@ -34,49 +33,25 @@ Polynomial(1.0)
3433
struct Polynomial{T <: Number} <: StandardBasisPolynomial{T}
3534
coeffs::Vector{T}
3635
var::Symbol
37-
function Polynomial{T}(coeffs::Vector{T}, var::Symbol) where {T <: Number}
36+
function Polynomial{T}(coeffs::AbstractVector{T}, var::Symbol) where {T <: Number}
37+
if Base.has_offset_axes(coeffs)
38+
@warn "ignoring the axis offset of the coefficient vector"
39+
end
3840
length(coeffs) == 0 && return new{T}(zeros(T, 1), var)
39-
last_nz = findlast(!iszero, coeffs)
41+
c = OffsetArrays.no_offset_view(coeffs) # ensure 1-based indexing
42+
last_nz = findlast(!iszero, c)
4043
last = max(1, last_nz === nothing ? 0 : last_nz)
41-
return new{T}(coeffs[1:last], var)
44+
return new{T}(c[1:last], var)
4245
end
4346
end
4447

4548
@register Polynomial
4649

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-
6350
"""
6451
(p::Polynomial)(x)
6552
6653
Evaluate the polynomial using [Horner's Method](https://en.wikipedia.org/wiki/Horner%27s_method), also known as synthetic division, as implemented in `evalpoly` of base `Julia`.
6754
68-
```@meta
69-
DocTestSetup = quote
70-
using Polynomials
71-
end
72-
```
73-
74-
```@meta
75-
DocTestSetup = quote
76-
using Polynomials
77-
end
78-
```
79-
8055
# Examples
8156
```jldoctest
8257
julia> p = Polynomial([1, 0, 3])

src/polynomials/SparsePolynomial.jl

Lines changed: 14 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -1,14 +1,14 @@
11
export SparsePolynomial
22

33
"""
4-
SparsePolynomial(coeffs::Dict, var)
4+
SparsePolynomial(coeffs::Dict, [var = :x])
55
66
Polynomials in the standard basis backed by a dictionary holding the
77
non-zero coefficients. For polynomials of high degree, this might be
88
advantageous. Addition and multiplication with constant polynomials
99
are treated as having no symbol.
1010
11-
Examples:
11+
# Examples:
1212
1313
```jldoctest
1414
julia> using Polynomials
@@ -42,37 +42,30 @@ 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::SymbolLike) where {T <: Number}
45+
function SparsePolynomial{T}(coeffs::AbstractDict{Int, T}, var::SymbolLike) where {T <: Number}
46+
c = Dict(coeffs)
4647
for (k,v) in coeffs
47-
iszero(v) && pop!(coeffs, k)
48+
iszero(v) && pop!(c, k)
4849
end
49-
new{T}(coeffs, var)
50+
new{T}(c, Symbol(var))
5051
end
5152
end
5253

5354
@register SparsePolynomial
5455

5556
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
60-
end
61-
end
62-
return SparsePolynomial{T}(D, var)
63-
end
57+
firstindex(coeffs) >= 0 || throw(ArgumentError("Use the `LaurentPolynomial` type for arrays with a negative first index"))
6458

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]
59+
if Base.has_offset_axes(coeffs)
60+
@warn "ignoring the axis offset of the coefficient vector"
7061
end
71-
SparsePolynomial{T}(D, var)
62+
c = OffsetArrays.no_offset_view(coeffs) # ensure 1-based indexing
63+
p = Dict{Int,T}(i - 1 => v for (i,v) in pairs(c))
64+
return SparsePolynomial{T}(p, var)
7265
end
7366

74-
function SparsePolynomial(coeffs::Dict{Int, T}, var::SymbolLike=:x) where {T <: Number}
75-
SparsePolynomial{T}(coeffs, Symbol(var))
67+
function SparsePolynomial(coeffs::AbstractDict{Int, T}, var::SymbolLike=:x) where {T <: Number}
68+
SparsePolynomial{T}(coeffs, var)
7669
end
7770

7871

test/ChebyshevT.jl

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -43,10 +43,14 @@ end
4343
@test iszero(p0)
4444
@test degree(p0) == -1
4545

46-
as = ones(3:4) # offsetvector
47-
bs = [0,0,0,1,1]
46+
as = ones(3:4)
47+
bs = parent(as)
4848
@test ChebyshevT(as) == ChebyshevT(bs)
4949
@test ChebyshevT{Float64}(as) == ChebyshevT{Float64}(bs)
50+
51+
a = [1,1]
52+
b = OffsetVector(a, axes(a))
53+
@test ChebyshevT(a) == ChebyshevT(b)
5054
end
5155

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

0 commit comments

Comments
 (0)