From de521312d549f5876afdd30d59cf25e818d99cc3 Mon Sep 17 00:00:00 2001 From: Neven Sajko Date: Wed, 12 Oct 2022 21:36:05 +0200 Subject: [PATCH] prevent run-time dispatch in roots(::AbstractVector) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit findlast does not necessarily return an Integer, it may also return Nothing, depending on the value of its arguments. The latter case needs to be handled to prevent run-time dispatch for the slicing operation. Additionally, the truncated coefficients slice is now a view instead of a copy. This should make things slightly more performant, apart from making the code simpler. JET warnings before this commit: ``` julia> @report_opt PolynomialRoots.roots(rand(5)) ═════ 1 possible error found ═════ ┌ @ /home/nsajko/.julia/packages/PolynomialRoots/1iZQh/src/PolynomialRoots.jl:610 PolynomialRoots.:(var"#roots#2")(PolynomialRoots.NaN, false, #self#, poly) │┌ @ /home/nsajko/.julia/packages/PolynomialRoots/1iZQh/src/PolynomialRoots.jl:617 1 PolynomialRoots.:(:) %28 ││ runtime dispatch detected: (1 PolynomialRoots.:(:) %28::Union{Nothing, Int64})::UnitRange{Int64} │└───────────────────────────────────────────────────────────────────────────────── ``` After this commit: ``` julia> using JET, PolynomialRoots julia> PolynomialRoots.roots(rand(5)) 4-element Vector{ComplexF64}: -31.46526111710419 - 1.6981590262940853e-16im -0.13020601478750077 - 0.8162149121849699im -0.13020601478750077 + 0.81621491218497im -0.48806837350069937 + 2.350988701644575e-38im julia> @report_opt PolynomialRoots.roots(rand(5)) No errors detected ``` --- src/PolynomialRoots.jl | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/src/PolynomialRoots.jl b/src/PolynomialRoots.jl index 1722992..a1e4d87 100644 --- a/src/PolynomialRoots.jl +++ b/src/PolynomialRoots.jl @@ -607,15 +607,20 @@ function roots(poly::AbstractVector{<:Number}, roots::AbstractVector{<:Number}; epsilon, degree, polish) end +# Returns a view into the polynomial that excludes any trailing zero +# coefficients. +# +# Errors if the coefficients are all zero. +function normalized_poly_coefs(poly::AbstractVector{<:Number}) + l = findlast(!iszero, poly) + isnothing(l) && error("zero polynomial") + @view poly[begin:l] +end + function roots(poly::AbstractVector{N}; epsilon::AbstractFloat=NaN, polish::Bool=false) where {N<:Number} # Before starting, truncate the polynomial if it has zeros in the trailing elements - last_nz = findlast(!iszero, poly) - if lastindex(poly) == last_nz - _poly = poly - else - _poly = poly[1:last_nz] - end + _poly = normalized_poly_coefs(poly) degree = length(_poly) - 1 roots!(zeros(Complex{real(float(N))}, degree), float.(complex(_poly)), epsilon, degree, polish)