On 30.09.2015 13:21, Pauli Virtanen wrote:
Juha Jeronen <juha.jeronen <at> jyu.fi> writes:
I recently developed a Cython-based, OpenMP-accelerated quartic (and cubic, quadratic) polynomial solver to address a personal research need for quickly solving a very large number of independent low-degree polynomial equations for both real and complex coefficients. My 2c in this context would be to also think how this best fits in how collections of polynomials are represented in Numpy.
AFAICS, Numpy's polynomial module supports evaluation of polynomial collections (cf. numpy.polynomial.polynomial.polyval), but the corresponding root finding routine (numpy.polynomial.polynomial.polyroots) only deals with one polynomial at a time.
The present code in principle could be used to speed up the latter routine after it is generalized to multiple polynomials. The general case is probably doable just by coding up the companion matrix approach using low-level routines (or maybe with the new vectorized np.linalg routines).
Some relevant code elsewhere for similar purposes can be found in scipy.inteprolate.PPoly/BPoly (and in the future BSpline).
https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.PPoly...
However, since it's piecewise, there's purposefully support only for real-valued roots.
Thanks. It sounds like numpy.polynomial.polynomial.polyroots() is a good candidate place for this, if as you suggest, it is first generalized to handle multiple polynomials. -J