Rewrite balancer to avoid overflows.

The previous balancer overflowed for large row/column norms.
Modified to prevent that.

Fixes #2273.

1 files changed, 26 insertions, 21 deletions

diff --git a/unsupported/Eigen/src/Polynomials/Companion.h b/unsupported/Eigen/src/Polynomials/Companion.h
index 6ab8f9714..59a15b098 100644
--- a/unsupported/Eigen/src/Polynomials/Companion.h
+++ b/unsupported/Eigen/src/Polynomials/Companion.h
@@ -20,12 +20,6 @@ namespace internal {
 
 #ifndef EIGEN_PARSED_BY_DOXYGEN
 
-template <typename T>
-T radix(){ return 2; }
-
-template <typename T>
-T radix2(){ return radix<T>()*radix<T>(); }
-
 template<int Size>
 struct decrement_if_fixed_size
 {
@@ -141,7 +135,10 @@ inline
 bool companion<_Scalar,_Deg>::balanced( RealScalar colNorm, RealScalar rowNorm,
     bool& isBalanced, RealScalar& colB, RealScalar& rowB )
 {
-  if( RealScalar(0) == colNorm || RealScalar(0) == rowNorm ){ return true; }
+  if( RealScalar(0) == colNorm || RealScalar(0) == rowNorm 
+      || !(numext::isfinite)(colNorm) || !(numext::isfinite)(rowNorm)){
+    return true;
+  }
   else
   {
     //To find the balancing coefficients, if the radix is 2,
@@ -149,33 +146,41 @@ bool companion<_Scalar,_Deg>::balanced( RealScalar colNorm, RealScalar rowNorm,
     // \f$ 2^{2\sigma-1} < rowNorm / colNorm \le 2^{2\sigma+1} \f$
     // then the balancing coefficient for the row is \f$ 1/2^{\sigma} \f$
     // and the balancing coefficient for the column is \f$ 2^{\sigma} \f$
-    rowB = rowNorm / radix<RealScalar>();
+    const RealScalar radix = RealScalar(2);
+    const RealScalar radix2 = RealScalar(4);
+    
+    rowB = rowNorm / radix;
     colB = RealScalar(1);
     const RealScalar s = colNorm + rowNorm;
 
-    while (colNorm < rowB)
+    // Find sigma s.t. rowNorm / 2 <= 2^(2*sigma) * colNorm
+    RealScalar scout = colNorm;
+    while (scout < rowB)
     {
-      colB *= radix<RealScalar>();
-      colNorm *= radix2<RealScalar>();
+      colB *= radix;
+      scout *= radix2;
     }
-
-    rowB = rowNorm * radix<RealScalar>();
-
-    while (colNorm >= rowB)
+    
+    // We now have an upper-bound for sigma, try to lower it.
+    // Find sigma s.t. 2^(2*sigma) * colNorm / 2 < rowNorm
+    scout = colNorm * (colB / radix) * colB;  // Avoid overflow.
+    while (scout >= rowNorm)
     {
-      colB /= radix<RealScalar>();
-      colNorm /= radix2<RealScalar>();
+      colB /= radix;
+      scout /= radix2;
     }
 
-    //This line is used to avoid insubstantial balancing
-    if ((rowNorm + colNorm) < RealScalar(0.95) * s * colB)
+    // This line is used to avoid insubstantial balancing.
+    if ((rowNorm + radix * scout) < RealScalar(0.95) * s * colB)
     {
       isBalanced = false;
       rowB = RealScalar(1) / colB;
       return false;
     }
-    else{
-      return true; }
+    else
+    {
+      return true;
+    }
   }
 }