Special support for roots when the coefficients are real?

[nsajko]

Empirically, it seems that roots(::Vector{Float64}) returns the roots such that they exactly conform to the form that follows from the complex conjugate root theorem. That is, a real root has its imaginary part as exactly zero, and complex roots come in complex conjugate pairs, where the complex conjugate operation is exact. Example:

julia> import AMRVW

julia> const A = AMRVW
AMRVW

julia> A.roots(rand(20))
19-element Vector{ComplexF64}:
  -1.2606264533481746 + 0.0im
  -0.9310479653026866 - 0.48344571715462364im
  -0.9310479653026866 + 0.48344571715462364im
  -0.7452679583988157 + 0.0im
 -0.48547663249542244 - 0.8197940002283542im
 -0.48547663249542244 + 0.8197940002283542im
 -0.30132852172591407 - 1.010361296238549im
 -0.30132852172591407 + 1.010361296238549im
  -0.2338481178786882 - 1.3556065317036028im
  -0.2338481178786882 + 1.3556065317036028im
  -0.1383958798698237 + 0.0im
  0.34088495469868396 - 1.0901157173978522im
  0.34088495469868396 + 1.0901157173978522im
  0.34435589931079486 - 0.8281753378613422im
  0.34435589931079486 + 0.8281753378613422im
    0.754774828425891 - 0.5684903527552175im
    0.754774828425891 + 0.5684903527552175im
   0.9427060520193593 - 0.3181694681365702im
   0.9427060520193593 + 0.3181694681365702im

I wonder if (or to what degree) is this guaranteed or (alternatively) an accident?

If it is guaranteed, I think it would make sense if this package exposed a new interface that would return the roots in a less redundant and more convenient form, as expected from the complex conjugate root theorem. So, instead of returning a Vector{Complex{<:AbstractFloat}}, this new interface would return an instance of a type like this, where the complex conjugate root pairs are stored as only the root whose imaginary part is positive:

struct ComplexConjugateRootTheoremRoots{F<:AbstractFloat}
  r::Vector{F}

  # The imaginary part of each element is guaranteed to be positive.
  c::Vector{Complex{F}}
end

How much sense am I making?

[jverzani]

It is guaranteed, the algorithm reduces the matrix until there is a linear or quadratic to evaluate. I'm not sure what is more convenient here though. Without thinking too much, this seems a bit complicated, but you could play around with it.

The key lines in the https://github.com/jverzani/AMRVW.jl/blob/master/src/AMRVW_algorithm.jl file are 50, 79, 102, and 112. Unless I missed something, here is where EIGS is written to.

[nsajko]

Without thinking too much, this seems a bit complicated, but you could play around with it.

The key lines in the https://github.com/jverzani/AMRVW.jl/blob/master/src/AMRVW_algorithm.jl file are 50, 79, 102, and 112. Unless I missed something, here is where EIGS is written to.

The linked PR implements the feature, and it's actually trivial. I suppose your comment discusses a hypothetical lower-level implementation that would be more efficient? Still, I suppose my proof-of-concept PR is good for a start.