Adding Cauchy based bound estimate to polysolve.

Added code for an original Cauchy bound estimate good for the lower root
and a Cauchy based upper bound method good for all roots. The latter is
currently hard coded on and this value is used over MAX_DISTANCE.

Longer term want to allow passed value and option to pick either of the two
methods added with this commit or use MAX_DISTANCE.

Author: wfpokorny

source/core/math/polynomialsolver.cpp

@@ -1513,11 +1513,10 @@ static int solve_quartic(const DBL *x, DBL *results)
 static int polysolve(int order, const DBL *Coeffs, DBL *roots)
 {
     polynomial sseq[MAX_ORDER+1];
-    DBL min_value, max_value;
+    DBL min_value, max_value, max_value2, Abs_Coeff_n;
     int i, nroots, np, atmin, atmax;
 
-    /* Put the coefficients into the top of the stack. */
-
+    // Reverse coefficients into order used herein.
     for (i = 0; i <= order; i++)
     {
         sseq[0].coef[order-i] = Coeffs[i];
@@ -1534,11 +1533,46 @@ static int polysolve(int order, const DBL *Coeffs, DBL *roots)
         return(0);
     }
 
-    /* Bracket the roots */
+    // Bracket the roots
 
     min_value = 0.0;
     max_value = MAX_DISTANCE;
 
+    // Optionally use Augustin-Louis Cauchy's methods to determine an upper bound for max_value.
+
+    // This is the lower bound or bound including, for certain, only the first root but
+    // perhaps others.
+    if (0)
+    {
+        max_value   = fabs(Coeffs[order]);
+        Abs_Coeff_n = fabs(Coeffs[0]);
+        for (i = 1; i < order; i++)
+        {
+            max_value = max(fabs(Coeffs[i]),max_value);
+        }
+        max_value /= Abs_Coeff_n + 1;
+        max_value = min(max_value,MAX_DISTANCE);
+    }
+
+    // Tighter upper bound found at:
+    //    https://en.wikipedia.org/wiki/Properties_of_polynomial_roots#Other_bounds
+    // which took it from:
+    //    Cohen, Alan M. (2009). "Bounds for the roots of polynomial equations".
+    //    Mathematical Gazette. 93: 87-88.
+    // NOTE: Had to use > 1.0 in max_value2 calculation in practice...
+    if (1)
+    {
+        Abs_Coeff_n = fabs(Coeffs[0]); // Solve_Polynomial() dumps leading zeroes.
+        max_value2  = 1.1 + fabs(Coeffs[1]/Abs_Coeff_n);
+        max_value   = fabs(Coeffs[2]);
+        for (i = 3; i <= order; i++)
+        {
+            max_value = max(fabs(Coeffs[i]),max_value);
+        }
+        max_value /= Abs_Coeff_n + EPSILON;
+        max_value = min(max(max_value,max_value2),MAX_DISTANCE);
+    }
+
     atmin = numchanges(np, sseq, min_value);
     atmax = numchanges(np, sseq, max_value);