search

giordano / PolynomialRoots.jl

Fast complex polynomial root finder, with support for arbitrary precision calculations
Other
54 stars 15 forks source link

Support for StaticArrays (MVector) #9

Open fgasdia opened 6 years ago

fgasdia commented 6 years ago

This package is strictly coded for Julia native Vectors. I can't tell just by looking at the source if there would be a performance speedup by generalizing the types to allow for S or MVectors of the StaticArrays package, but it's the preferred Julia style to allow as general type arguments as possible.

giordano commented 6 years ago

I don't think this will be possible. The roots function internally calls a roots! function which operates in-place, so a StaticVector wouldn't work. In terms of performance, a function acting in-place is already good.

fgasdia commented 6 years ago

The SVector is immutable, but the MVector type is only static in the sense of its size (it's mutable) and are extremely fast with mutating methods. Though I admit that generalizing the types to e.g. AbstractVector may result in less graceful errors if incompatible types are used...

giordano commented 6 years ago

I've opened #10 to replace all Vectors with AbstractVectors. It's good in general to have a more general signature, but I'm afraid this is not going to help to use SVector for the above-mentioned reasons

fgasdia commented 6 years ago

Thanks. MVector is currently being reworked so that they stay mutable between operations (they can currently return SVectors), see https://github.com/JuliaArrays/StaticArrays.jl/pull/536. So if the zeros() in roots() returns not just the Complex float promoted zero Vector for the initial roots, but a vector of the same AbstractVector type as the input poly is, then it looks like that may now cascade automatically through the rest of the code since it is now all AbstractVector.

In other words, if poly is of type Vector, then roots would also be of type Vector, but if poly is of type MVector (or whatever AbstractVector), then roots would also be of type MVector.

It's not obvious how to get zeros() to return the promoted eltype with the correct container type though.

giordano commented 5 years ago 
fill!(similar(float.(complex(poly))), 0)

should do the trick

giordano commented 5 years ago

That would actually need to be

fill!(similar(float.(complex(poly)), (degree,)), 0)

but doesn't work because

similar(float.(complex(poly)), (degree,))

gives a Vector{Complex}, instead of MVector{Complex}, contrary to 

similar(float.(complex(poly)))

which gives an MVector as expected. I don't know if this is a bug or an intended behavior. Do you mind opening a ticket there?

fgasdia commented 5 years ago

This was actually just brought up by someone else last week https://github.com/JuliaArrays/StaticArrays.jl/issues/545. The current suggestion is to use similar(float.(complex(poly)), Size(degree,)), where Size() is a function from StaticArrays, but that's not a good solution here because it produces a SizedArray and requires importing StaticArrays.