Non-deterministic behavior

[waruyama]

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)

[Feirell]

Hi @waruyama even though I am neither an author of this library nor am I sure that my suggestion is fail poof, so take this with a grain of salt, I think I have an idea. You can make this determenistic by supplying the initial guesses. As you wrote the library uses Math.random for the guesses which results in your problem.

You could pick the initials by using the same conditions as was done in the wikipedia article:

(0.4 + 0.9i)^0 ==  1          + 0i
(0.4 + 0.9i)^1 ==  0.4        + 0.9i
(0.4 + 0.9i)^2 == -0.65       + 0.72i
(0.4 + 0.9i)^3 == -0.908      - 0.297i
(0.4 + 0.9i)^4 == -0.0959     - 0.936i
(0.4 + 0.9i)^5 ==  0.80404    - 0.46071i
(0.4 + 0.9i)^6 ==  0.736255   + 0.539353i
(0.4 + 0.9i)^7 == -0.1909148  + 0.8783703i
(0.4 + 0.9i)^8 == -0.86689919 + 0.1795248i

So you would call findRoots like this:

const guessR = [
     1,
     0.4,
    -0.65,
    -0.908,
    -0.0959,
     0.80404,
     0.736255,
    -0.1909148,
    -0.86689919
];

const guessI = [
     0,
     0.9,
     0.72,
    -0.297,
    -0.936,
    -0.46071,
     0.539353,
     0.8783703,
     0.1795248
];

const coeffR = [1, 9, 36, 69, 36, -2043, 3747, -7092, 5328, -3520];

findRootsNpmPack(coeffR, undefined, undefined, undefined, guessR, guessI);

Which would result in the following, rounded, results:

 0.4095 + 0.0000i
 0.5638 + 0.8316i
-0.0407 + 0.2083i
-0.0407 - 0.2083i
 0.1703 - 0.7626i
 0.5638 - 0.8316i
-0.1413 - 0.0709i
 0.1703 + 0.7626i
-0.1413 + 0.0709i

As I said I am not sure if this would hold up for all possible coefficients but it at least gives you a function which either fails or works correctly. Another idea would be to create two or three sets of initial guesses which could be used when the previous failed.

For your third point you could take a look here https://en.wikipedia.org/wiki/Root-finding_algorithms#Finding_all_roots_at_once because I don't think that there is a closed form for a polynomial of grade 9 so you would always rely on an iterative approximation.