PolynomialSolver infinite loop for overflowed coefficients

Summary

Companion matrix balancing get stuck at infinite loop if coefficients are greater than 1e+308/infinity (x64).

Environment

  • Operating System : Linux
  • Architecture : x64
  • Eigen Version : 3.3.7
  • Compiler Version : CLANG 8.0

Minimal Example

    double const a = 1e+308; // infinite loop
//    auto const a = std::numeric_limits<double>::infinity();  // infinite loop
//    double const a = 1e+307;  // doesn't get stuck
    double const b = -1;
    double const c = 1;

    auto const coeffs = Eigen::Vector4d(a, b, c, 1).eval();  // infinite loop

    auto psolver = Eigen::PolynomialSolver<double, 3>(coeffs);  // performs balancing at construction

    psolver.roots(); // unreachable code

What is the current bug behavior?

At Companion.h (https://gitlab.com/libeigen/eigen/-/blob/master/unsupported/Eigen/src/Polynomials/Companion.h#L164), there is a second loop that has the following code:

    while (colNorm >= rowB)
    {
      colB /= radix<Scalar>();
      colNorm /= radix2<Scalar>();
    }

It is expected to have infinite values for colNorm and rowB due the given coefficients. However, in this case, the evaluation of inf >= inf gives true, in opposite of what happens with inf > inf, which gives false.

In my particular case, I was avoiding this condition by checking each coefficient with std::isfinite. However, as shown in my code, when a=1e+308, it scapes from the checking and the same problem happens.

What is the expected correct behavior?

Independently of the result (probably garbage), the balancing should be able to avoid infinite loop. A quick fix is to use colNorm > rowB (instead of >=), so the inf comparison problem goes away. As this is about floating points, I think that the equality doesn't play a major role there anyway, but I'm not very confident about it.

  • Have a plan to fix this issue.
Edited by Danilo Silva