Skip to content

Commit

Permalink
Merge pull request #16 from andreasnoack/an/polynomial
Browse files Browse the repository at this point in the history
Make Polynomials a dependency and stop using Requires
  • Loading branch information
andreasnoack authored Apr 18, 2020
2 parents 82a4f42 + 49a25e6 commit 99d8a2f
Show file tree
Hide file tree
Showing 4 changed files with 16 additions and 25 deletions.
9 changes: 4 additions & 5 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -4,16 +4,15 @@ version = "1.0.0"

[deps]
LibAMVW_jll = "3420f54a-cca3-5752-a4b4-c8f5453f04ac"
Requires = "ae029012-a4dd-5104-9daa-d747884805df"
Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"

[compat]
julia = "1.3"
LibAMVW_jll = "1"
Requires = "1"
Polynomials = "0.7"
julia = "1.3"

[extras]
Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[targets]
test = ["Polynomials", "Test"]
test = ["Test"]
13 changes: 5 additions & 8 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -7,21 +7,19 @@ This package is a Julia wrapper of the Fortran programs accompanying [Fast and B
## Usage

The package provides the unexported function `FastPolynomialRoots.rootsFastPolynomialRoots(p::Vector{<:Union{Float64,Complex{Float64}}})`
which computes the roots of the polynomial `p[1] + p[2]*x + p[3]*x^2 + ... + p[k]*x^(k-1)`. If the
`Polynomials` packages is loaded, the `roots(::Poly)` methods for `Float64` and `Complex{Float64}` will
be overwritten to that fast version provided by this package. See the examples below.
which computes the roots of the polynomial `p[1] + p[2]*x + p[3]*x^2 + ... + p[k]*x^(k-1)`. The package also overwrites the `roots(::Polynomial)` methods in the `Polynomials` package for `Float64` and `Complex{Float64}` elements with the fast versions provided by this package. See the examples below.

## Example 1: Speed up `roots`
```julia
julia> using Polynomials, BenchmarkTools

julia> @btime roots(p) setup=(p = Polynomial(randn(500)));
408.564 ms (54 allocations: 3.99 MiB)
223.135 ms (23 allocations: 3.97 MiB)

julia> using FastPolynomialRoots

julia> @btime roots(p) setup=(p = Polynomial(randn(500)));
46.507 ms (7 allocations: 26.41 KiB)
30.786 ms (7 allocations: 26.41 KiB)
```

## Example 2: Roots of a polynomial of degree 10,000
Expand All @@ -30,11 +28,10 @@ but can be handled by FastPolynomialRoots.
```julia
julia> using Polynomials, BenchmarkTools, FastPolynomialRoots

julia> n = 10000
10000
julia> n = 10000;

julia> r = @btime roots(p) setup=(p = Polynomial(randn(n + 1)));
15.715 s (13 allocations: 508.38 KiB)
10.290 s (13 allocations: 508.38 KiB)

julia> sum(isreal, r)
7
Expand Down
11 changes: 3 additions & 8 deletions src/FastPolynomialRoots.jl
Original file line number Diff line number Diff line change
@@ -1,14 +1,9 @@
module FastPolynomialRoots

using LibAMVW_jll, Requires
using LibAMVW_jll, Polynomials

function __init__()

@require Polynomials="f27b6e38-b328-58d1-80ce-0feddd5e7a45" begin
Polynomials.roots(p::Union{Polynomials.Poly{Float64},Polynomials.Poly{Complex{Float64}}}) = rootsFastPolynomialRoots(p.a)
Polynomials.roots(p::Polynomials.Poly{<:Integer}) = rootsFastPolynomialRoots(convert(Polynomials.Poly{Float64}, p))
end
end
Polynomials.roots(p::Union{Polynomial{Float64},Polynomial{Complex{Float64}}}) = rootsFastPolynomialRoots(coeffs(p))
Polynomials.roots(p::Polynomial{<:Integer}) = rootsFastPolynomialRoots(convert(Polynomial{Float64}, p))

function rootsFastPolynomialRoots(a::Vector{Float64})

Expand Down
8 changes: 4 additions & 4 deletions test/runtests.jl
Original file line number Diff line number Diff line change
@@ -1,21 +1,21 @@
using Test, FastPolynomialRoots, Polynomials

@testset "Standard normal coefficients" begin
p = Poly(randn(50))
p = Polynomial(randn(50))
@test sort(abs.(roots(p))) sort(abs.(roots(p)))
end
@testset "Standard normal complex coefficient" begin
p = Poly(complex.(randn(50), randn(50)))
p = Polynomial(complex.(randn(50), randn(50)))
@test sort(abs.(roots(p))) sort(abs.(roots(p)))
end

@testset "Large polynomial" begin
p = Poly(randn(5000))
p = Polynomial(randn(5000))
@time roots(p)

@info "Possible to calculate roots of large polynomial"
@show rts = 1:100.0
p = mapreduce(t -> Poly([1, -t]), *, rts, init=Poly(Float64[1]))
p = mapreduce(t -> Polynomial([1, -t]), *, rts, init=Polynomial(Float64[1]))
@info "But polynomial root finding is ill conditioned"
@show sum(abs2, sort(abs.(roots(p))) - rts)
end

2 comments on commit 99d8a2f

@andreasnoack
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request updated: JuliaRegistries/General/13212

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v1.0.0 -m "<description of version>" 99d8a2f2fe4681bfb18b4e93b551fdb84acd97b9
git push origin v1.0.0

Please sign in to comment.