Merge pull request #251 from hurak/master

Added roots() for LaurentPolynomial type.

Author: jverzani

Project.toml

@@ -2,7 +2,7 @@ name = "Polynomials"
 uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 license = "MIT"
 author = "JuliaMath"
-version = "1.1.4"
+version = "1.1.5"
 
 [deps]
 Intervals = "d8418881-c3e1-53bb-8760-2df7ec849ed5"

src/polynomials/LaurentPolynomial.jl

@@ -8,7 +8,7 @@ A [Laurent](https://en.wikipedia.org/wiki/Laurent_polynomial) polynomial is of t
 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.
 
 Laurent polynomials and standard basis polynomials  promote to  Laurent polynomials. Laurent polynomials may be  converted to a standard basis  polynomial when `m >= 0`
-. 
+.
 
 Integration will fail if there is a `x⁻¹` term in the polynomial.
 
@@ -203,7 +203,7 @@ function Base.setindex!(p::LaurentPolynomial{T}, value::Number, idx::Int) where
     p.coeffs[i] = value
 
     return p
-    
+
 end
 
 Base.firstindex(p::LaurentPolynomial) = p.m[]
@@ -215,7 +215,7 @@ Base.eachindex(p::LaurentPolynomial) = range(p)
 function chop!(p::P;
                rtol::Real = Base.rtoldefault(real(T)),
                atol::Real = 0,) where {T, P <: LaurentPolynomial{T}}
-    
+
     m0,n0 = m,n = extrema(p)
     for k in n:-1:m
         if isapprox(p[k], zero(T); rtol = rtol, atol = atol)
@@ -246,7 +246,7 @@ end
 function truncate!(p::LaurentPolynomial{T};
                   rtol::Real = Base.rtoldefault(real(T)),
                   atol::Real = 0,) where {T}
-    
+
     max_coeff = maximum(abs, coeffs(p))
     thresh = max_coeff * rtol + atol
 
@@ -259,7 +259,7 @@ function truncate!(p::LaurentPolynomial{T};
     chop!(p)
 
     return p
-    
+
 end
 
 # use unicode exponents. XXX modify printexponent to always use these?
@@ -314,7 +314,7 @@ end
 
 [cf.](https://ccrma.stanford.edu/~jos/filters/Paraunitary_FiltersC_3.html)
 
-Call `p̂ = paraconj(p)` and `p̄` = conj(p)`, then this satisfies 
+Call `p̂ = paraconj(p)` and `p̄` = conj(p)`, then this satisfies
 `conj(p(z)) = p̂(1/conj(z))` or `p̂(z) = p̄(1/z) = (conj ∘ p ∘ conj ∘ inf)(z)`.
 
 Examples:
@@ -364,7 +364,7 @@ julia> s = 2im
 julia> p = LaurentPolynomial([im,-1, -im, 1], 1:2, :s)
 LaurentPolynomial(im*s - s² - im*s³ + s⁴)
 
-julia> Polynomials.cconj(p)(s) ≈ conj(p(s)) 
+julia> Polynomials.cconj(p)(s) ≈ conj(p(s))
 true
 
 julia> a = LaurentPolynomial([-0.12, -0.29, 1],:s)
@@ -415,8 +415,8 @@ function (p::LaurentPolynomial{T})(x::S) where {T,S}
         l + r - mid
     end
 end
-                 
-        
+
+
 
 # scalar operattoinis
 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
     elseif isconstant(p2)
         p2 = P2(p2.coeffs, range(p2), p1.var)
     end
-    
+
     p1.var != p2.var && error("LaurentPolynomials must have same variable")
 
     R = promote_type(T,S)
@@ -465,7 +465,7 @@ function Base.:+(p1::P1, p2::P2) where {T,P1<:LaurentPolynomial{T}, S, P2<:Laure
     chop!(q)
 
     return q
-    
+
 end
 
 function Base.:*(p1::LaurentPolynomial{T}, p2::LaurentPolynomial{S}) where {T,S}
@@ -495,11 +495,45 @@ function Base.:*(p1::LaurentPolynomial{T}, p2::LaurentPolynomial{S}) where {T,S}
     return p
 end
 
+##
+## roots
+##
+"""
+    roots(p)
+
+Compute the roots of the Laurent polynomial `p`.
+
+The roots of a function (Laurent polynomial in this case) `a(z)` are the values of `z` for which the function vanishes. A Laurent polynomial ``a(z) = a_m z^m + a_{m+1} z^{m+1} + ... + a_{-1} z^{-1} + a_0 + a_1 z + ... + a_{n-1} z^{n-1} + a_n z^n`` can equivalently be viewed as a rational function with a multiple singularity (pole) at the origin. The roots are then the roots of the numerator polynomial. For example, ``a(z) = 1/z + 2 + z`` can be written as ``a(z) = (1+2z+z^2) / z`` and the roots of `a` are the roots of ``1+2z+z^2``.
+
+# Example
+
+```julia
+julia> p = LaurentPolynomial([24,10,-15,0,1],-2:1,:z)
+LaurentPolynomial(24*z⁻² + 10*z⁻¹ - 15 + z²)
+
+julia> roots(a)
+4-element Array{Float64,1}:
+ -3.999999999999999
+ -0.9999999999999994
+  1.9999999999999998
+  2.9999999999999982
+```
+"""
+function  roots(p::P; kwargs...)  where  {T, P <: LaurentPolynomial{T}}
+    c = coeffs(p)
+    r = range(p)
+    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).
+    z = zeros(T,d)                # Reserves space for the coefficient vector.
+    z[max(0,r[1])+1:end] = c    # Leaves the coeffs of the lower powers as zeros.
+    a = Polynomial(z,p.var)     # The root is then the root of the numerator polynomial.
+    return roots(a; kwargs...)
+end
+
 ##
 ## d/dx, ∫
 ##
 function derivative(p::P, order::Integer = 1) where {T, P<:LaurentPolynomial{T}}
-    
+
     order < 0 && error("Order of derivative must be non-negative")
     order == 0 && return p
 
@@ -518,9 +552,9 @@ function derivative(p::P, order::Integer = 1) where {T, P<:LaurentPolynomial{T}}
             as[idx] = reduce(*, (k - order + 1):k, init = p[k])
         end
     end
-    
+
     chop!(LaurentPolynomial(as, m:n, p.var))
-    
+
 end
 
 
@@ -546,15 +580,15 @@ function integrate(p::P, k::S) where {T, P<: LaurentPolynomial{T}, S<:Number}
         m = 0
     end
     as = zeros(R,  length(m:n))
-    
+
     for k in eachindex(p)
         as[1 + k+1-m]  =  p[k]/(k+1)
     end
 
     as[1-m] = k
 
     return ⟒(P)(as, m:n, p.var)
-    
+
 end
 
 
@@ -568,7 +602,7 @@ function Base.gcd(p::LaurentPolynomial{T}, q::LaurentPolynomial{T}, args...; kwa
     degree(p) == 0 && return iszero(p) ? q : one(q)
     degree(q) == 0 && return iszero(q) ? p : one(p)
     check_same_variable(p,q) || throw(ArgumentError("p and q have different symbols"))
-    
+
     pp, qq = convert(Polynomial, p), convert(Polynomial, q)
     u = gcd(pp, qq, args..., kwargs...)
     return LaurentPolynomial(coeffs(u), p.var)

test/StandardBasis.jl

@@ -247,8 +247,8 @@ end
         @test (P([NaN]) ≈ NaN)                    == (false)
         @test (P([Inf]) ≈ P([Inf]))               == ([Inf] ≈ [Inf]) # true
         @test (P([Inf]) ≈ Inf)                    == (true)
-        @test (P([1,Inf]) ≈ P([0,Inf]))           == ([1,Inf] ≈ [0,Inf]) # false 
-        @test (P([1,NaN,Inf]) ≈ P([0,NaN, Inf]))  == ([1,NaN,Inf] ≈ [0,NaN, Inf]) #false  
+        @test (P([1,Inf]) ≈ P([0,Inf]))           == ([1,Inf] ≈ [0,Inf]) # false
+        @test (P([1,NaN,Inf]) ≈ P([0,NaN, Inf]))  == ([1,NaN,Inf] ≈ [0,NaN, Inf]) #false
         @test (P([eps(), eps()]) ≈ P([0,0]))      == ([eps(), eps()] ≈ [0,0]) # false
         @test (P([1,eps(), 1]) ≈ P([1,0,1]))      == ([1,eps(), 1] ≈ [1,0,1]) # true
         @test (P([1,2]) ≈ P([1,2,eps()]))         == ([1,2,0] ≈ [1,2,eps()])
@@ -431,6 +431,8 @@ end
         @test roots(pNULL) == []
         @test sort(roots(pR)) == [1 // 2, 3 // 2]
 
+        @test sort(roots(LaurentPolynomial([24,10,-15,0,1],-2:1,:z)))≈[-4.0,-1.0,2.0,3.0]
+
         A = [1 0; 0 1]
         @test fromroots(A) == Polynomial(Float64[1, -2, 1])
         p = fromroots(P, A)
@@ -511,7 +513,7 @@ end
         q   = [3, p1]
         if P !=  ImmutablePolynomial
             @test q isa Vector{typeof(p1)}
-            @test p isa Vector{typeof(p2)} 
+            @test p isa Vector{typeof(p2)}
         else
             @test q isa Vector{P{eltype(p1)}} # ImmutablePolynomial{Int64,N} where {N}, different  Ns
             @test p isa Vector{P{eltype(p2)}} # ImmutablePolynomial{Int64,N} where {N}, different  Ns
@@ -651,7 +653,7 @@ end
         end
 
         @test eltype(p1) == Int
-        for P in Ps 
+        for P in Ps
             p1 = P([1,2,0,3])
             @test eltype(collect(p1)) <: P{Int}
             @test eltype(collect(P{Float64}, p1)) <: P{Float64}
@@ -754,7 +756,7 @@ end
         for n in (2,5,10,20,50, 100)
             p = (x-1)^n * (x-2)^n * (x-3)
             q = (x-1) * (x-2) * (x-4)
-            @test degree(gcd(p,q, method=:numerical)) == 2  
+            @test degree(gcd(p,q, method=:numerical)) == 2
         end
     end
 end
@@ -889,6 +891,6 @@ end
             @test typeof(p₁) == P{T,5} # conversion works
             @test_throws InexactError  P{T}(rand(5))
         end
-    end        
+    end
 
 end