Skip to content

Commit

Permalink
prevent run-time dispatch in roots(::AbstractVector)
Browse files Browse the repository at this point in the history
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
```
  • Loading branch information
nsajko committed Oct 12, 2022
1 parent ccc2222 commit de52131
Showing 1 changed file with 11 additions and 6 deletions.
17 changes: 11 additions & 6 deletions src/PolynomialRoots.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down

0 comments on commit de52131

Please sign in to comment.