Merge brent root finding algorithm

src/Numerics/Numerics.csproj

@@ -109,6 +109,9 @@
     <Compile Include="Financial\AbsoluteRiskMeasures.cs" />
     <Compile Include="LinearAlgebra\Generic\Matrix.BCL.cs" />
     <Compile Include="LinearAlgebra\Generic\Vector.BCL.cs" />
+    <Compile Include="RootFinding\Bracketing.cs" />
+    <Compile Include="RootFinding\FindRoots.cs" />
+    <Compile Include="RootFinding\RootFindingException.cs" />
     <Compile Include="SpecialFunctions\Evaluate.cs" />
     <Compile Include="SpecialFunctions\ModifiedStruve.cs" />
     <Compile Include="SpecialFunctions\ModifiedBessel.cs" />

src/Numerics/Properties/Resources.Designer.cs

@@ -60,6 +60,15 @@ namespace MathNet.Numerics.Properties {
             }
         }
         
+        /// <summary>
+        ///   Looks up a localized string similar to The accuracy couldn&apos;t be reached with the specified number of iterations..
+        /// </summary>
+        public static string AccuracyNotReached {
+            get {
+                return ResourceManager.GetString("AccuracyNotReached", resourceCulture);
+            }
+        }
+        
         /// <summary>
         ///   Looks up a localized string similar to The array arguments must have the same length..
         /// </summary>
@@ -744,6 +753,15 @@ namespace MathNet.Numerics.Properties {
             }
         }
         
+        /// <summary>
+        ///   Looks up a localized string similar to The algorithm ended without root in the range..
+        /// </summary>
+        public static string RootNotFound {
+            get {
+                return ResourceManager.GetString("RootNotFound", resourceCulture);
+            }
+        }
+        
         /// <summary>
         ///   Looks up a localized string similar to The number of rows must greater than or equal to the number of columns..
         /// </summary>

src/Numerics/Properties/Resources.resx

@@ -369,10 +369,16 @@
   <data name="TooManyElements" xml:space="preserve">
     <value>We only support sparse matrix with less than int.MaxValue elements.</value>
   </data>
-  <data name="Interpolation_Initialize_SamplePointsNotStrictlyAscendingOrder">
-    <value xml:space="preserve">Sample points should be sorted in strictly ascending order</value>
+  <data name="Interpolation_Initialize_SamplePointsNotStrictlyAscendingOrder" xml:space="preserve">
+    <value>Sample points should be sorted in strictly ascending order</value>
   </data>
-  <data name="Interpolation_Initialize_SamplePointsNotUnique">
-    <value xml:space="preserve">All sample points should be unique.</value>
+  <data name="Interpolation_Initialize_SamplePointsNotUnique" xml:space="preserve">
+    <value>All sample points should be unique.</value>
+  </data>
+  <data name="AccuracyNotReached" xml:space="preserve">
+    <value>The accuracy couldn't be reached with the specified number of iterations.</value>
+  </data>
+  <data name="RootNotFound" xml:space="preserve">
+    <value>The algorithm ended without root in the range.</value>
   </data>
 </root>

src/Numerics/RootFinding/Bracketing.cs

@@ -0,0 +1,48 @@
+using System;
+using MathNet.Numerics.Properties;
+
+namespace MathNet.Numerics.RootFinding
+{
+    public static class Bracketing
+    {
+        /// <summary>Detect a range containing at least one root.</summary>
+        /// <param name="f">The function to detect roots from.</param>
+        /// <param name="xmin">Lower value of the range.</param>
+        /// <param name="xmax">Upper value of the range</param>
+        /// <param name="factor">The growing factor of research. Usually 1.6.</param>
+        /// <param name="maxIterations">Maximum number of iterations. Usually 50.</param>
+        /// <returns>True if the bracketing operation succeeded, false otherwise.</returns>
+        /// <remarks>This iterative methods stops when two values with opposite signs are found.</remarks>
+        public static bool SearchOutward(Func<double, double> f, ref double xmin, ref double xmax, double factor = 1.6, int maxIterations = 50)
+        {
+            if (xmin >= xmax)
+            {
+                throw new ArgumentOutOfRangeException("xmax", string.Format(Resources.ArgumentOutOfRangeGreater, "xmax", "xmin"));
+            }
+
+            double fmin = f(xmin);
+            double fmax = f(xmax);
+
+            for(int i=0;i<maxIterations; i++)
+            {
+                if (Math.Sign(fmin) != Math.Sign(fmax))
+                {
+                    return true;
+                }
+
+                if (Math.Abs(fmin) < Math.Abs(fmax))
+                {
+                    xmin += factor * (xmin - xmax);
+                    fmin = f(xmin);
+                }
+                else
+                {
+                    xmax += factor * (xmax - xmin);
+                    fmax = f(xmax);
+                }
+            }
+
+            return false;
+        }
+    }
+}

src/Numerics/RootFinding/FindRoots.cs

