Skip to content

Commit 5ce8bab

Browse files
jverzanijverzani
committed
use m (not m:n) to specify polynomial; deprecate extrema, range; close #254 
1 parent cf60052 commit 5ce8bab

File tree

2 files changed

+90
-64
lines changed

2 files changed

+90
-64
lines changed

‎src/polynomials/LaurentPolynomial.jl

Lines changed: 89 additions & 63 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:n`, if left empty, `0:length(coeffs)-1` is used (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` (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.
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
.
@@ -19,7 +19,7 @@ julia> using Polynomials
1919
julia> P = LaurentPolynomial
2020
LaurentPolynomial
2121
22-
julia> p = P([1,1,1], -1:1)
22+
julia> p = P([1,1,1], -1)
2323
LaurentPolynomial(x⁻¹ + 1 + x)
2424
2525
julia> q = P([1,1,1])
@@ -49,7 +49,7 @@ LaurentPolynomial(1.0*x + 0.5*x² + 0.3333333333333333*x³)
4949
julia> integrate(p) # x⁻¹ term is an issue
5050
ERROR: ArgumentError: Can't integrate Laurent polynomial with `x⁻¹` term
5151
52-
julia> integrate(P([1,1,1], -5:-3))
52+
julia> integrate(P([1,1,1], -5))
5353
LaurentPolynomial(-0.25*x⁻⁴ - 0.3333333333333333*x⁻³ - 0.5*x⁻²)
5454
5555
julia> x⁻¹ = inv(variable(LaurentPolynomial)) # `inv` defined on monomials
@@ -71,20 +71,18 @@ struct LaurentPolynomial{T <: Number} <: StandardBasisPolynomial{T}
7171
m::Base.RefValue{Int64}
7272
n::Base.RefValue{Int64}
7373
function LaurentPolynomial{T}(coeffs::AbstractVector{T},
74-
rng::UnitRange{Int64}=0:length(coeffs)-1,
74+
m::Int,
7575
var::Symbol=:x) where {T <: Number}
7676

77-
m,n = first(rng), last(rng)
78-
77+
7978
# trim zeros from front and back
8079
lnz = findlast(!iszero, coeffs)
8180
fnz = findfirst(!iszero, coeffs)
8281
(lnz == nothing || length(coeffs) == 0) && return new{T}(zeros(T,1), var, Ref(0), Ref(0))
83-
if lnz != length(rng) || fnz != 1
84-
coeffs = coeffs[fnz:lnz]
85-
m = m + fnz - 1
86-
n = m + (lnz-fnz)
87-
end
82+
coeffs = coeffs[fnz:lnz]
83+
m = m + fnz - 1
84+
n = m + (lnz-fnz)
85+
8886
(n-m+1 == length(coeffs)) || throw(ArgumentError("Lengths do not match"))
8987

9088
new{T}(coeffs, var, Ref(m), Ref(n))
@@ -94,6 +92,27 @@ end
9492

9593
@register LaurentPolynomial
9694

95+
## constructors
96+
function LaurentPolynomial{T}(coeffs::AbstractVector{S}, m::Int, var::Symbol=:x) where {
97+
T <: Number, S <: Number}
98+
LaurentPolynomial{T}(T.(coeffs), m, var)
99+
end
100+
101+
function LaurentPolynomial{T}(coeffs::AbstractVector{S}, var::Symbol=: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))
108+
end
109+
110+
function LaurentPolynomial(coeffs::AbstractVector{T}, var::SymbolLike=:x) where {T <: Number}
111+
LaurentPolynomial{T}(coeffs, 0, Symbol(var))
112+
end
113+
114+
115+
97116
# Add interface for OffsetArray
98117
function LaurentPolynomial{T}(coeffs::OffsetArray{S, 1, Array{S,1}}, var::SymbolLike=:x) where {T, S}
99118
m,n = axes(coeffs, 1)
@@ -104,22 +123,19 @@ function LaurentPolynomial(coeffs::OffsetArray{S, 1, Array{S,1}}, var::SymbolLi
104123
end
105124

106125

107-
126+
## Alternate with range specified
108127
function LaurentPolynomial{T}(coeffs::AbstractVector{S},
109-
rng::UnitRange{Int64}=0:length(coeffs)-1,
128+
rng::UnitRange{Int64},
110129
var::Symbol=:x) where {T <: Number, S <: Number}
111-
LaurentPolynomial{T}(T.(coeffs), rng, var)
130+
LaurentPolynomial{T}(T.(coeffs), first(rng), var)
112131
end
113132

114133
function LaurentPolynomial(coeffs::AbstractVector{T}, rng::UnitRange, var::SymbolLike=:x) where {T <: Number}
115134
LaurentPolynomial{T}(coeffs, rng, Symbol(var))
116135
end
117136

118-
function LaurentPolynomial(coeffs::AbstractVector{T}, var::SymbolLike=:x) where {T <: Number}
119-
LaurentPolynomial{T}(coeffs, 0:length(coeffs)-1, Symbol(var))
120-
end
121137

122-
## Alternate interface
138+
## Alternate interface for Polynomial
123139
Polynomial(coeffs::OffsetArray{T,1,Array{T,1}}, var::SymbolLike=:x) where {T <: Number} =
124140
LaurentPolynomial{T}(coeffs, var)
125141

@@ -138,37 +154,44 @@ Base.promote_rule(::Type{Q},::Type{P}) where {T, P <: LaurentPolynomial{T}, S, Q
138154
LaurentPolynomial{promote_type(T, S)}
139155

140156
function Base.convert(P::Type{<:Polynomial}, q::LaurentPolynomial)
141-
m,n = extrema(q)
157+
m,n = (extrema∘degreerange)(q)
142158
m < 0 && throw(ArgumentError("Can't convert a Laurent polynomial with m < 0"))
143159
P([q[i] for i in 0:n], q.var)
144160
end
145161

146162
function Base.convert(::Type{P}, q::StandardBasisPolynomial{S}) where {T, P <:LaurentPolynomial{T},S}
147163
d = degree(q)
148-
P([q[i] for i in 0:d], 0:degree(q), q.var)
164+
P([q[i] for i in 0:d], 0, q.var)
149165
end
150166

151167
##
152168
## generic functions
153169
##
154-
Base.extrema(p::LaurentPolynomial) = (p.m[],p.n[])
155-
Base.range(p::LaurentPolynomial) = p.m[]:p.n[]
170+
function Base.extrema(p::LaurentPolynomial)
171+
Base.depwarn("`extrema(::LaurentPolynomial)` is deprecated. Use `(firstindex(p), lastindex(p))`", :extrema)
172+
(p.m[], p.n[])
173+
end
174+
function Base.range(p::LaurentPolynomial)
175+
Base.depwarn("`range(::LaurentPolynomial)` is deprecated. Use `firstindex(p):lastindex(p)`", :range)
176+
p.m[]:p.n[]
177+
end
178+
156179
function Base.inv(p::LaurentPolynomial)
157-
m,n = extrema(p)
180+
m,n = degreerange(p)
158181
m != n && throw(ArgumentError("Only monomials can be inverted"))
159-
LaurentPolynomial([1/p for p in p.coeffs], -m:-n, p.var)
182+
LaurentPolynomial([1/p for p in p.coeffs], -m, p.var)
160183
end
161184

162185
##
163186
## changes to common.jl mostly as the range in the type is different
164187
##
165-
Base.copy(p::P) where {P <: LaurentPolynomial} = P(copy(coeffs(p)), range(p), p.var)
188+
Base.copy(p::P) where {P <: LaurentPolynomial} = P(copy(coeffs(p)), degreerange(p), p.var)
166189
Base.:(==)(p1::LaurentPolynomial, p2::LaurentPolynomial) =
167-
check_same_variable(p1, p2) && (range(p1) == range(p2)) && (coeffs(p1) == coeffs(p2))
168-
Base.hash(p::LaurentPolynomial, h::UInt) = hash(p.var, hash(range(p), hash(coeffs(p), h)))
190+
check_same_variable(p1, p2) && (degreerange(p1) == degreerange(p2)) && (coeffs(p1) == coeffs(p2))
191+
Base.hash(p::LaurentPolynomial, h::UInt) = hash(p.var, hash(degreerange(p), hash(coeffs(p), h)))
169192

170193
degree(p::LaurentPolynomial) = p.n[]
171-
isconstant(p::LaurentPolynomial) = range(p) == 0:0
194+
isconstant(p::LaurentPolynomial) = degreerange(p) == 0:0
172195
basis(P::Type{<:LaurentPolynomial{T}}, n::Int, var::SymbolLike=:x) where{T} = LaurentPolynomial(ones(T,1), n:n, var)
173196
basis(P::Type{LaurentPolynomial}, n::Int, var::SymbolLike=:x) = LaurentPolynomial(ones(Float64, 1), n:n, var)
174197

@@ -179,7 +202,7 @@ Base.zero(p::P, var=Symbollike=:x) where {P <: LaurentPolynomial} = zero(P, var
179202

180203
# get/set index. Work with offset
181204
function Base.getindex(p::LaurentPolynomial{T}, idx::Int) where {T <: Number}
182-
m,n = extrema(p)
205+
m,n = (extrema ∘ degreerange)(p)
183206
i = idx - m + 1
184207
(i < 1 || i > (n-m+1)) && return zero(T)
185208
p.coeffs[i]
@@ -188,7 +211,7 @@ end
188211
# extend if out of bounds
189212
function Base.setindex!(p::LaurentPolynomial{T}, value::Number, idx::Int) where {T}
190213

191-
m,n = extrema(p)
214+
m,n = (extrema ∘ degreerange)(p)
192215
if idx > n
193216
append!(p.coeffs, zeros(T, idx-n))
194217
n = idx
@@ -208,15 +231,16 @@ end
208231

209232
Base.firstindex(p::LaurentPolynomial) = p.m[]
210233
Base.lastindex(p::LaurentPolynomial) = p.n[]
211-
Base.eachindex(p::LaurentPolynomial) = range(p)
234+
Base.eachindex(p::LaurentPolynomial) = degreerange(p)
235+
degreerange(p::LaurentPolynomial) = firstindex(p):lastindex(p)
212236

213237
## chop/truncation
214238
# trim from *both* ends
215239
function chop!(p::P;
216240
rtol::Real = Base.rtoldefault(real(T)),
217241
atol::Real = 0,) where {T, P <: LaurentPolynomial{T}}
218242

219-
m0,n0 = m,n = extrema(p)
243+
m0,n0 = m,n = (extrema ∘ degreerange)(p)
220244
for k in n:-1:m
221245
if isapprox(p[k], zero(T); rtol = rtol, atol = atol)
222246
n -= 1
@@ -304,8 +328,8 @@ true
304328
"""
305329
function LinearAlgebra.conj(p::P) where {P <: LaurentPolynomial}
306330
ps = coeffs(p)
307-
m,n = extrema(p)
308-
⟒(P)(conj(ps),m:n, p.var)
331+
m = firstindex(p)
332+
⟒(P)(conj(ps), m, p.var)
309333
end
310334

311335

@@ -341,8 +365,8 @@ true
341365
function paraconj(p::LaurentPolynomial)
342366
cs = p.coeffs
343367
ds = adjoint.(cs)
344-
m,n = extrema(p)
345-
LaurentPolynomial(reverse(ds), -n:-m, p.var)
368+
n = degree(p)
369+
LaurentPolynomial(reverse(ds), -n, p.var)
346370
end
347371

348372
"""
@@ -383,13 +407,13 @@ true
383407
"""
384408
function cconj(p::LaurentPolynomial)
385409
ps = conj.(coeffs(p))
386-
m,n = extrema(p)
410+
m,n = (extrema ∘ degreerange)(p)
387411
for i in m:n
388412
if isodd(i)
389413
ps[i+1-m] *= -1
390414
end
391415
end
392-
LaurentPolynomial(ps, m:n, p.var)
416+
LaurentPolynomial(ps, m, p.var)
393417
end
394418

395419

@@ -401,7 +425,7 @@ end
401425

402426
# evaluation uses `evalpoly`
403427
function (p::LaurentPolynomial{T})(x::S) where {T,S}
404-
m,n = extrema(p)
428+
m,n = (extrema ∘ degreerange)(p)
405429
m == n == 0 && return p[0] * _one(S)
406430
if m >= 0
407431
evalpoly(x, NTuple{n+1,T}(p[i] for i in 0:n))
@@ -419,22 +443,22 @@ end
419443

420444

421445
# scalar operattoinis
422-
Base.:-(p::P) where {P <: LaurentPolynomial} = P(-coeffs(p), range(p), p.var)
446+
Base.:-(p::P) where {P <: LaurentPolynomial} = P(-coeffs(p), firstindex(p), p.var)
423447

424448
function Base.:+(p::LaurentPolynomial{T}, c::S) where {T, S <: Number}
425-
q = LaurentPolynomial([c], 0:0, p.var)
449+
q = LaurentPolynomial([c], 0, p.var)
426450
p + q
427451
end
428452

429453
function Base.:*(p::P, c::S) where {T,P <: LaurentPolynomial, S <: Number}
430454
as = c * copy(coeffs(p))
431-
return ⟒(P)(as, range(p), p.var)
455+
return ⟒(P)(as, firstindex(p), p.var)
432456
end
433457

434458

435459
function Base.:/(p::P, c::S) where {T,P <: LaurentPolynomial{T},S <: Number}
436460
R = promote_type(P, eltype(one(T) / one(S)))
437-
return R(coeffs(p) ./ c, range(p), p.var)
461+
return R(coeffs(p) ./ c, firstindex(p), p.var)
438462
end
439463

440464
##
@@ -443,25 +467,25 @@ end
443467
function Base.:+(p1::P1, p2::P2) where {T,P1<:LaurentPolynomial{T}, S, P2<:LaurentPolynomial{S}}
444468

445469
if isconstant(p1)
446-
p1 = P1(p1.coeffs, range(p1), p2.var)
470+
p1 = P1(p1.coeffs, firstindex(p1), p2.var)
447471
elseif isconstant(p2)
448-
p2 = P2(p2.coeffs, range(p2), p1.var)
472+
p2 = P2(p2.coeffs, firstindex(p2), p1.var)
449473
end
450474

451475
p1.var != p2.var && error("LaurentPolynomials must have same variable")
452476

453477
R = promote_type(T,S)
454478

455-
m1,n1 = extrema(p1)
456-
m2,n2 = extrema(p2)
479+
m1,n1 = (extrema ∘ degreerange)(p1)
480+
m2,n2 = (extrema ∘ degreerange)(p2)
457481
m,n = min(m1,m2), max(n1, n2)
458482

459483
as = zeros(R, length(m:n))
460484
for i in m:n
461485
as[1 + i-m] = p1[i] + p2[i]
462486
end
463487

464-
q = LaurentPolynomial{R}(as, m:n, p1.var)
488+
q = LaurentPolynomial{R}(as, m, p1.var)
465489
chop!(q)
466490

467491
return q
@@ -477,8 +501,8 @@ function Base.:*(p1::LaurentPolynomial{T}, p2::LaurentPolynomial{S}) where {T,S}
477501

478502
R = promote_type(T,S)
479503

480-
m1,n1 = extrema(p1)
481-
m2,n2 = extrema(p2)
504+
m1,n1 = (extrema ∘ degreerange)(p1)
505+
m2,n2 = (extrema ∘ degreerange)(p2)
482506
m,n = m1 + m2, n1+n2
483507

484508
as = zeros(R, length(m:n))
@@ -489,7 +513,7 @@ function Base.:*(p1::LaurentPolynomial{T}, p2::LaurentPolynomial{S}) where {T,S}
489513
end
490514
end
491515

492-
p = LaurentPolynomial(as, m:n, p1.var)
516+
p = LaurentPolynomial(as, m, p1.var)
493517
chop!(p)
494518

495519
return p
@@ -521,11 +545,13 @@ julia> roots(a)
521545
"""
522546
function roots(p::P; kwargs...) where {T, P <: LaurentPolynomial{T}}
523547
c = coeffs(p)
524-
r = range(p)
525-
d = r[end]-min(0,r[1])+1 # Length of the coefficient vector, taking into consideration the case when the lower degree is strictly positive (like p=3z^2).
526-
z = zeros(T,d) # Reserves space for the coefficient vector.
527-
z[max(0,r[1])+1:end] = c # Leaves the coeffs of the lower powers as zeros.
528-
a = Polynomial(z,p.var) # The root is then the root of the numerator polynomial.
548+
r = degreerange(p)
549+
d = r[end] - min(0, r[1]) + 1 # Length of the coefficient vector, taking into consideration
550+
# the case when the lower degree is strictly positive
551+
# (like p=3z^2).
552+
z = zeros(T, d) # Reserves space for the coefficient vector.
553+
z[max(0, r[1]) + 1:end] = c # Leaves the coeffs of the lower powers as zeros.
554+
a = Polynomial(z, p.var) # The root is then the root of the numerator polynomial.
529555
return roots(a; kwargs...)
530556
end
531557

@@ -539,7 +565,7 @@ function derivative(p::P, order::Integer = 1) where {T, P<:LaurentPolynomial{T}}
539565

540566
hasnan(p) && return ⟒(P)(T[NaN], 0:0, p.var)
541567

542-
m,n = extrema(p)
568+
m,n = (extrema ∘ degreerange)(p)
543569
m = m - order
544570
n = n - order
545571
as = zeros(T, length(m:n))
@@ -553,7 +579,7 @@ function derivative(p::P, order::Integer = 1) where {T, P<:LaurentPolynomial{T}}
553579
end
554580
end
555581

556-
chop!(LaurentPolynomial(as, m:n, p.var))
582+
chop!(LaurentPolynomial(as, m, p.var))
557583

558584
end
559585

@@ -564,11 +590,11 @@ function integrate(p::P, k::S) where {T, P<: LaurentPolynomial{T}, S<:Number}
564590
R = eltype((one(T)+one(S))/1)
565591

566592
if hasnan(p) || isnan(k)
567-
return P([NaN], 0:0, p.var) # not R(NaN)!! don't like XXX
593+
return P([NaN], 0, p.var) # not R(NaN)!! don't like XXX
568594
end
569595

570596

571-
m,n = extrema(p)
597+
m,n = (extrema ∘ degreerange)(p)
572598
if n < 0
573599
n = 0
574600
else
@@ -587,14 +613,14 @@ function integrate(p::P, k::S) where {T, P<: LaurentPolynomial{T}, S<:Number}
587613

588614
as[1-m] = k
589615

590-
return ⟒(P)(as, m:n, p.var)
616+
return ⟒(P)(as, m, p.var)
591617

592618
end
593619

594620

595621
function Base.gcd(p::LaurentPolynomial{T}, q::LaurentPolynomial{T}, args...; kwargs...) where {T}
596-
mp, Mp = extrema(p)
597-
mq, Mq = extrema(q)
622+
mp, Mp = (extrema ∘ degreerange)(p)
623+
mq, Mq = (extrema ∘ degreerange)(q)
598624
if mp < 0 || mq < 0
599625
throw(ArgumentError("GCD is not defined when there are `x⁻¹` terms"))
600626
end

0 commit comments

Comments
 (0)