Skip to content

Commit d11546a

Browse files
committed
Added roots() for LaurentPolynomial type.
1 parent ad67b84 commit d11546a

File tree

1 file changed

+40
-17
lines changed

1 file changed

+40
-17
lines changed

src/polynomials/LaurentPolynomial.jl

Lines changed: 40 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@ A [Laurent](https://en.wikipedia.org/wiki/Laurent_polynomial) polynomial is of t
88
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.
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`
11-
.
11+
.
1212
1313
Integration will fail if there is a `x⁻¹` term in the polynomial.
1414
@@ -203,7 +203,7 @@ function Base.setindex!(p::LaurentPolynomial{T}, value::Number, idx::Int) where
203203
p.coeffs[i] = value
204204

205205
return p
206-
206+
207207
end
208208

209209
Base.firstindex(p::LaurentPolynomial) = p.m[]
@@ -215,7 +215,7 @@ Base.eachindex(p::LaurentPolynomial) = range(p)
215215
function chop!(p::P;
216216
rtol::Real = Base.rtoldefault(real(T)),
217217
atol::Real = 0,) where {T, P <: LaurentPolynomial{T}}
218-
218+
219219
m0,n0 = m,n = extrema(p)
220220
for k in n:-1:m
221221
if isapprox(p[k], zero(T); rtol = rtol, atol = atol)
@@ -246,7 +246,7 @@ end
246246
function truncate!(p::LaurentPolynomial{T};
247247
rtol::Real = Base.rtoldefault(real(T)),
248248
atol::Real = 0,) where {T}
249-
249+
250250
max_coeff = maximum(abs, coeffs(p))
251251
thresh = max_coeff * rtol + atol
252252

@@ -259,7 +259,7 @@ function truncate!(p::LaurentPolynomial{T};
259259
chop!(p)
260260

261261
return p
262-
262+
263263
end
264264

265265
# use unicode exponents. XXX modify printexponent to always use these?
@@ -314,7 +314,7 @@ end
314314
315315
[cf.](https://ccrma.stanford.edu/~jos/filters/Paraunitary_FiltersC_3.html)
316316
317-
Call `p̂ = paraconj(p)` and `p̄` = conj(p)`, then this satisfies
317+
Call `p̂ = paraconj(p)` and `p̄` = conj(p)`, then this satisfies
318318
`conj(p(z)) = p̂(1/conj(z))` or `p̂(z) = p̄(1/z) = (conj ∘ p ∘ conj ∘ inf)(z)`.
319319
320320
Examples:
@@ -364,7 +364,7 @@ julia> s = 2im
364364
julia> p = LaurentPolynomial([im,-1, -im, 1], 1:2, :s)
365365
LaurentPolynomial(im*s - s² - im*s³ + s⁴)
366366
367-
julia> Polynomials.cconj(p)(s) ≈ conj(p(s))
367+
julia> Polynomials.cconj(p)(s) ≈ conj(p(s))
368368
true
369369
370370
julia> a = LaurentPolynomial([-0.12, -0.29, 1],:s)
@@ -415,8 +415,8 @@ function (p::LaurentPolynomial{T})(x::S) where {T,S}
415415
l + r - mid
416416
end
417417
end
418-
419-
418+
419+
420420

421421
# scalar operattoinis
422422
Base.:-(p::P) where {P <: LaurentPolynomial} = P(-coeffs(p), range(p), p.var)
@@ -447,7 +447,7 @@ function Base.:+(p1::P1, p2::P2) where {T,P1<:LaurentPolynomial{T}, S, P2<:Laure
447447
elseif isconstant(p2)
448448
p2 = P2(p2.coeffs, range(p2), p1.var)
449449
end
450-
450+
451451
p1.var != p2.var && error("LaurentPolynomials must have same variable")
452452

453453
R = promote_type(T,S)
@@ -465,7 +465,7 @@ function Base.:+(p1::P1, p2::P2) where {T,P1<:LaurentPolynomial{T}, S, P2<:Laure
465465
chop!(q)
466466

467467
return q
468-
468+
469469
end
470470

471471
function Base.:*(p1::LaurentPolynomial{T}, p2::LaurentPolynomial{S}) where {T,S}
@@ -495,11 +495,34 @@ function Base.:*(p1::LaurentPolynomial{T}, p2::LaurentPolynomial{S}) where {T,S}
495495
return p
496496
end
497497

498+
##
499+
## roots
500+
##
501+
"""
502+
roots(p)
503+
504+
Compute the roots of the Laurent polynomial `p`.
505+
506+
# Example
507+
508+
```julia
509+
julia> p = LaurentPolynomial([24,10,-15,0,1],-2:1,:z)
510+
LaurentPolynomial(24*z⁻² + 10*z⁻¹ - 15 + z²)
511+
512+
513+
```
514+
"""
515+
function roots(p::P; kwargs...) where {T, P <: LaurentPolynomial{T}}
516+
c = coeffs(p)
517+
a = Polynomial(c,p.var)
518+
return roots(p; kwargs...)
519+
end
520+
498521
##
499522
## d/dx, ∫
500523
##
501524
function derivative(p::P, order::Integer = 1) where {T, P<:LaurentPolynomial{T}}
502-
525+
503526
order < 0 && error("Order of derivative must be non-negative")
504527
order == 0 && return p
505528

@@ -518,9 +541,9 @@ function derivative(p::P, order::Integer = 1) where {T, P<:LaurentPolynomial{T}}
518541
as[idx] = reduce(*, (k - order + 1):k, init = p[k])
519542
end
520543
end
521-
544+
522545
chop!(LaurentPolynomial(as, m:n, p.var))
523-
546+
524547
end
525548

526549

@@ -546,15 +569,15 @@ function integrate(p::P, k::S) where {T, P<: LaurentPolynomial{T}, S<:Number}
546569
m = 0
547570
end
548571
as = zeros(R, length(m:n))
549-
572+
550573
for k in eachindex(p)
551574
as[1 + k+1-m] = p[k]/(k+1)
552575
end
553576

554577
as[1-m] = k
555578

556579
return (P)(as, m:n, p.var)
557-
580+
558581
end
559582

560583

@@ -568,7 +591,7 @@ function Base.gcd(p::LaurentPolynomial{T}, q::LaurentPolynomial{T}, args...; kwa
568591
degree(p) == 0 && return iszero(p) ? q : one(q)
569592
degree(q) == 0 && return iszero(q) ? p : one(p)
570593
check_same_variable(p,q) || throw(ArgumentError("p and q have different symbols"))
571-
594+
572595
pp, qq = convert(Polynomial, p), convert(Polynomial, q)
573596
u = gcd(pp, qq, args..., kwargs...)
574597
return LaurentPolynomial(coeffs(u), p.var)

0 commit comments

Comments
 (0)