I ran a performance test on this package and incidentally found that the results of the root finding are not deterministic (the performance results are great, btw).

If I run the code below 200000 times there are usually one or two cases when the algorithm fails to find the real root, Sometimes none of the found roots is not even close to the correct result, which is at 0.409458563186124. Changing tolerance values did not help so far.

for (let i = 0; i< 200000; i++) {
     window.r = [];
     let r = findRoots([1, 9, 36, 69, 36, -2043, 3747, -7092, 5328, -3520]);
     let realR = [];
     for (let j = 0; j < r[0].length; j++) {
         if (Math.abs(r[1][j]) < 1e-10) {
             realR.push(r[0][j]);
         }
     }
     if (realR.length == 0) {
         console.log("failed", r, window.r);

     }
 }
 console.log("done")

This non-deterministic if caused by the use of Math.random(). I recorded the generated random values to see if they are special in the failing cases, but I could not see a pattern at all.

My question are

1. Can I do something to make the algorithm fail less often?
2. Is it possible to make the algorithm deterministic?
3. Is there a better algorithm for my case (only real roots needed, always polynomials of degree 9)

I ran a couple of tests against another implementation of the Weierstrass / Durand-Kerner method and noticed non-convergence in the JS code. The JS code computes the reciprocal of (a,b) for the denominator:

      //Compute reciprocal
      k1 = a*a + b*b
      if(abs(k1) > EPSILON) {
        a /=  k1
        b /= -k1
      } else {
        a = 1.0
        b = 0.0
      }

This operation when combined with the multiplication by the numerator (na,nb) can lead to loss of precision! To avoid this:

        // compute (qa,qb)=(na,nb)/(a,b) without loss of precision, overflow/underflow
        if (abs(a) >= abs(b)) {
          k1 = b/a
          k2 = a + k1*b
          qa = (na + k1*nb)/k2
          qb = (nb - k1*na)/k2
        } else {
          k1 = a/b
          k2 = k1*a + b
          qa = (k1*na + nb)/k2
          qb = (k1*nb - na)/k2
        }

To verify and compare the correction, use a polynomial of high degree with small coefficients. For example, a 100 degree polynomial with all coefficients set to 1 (that is, zr[] has 101 ones and zi[] has 101 zeros).

To complete this change, it may still be good to check that a*a + b*b <= EPSILON by adding this guard:

      k1 = a*a + b*b
      if(abs(k1) > EPSILON) {
        // compute (qa,qb)=(na,nb)/(a,b) without loss of precision, overflow/underflow
        if (abs(a) >= abs(b)) {
          k1 = b/a
          k2 = a + k1*b
          qa = (na + k1*nb)/k2
          qb = (nb - k1*na)/k2
        } else {
          k1 = a/b
          k2 = k1*a + b
          qa = (k1*na + nb)/k2
          qb = (k1*nb - na)/k2
        }
      } else {
        qa = na
        qb = nb
      }

Hello,

I've been using this package to compute roots with varying succes
(x+1)^3 for instance did not always return {-1,-1,-1} but had {-0.9995-0.005403i, -1.004+0.001627i,-0.9948+0.003137i} for 100nn iterations and tolerance of 1e-9.

By increasing the EPSILON from 1e-8 tot 2e-16 (the machine epsilon of Javascript) inside the package itself, I've managed to get {1,1,1}.

Can you perhaps explain why EPSILON has this effect and if I should do this a different way?