RpolyPlusPlus icon indicating copy to clipboard operation
RpolyPlusPlus copied to clipboard

Incorrect roots found

Open advanpix opened this issue 5 years ago • 3 comments

I have been testing the library on some polynomials with integer coefficients. On many of them incorrect roots are found. One of the examples:

Coefficients: [1,0,-1,0,1,1,1,1,1,1,1,0,-1,0,1]

Roots found:

  1.1910 + 0.6054i
  1.1910 - 0.6054i
 -1.2766 + 0.3311i
 -1.2766 - 0.3311i
 -0.9114 + 0.8715i
 -0.9114 - 0.8715i
 -0.2781 + 1.1229i
 -0.2781 - 1.1229i
  0.4173 + 0.9860i
  0.4173 - 0.9860i
  0.6739 + 0.3241i
  0.6739 - 0.3241i
  0.1838 + 0.2115i
  0.1838 - 0.2115i

True roots:

  1.2051 + 0.5795i
  1.2051 - 0.5795i
 -0.5464 + 0.8375i
 -0.5464 - 0.8375i
 -0.9793 + 0.2023i
 -0.9793 - 0.2023i
 -0.8508 + 0.5254i
 -0.8508 - 0.5254i
  0.5333 + 0.8459i
  0.5333 - 0.8459i
 -0.0358 + 0.9994i
 -0.0358 - 0.9994i
  0.6739 + 0.3241i
  0.6739 - 0.3241i

Roots are sorted by magnitude.

advanpix avatar Apr 11 '20 05:04 advanpix

I cannot reproduce the issue. This is the result I get:

     (0.673943,0.3241)
    (0.673943,-0.3241)
   (-0.85084,0.525425)
  (-0.85084,-0.525425)
  (-0.979328,0.202277)
 (-0.979328,-0.202277)
   (0.533286,0.845935)
  (0.533286,-0.845935)
 (-0.0357928,0.999359)
(-0.0357928,-0.999359)
  (-0.546374,0.837541)
 (-0.546374,-0.837541)
    (1.20511,0.579537)
   (1.20511,-0.579537)

sergiud avatar Apr 11 '20 13:04 sergiud

I did few other tests, accuracy varies greatly depending on parameters (tolerances, number of shifts, etc.). What parameters do you use? What are residuals for each root? I use defaults (Intel C++ on Windows).

I work with reciprocal polynomials with integer coefficients, most of them have ill-conditioned roots (especially cyclotomic polys).

So far I was not able to get consistent accuracy with RpolyPlusPlus (one set of parameters work for one polynomial but not for the other). I compare it with DHSEQR from MKL/LAPACK, which is very stable (but slow).

One more thing, it is more numerically robust to use std::hypot to compute magnitude of a complex number. It takes care of the overflow/underflow and other accuracy losses.

advanpix avatar Apr 12 '20 02:04 advanpix

I used the default parameters. Did you enable any particular optimization flags when compiling using Intel C++?

sergiud avatar Apr 12 '20 08:04 sergiud