@@ -0,0 +1,125 @@
+using System;
+using MathNet.Numerics.Properties;
+
+namespace MathNet.Numerics.RootFinding
+{
+	public static class FindRoots
+	{
+	    /// <summary>Find a solution of the equation f(x)=0.</summary>
+	    /// <param name="f">The function to find roots from.</param>
+	    /// <param name="xmin">The low value of the range where the root is supposed to be.</param>
+	    /// <param name="xmax">The high value of the range where the root is supposed to be.</param>
+	    /// <param name="accuracy">Desired accuracy. The root will be refined until the accuracy or the maximum number of iterations is reached.</param>
+	    /// <param name="maxIterations">Maximum number of iterations. Usually 100.</param>
+	    /// <returns>Returns the root with the specified accuracy.</returns>
+        /// <remarks>
+        /// Algorithm by by Brent, Van Wijngaarden, Dekker et al.
+        /// Implementation inspired by Press, Teukolsky, Vetterling, and Flannery, "Numerical Recipes in C", 2nd edition, Cambridge University Press
+        /// </remarks>
+	    public static double BrentMethod(Func<double, double> f, double xmin, double xmax, double accuracy = 1e-8, int maxIterations = 100)
+        {
+	        double xMid = 0;
+	        double d = 0.0, e = 0.0;
+
+            double fxmin = f(xmin);
+            double fxmax = f(xmax);
+            double root = xmax;
+            double froot = fxmax;
+
+            for (int i = 0; i <= maxIterations; i++)
+            {
+                // adjust bounds
+                if (Math.Sign(froot) == Math.Sign(fxmax))
+                {
+                    xmax = xmin;
+                    fxmax = fxmin;
+                    e = d = root - xmin;
+                }
+
+                if (Math.Abs(fxmax) < Math.Abs(froot))
+                {
+                    xmin = root;
+                    root = xmax;
+                    xmax = xmin;
+                    fxmin = froot;
+                    froot = fxmax;
+                    fxmax = fxmin;
+                }
+
+                // convergence check
+                double xAcc1 = 2.0 * Precision.DoubleMachinePrecision * Math.Abs(root) + 0.5 * accuracy;
+                xMid = (xmax - root) / 2.0;
+                if (Math.Abs(xMid) <= xAcc1 || froot.AlmostEqualWithAbsoluteError(0, froot, accuracy))
+                {
+                    return root;
+                }
+
+                if (Math.Abs(e) >= xAcc1 && Math.Abs(fxmin) > Math.Abs(froot))
+                {
+                    // Attempt inverse quadratic interpolation
+                    double s = froot / fxmin;
+                    double p;
+                    double q;
+                    if (xmin.AlmostEqual(xmax))
+                    {
+                        p = 2.0 * xMid * s;
+                        q = 1.0 - s;
+                    }
+                    else
+                    {
+                        q = fxmin / fxmax;
+                        double r = froot / fxmax;
+                        p = s * (2.0 * xMid * q * (q - r) - (root - xmin) * (r - 1.0));
+                        q = (q - 1.0) * (r - 1.0) * (s - 1.0);
+                    }
+
+                    if (p > 0.0)
+                    {
+                        // Check whether in bounds
+                        q = -q;
+                    }
+                    p = Math.Abs(p);
+                    if (2.0 * p < Math.Min(3.0 * xMid * q - Math.Abs(xAcc1 * q), Math.Abs(e * q)))
+                    {
+                        // Accept interpolation
+                        e = d;
+                        d = p / q;
+                    }
+                    else
+                    {
+                        // Interpolation failed, use bisection
+                        d = xMid;
+                        e = d;
+                    }
+                }
+                else
+                {
+                    // Bounds decreasing too slowly, use bisection
+                    d = xMid;
+                    e = d;
+                }
+                xmin = root;
+                fxmin = froot;
+                if (Math.Abs(d) > xAcc1)
+                {
+                    root += d;
+                }
+                else
+                {
+                    root += Sign(xAcc1, xMid);
+                }
+                froot = f(root);
+            }
+
+            // The algorithm has exceeded the number of iterations allowed
+            throw new RootFindingException(Resources.AccuracyNotReached, maxIterations, xmin, xmax, Math.Abs(xMid));
+        }
+
+        /// <summary>Helper method useful for preventing rounding errors.</summary>
+        /// <returns>a*sign(b)</returns>
+        static double Sign(double a, double b)
+        {
+            return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);
+        }
+	}
+}

src/Numerics/RootFinding/RootFindingException.cs

@@ -0,0 +1,21 @@
+using System;
+
+namespace MathNet.Numerics.RootFinding
+{
+    public class RootFindingException : Exception
+    {
+        public RootFindingException(string message, int iteration, double rangeMin, double rangeMax, double accuracy)
+            : base(message)
+        {
+            Iteration = iteration;
+            RangeMin = rangeMin;
+            RangeMax = rangeMax;
+            Accuracy = accuracy;
+        }
+
+        public int Iteration { get; set; }
+        public double RangeMin { get; set; }
+        public double RangeMax { get; set; }
+        public double Accuracy { set; get; }
+    }
+}

src/UnitTests/RootFindingTests/BrentTest.cs

@@ -0,0 +1,16 @@
+using MathNet.Numerics.RootFinding;
+using NUnit.Framework;
+
+namespace MathNet.Numerics.UnitTests.RootFindingTests
+{
+    [TestFixture]
+    public class BrentTest
+    {
+        [Test]
+        public void MultipleRoots()
+        {
+            double root = FindRoots.BrentMethod(x => x*x - 4, -5, 5, 1e-14, 100);
+            Assert.AreEqual(0, root*root - 4);
+        }
+    }
+}

src/UnitTests/UnitTests.csproj

@@ -770,6 +770,7 @@
     <Compile Include="Random\WH1982Tests.cs" />
     <Compile Include="Random\WH2006Tests.cs" />
     <Compile Include="Random\XorshiftTests.cs" />
+    <Compile Include="RootFindingTests\BrentTest.cs" />
     <Compile Include="SortingTests.cs" />
     <Compile Include="SpecialFunctionsTests\ModifiedStruveTests.cs" />
     <Compile Include="SpecialFunctionsTests\ModifiedBesselTests.cs" />