Find Roots: extend bracketing in derivative-free algorithms to also reduce/subdevision after trying to expand

src/Numerics/FindRoots.cs

@@ -49,7 +49,7 @@ namespace MathNet.Numerics
         {
             double root;
 
-            if (!ZeroCrossingBracketing.Expand(f, ref lowerBound, ref upperBound, 1.6, 100))
+            if (!ZeroCrossingBracketing.ExpandReduce(f, ref lowerBound, ref upperBound, 1.6, maxIterations, maxIterations*10))
             {
                 throw new NonConvergenceException(Resources.RootFindingFailed);
             }

src/Numerics/RootFinding/Bisection.cs

@@ -49,7 +49,7 @@ namespace MathNet.Numerics.RootFinding
         /// <exception cref="NonConvergenceException"></exception>
         public static double FindRootExpand(Func<double, double> f, double guessLowerBound, double guessUpperBound, double accuracy = 1e-8, int maxIterations = 100, double expandFactor = 1.6, int maxExpandIteratons = 100)
         {
-            ZeroCrossingBracketing.Expand(f, ref guessLowerBound, ref guessUpperBound, expandFactor, maxExpandIteratons);
+            ZeroCrossingBracketing.ExpandReduce(f, ref guessLowerBound, ref guessUpperBound, expandFactor, maxExpandIteratons, maxExpandIteratons*10);
             return FindRoot(f, guessLowerBound, guessUpperBound, accuracy, maxIterations);
         }
 

src/Numerics/RootFinding/Brent.cs

@@ -50,7 +50,7 @@ namespace MathNet.Numerics.RootFinding
         /// <exception cref="NonConvergenceException"></exception>
         public static double FindRootExpand(Func<double, double> f, double guessLowerBound, double guessUpperBound, double accuracy = 1e-8, int maxIterations = 100, double expandFactor = 1.6, int maxExpandIteratons = 100)
         {
-            ZeroCrossingBracketing.Expand(f, ref guessLowerBound, ref guessUpperBound, expandFactor, maxExpandIteratons);
+            ZeroCrossingBracketing.ExpandReduce(f, ref guessLowerBound, ref guessUpperBound, expandFactor, maxExpandIteratons, maxExpandIteratons*10);
             return FindRoot(f, guessLowerBound, guessUpperBound, accuracy, maxIterations);
         }
 

src/Numerics/RootFinding/RobustNewtonRaphson.cs

@@ -2,9 +2,9 @@
 // Math.NET Numerics, part of the Math.NET Project
 // http://numerics.mathdotnet.com
 // http://github.com/mathnet/mathnet-numerics
-// 
+//
 // Copyright (c) 2009-2013 Math.NET
-// 
+//
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
 // files (the "Software"), to deal in the Software without
@@ -13,10 +13,10 @@
 // copies of the Software, and to permit persons to whom the
 // Software is furnished to do so, subject to the following
 // conditions:
-// 
+//
 // The above copyright notice and this permission notice shall be
 // included in all copies or substantial portions of the Software.
-// 
+//
 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
 // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND

src/Numerics/RootFinding/ZeroCrossingBracketing.cs

@@ -35,7 +35,7 @@ namespace MathNet.Numerics.RootFinding
 {
     public static class ZeroCrossingBracketing
     {
-        public static IEnumerable<Tuple<double, double>> FindIntervalsWithin(Func<double, double> f, double lowerBound, double upperBound, int parts)
+        public static IEnumerable<Tuple<double, double>> FindIntervalsWithin(Func<double, double> f, double lowerBound, double upperBound, int subdivisions)
         {
             // TODO: Consider binary-style search instead of linear scan
             double fmin = f(lowerBound);
@@ -47,11 +47,11 @@ namespace MathNet.Numerics.RootFinding
                 yield break;
             }
 
-            double subdiv = (upperBound - lowerBound)/parts;
+            double subdiv = (upperBound - lowerBound)/subdivisions;
             double smin = lowerBound;
             int sign = Math.Sign(fmin);
 
-            for (int k = 0; k < parts; k++)
+            for (int k = 0; k < subdivisions; k++)
             {
                 double smax = smin + subdiv;
                 double sfmax = f(smax);
@@ -116,5 +116,59 @@ namespace MathNet.Numerics.RootFinding
             upperBound = originalUpperBound;
             return false;
         }
+
+        public static bool Reduce(Func<double, double> f, ref double lowerBound, ref double upperBound, int subdivisions = 1000)
+        {
+            double originalLowerBound = lowerBound;
+            double originalUpperBound = upperBound;
+
+            if (lowerBound >= upperBound)
+            {
+                throw new ArgumentOutOfRangeException("upperBound", string.Format(Resources.ArgumentOutOfRangeGreater, "xmax", "xmin"));
+            }
+
+            // TODO: Consider binary-style search instead of linear scan
+            double fmin = f(lowerBound);
+            double fmax = f(upperBound);
+
+            if (Math.Sign(fmin) != Math.Sign(fmax))
+            {
+                return true;
+            }
+
+            double subdiv = (upperBound - lowerBound) / subdivisions;
+            double smin = lowerBound;
+            int sign = Math.Sign(fmin);
+
+            for (int k = 0; k < subdivisions; k++)
+            {
+                double smax = smin + subdiv;
+                double sfmax = f(smax);
+                if (double.IsInfinity(sfmax))
+                {
+                    // expand interval to include pole
+                    smin = smax;
+                    continue;
+                }
+
+                if (Math.Sign(sfmax) != sign)
+                {
+                    lowerBound = smin;
+                    upperBound = smax;
+                    return true;
+                }
+
+                smin = smax;
+            }
+
+            lowerBound = originalLowerBound;
+            upperBound = originalUpperBound;
+            return false;
+        }
+
+        public static bool ExpandReduce(Func<double, double> f, ref double lowerBound, ref double upperBound, double expansionFactor = 1.6, int expansionMaxIterations = 50, int reduceSubdivisions = 100)
+        {
+            return Expand(f, ref lowerBound, ref upperBound, expansionFactor, expansionMaxIterations) || Reduce(f, ref lowerBound, ref upperBound, reduceSubdivisions);
+        }
     }
 }

src/UnitTests/RootFindingTests/FindRootsTest.cs

@@ -160,6 +160,16 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests
             Assert.AreEqual(0, f1(x), 1e-5);
         }
 
+        [Test]
+        public void StackOverflow39935588()
+        {
+            // Roots at -2, 2
+            Func<double, double> f1 = x => (x - 3.0)*(x - 4.0);
+            Assert.AreEqual(3.0, FindRoots.OfFunction(f1, -2.0, 3.5), 1e-10);
+            Assert.AreEqual(4.0, FindRoots.OfFunction(f1, 3.5, 5.5), 1e-10);
+            Assert.AreEqual(0.0, f1(FindRoots.OfFunction(f1, -2.0, 5.5, 1e-14)), 1e-14);
+        }
+
         void AssertComplexEqual(Complex expected, Complex actual, double delta)
         {
             Assert.AreEqual(expected.Real, actual.Real, delta);