Hi Steven,

2016-04-28 14:58 GMT+02:00 Steven G. Johnson <[email protected]>:
> On Thursday, April 28, 2016 at 8:46:06 AM UTC-4, Steven G. Johnson wrote:
>>
>> Out of curiosity, how does the performance and accuracy compare to the
>> "roots" function in the Polynomials.jl package, which uses the standard
>> companion-matrix approach?
>
>
> Just checked: it is much faster than Polynomials.roots, which I guess is not
> surprising for a Newton-like method compared to an eigenvalue method,

Yes, as I wrote in my first mail I found my package to be ~20 times
faster than Polynomials.jl and usually more precise (**very** roughly,
evaluating abs2(evaluation of polynomial at root) result in ~10^-29
for Polynomials.jl and ~10^-30 for CmplxRoots.jl).  Probably I
wouldn't have published the package if it had been both slower and less precise.

> although perhaps the latter is more reliable in some cases?
>
> I'm not sure what the state of the art here is, though.

Most of the methods I found are based on finding eigenvalues of
companion matrix.  A somewhat popular method is Jenkins–Traub
algorithm (https://en.wikipedia.org/wiki/Jenkins%E2%80%93Traub_algorithm),
but I didn't find a FLOSS implementation of that algorithm in any
language, to start with.

In my test suite I added a couple of tests of random polynomials, to
catch major flaws of root finding algorithm.

> (I find it a little bit worrying that the Skowron and Gould paper justify
> their method on account of it being "1.6 to 3 times" faster than the roots
> function in Numerical Recipes (NR), and that NR was the inspiration for
> their algorithm.  NR routines typically have deeply suboptimal performance,
> and that book is terrible starting point for developing state-of-the-art
> algorithms because NR's approaches are typically antiquated.)

I know that algorithm showed in Numerical Recipes are often far from
being the state of art (in addition, their routines aren't FLOSS), but
that series of book is fairly common in astrophysics even today, maybe
because at least the first two authors are well known astrophysicists.

I started from Skowron & Gould package because their package is FLOSS,
it's in Fortran and it's easy to be called from Julia (indeed I first
wrote a wrapper around the Fortran shared object, that took me a few
minutes), and I used it often, without major problems.

Bye,
Mosè

Reply via email to