From ddac0a0b615b62caf276ee029759424aa8af8646 Mon Sep 17 00:00:00 2001
From: Christoph Ruegg <git@cdrnet.ch>
Date: Wed, 1 May 2013 19:26:24 +0200
Subject: [PATCH] RootFinding: algorithm cosmetics

---
 src/Numerics/RootFinding/FindRoots.cs | 55 +++++++++++++++------------
 1 file changed, 30 insertions(+), 25 deletions(-)

diff --git a/src/Numerics/RootFinding/FindRoots.cs b/src/Numerics/RootFinding/FindRoots.cs
index adf93d84..5604cabf 100644
--- a/src/Numerics/RootFinding/FindRoots.cs
+++ b/src/Numerics/RootFinding/FindRoots.cs
@@ -18,11 +18,9 @@ namespace MathNet.Numerics.RootFinding
         /// </remarks>
 	    public static double BrentMethod(Func<double, double> f, double xmin, double xmax, double accuracy = 1e-8, int maxIterations = 100)
         {
-            double min1, min2;
-            double p, q, r, s, xAcc1, xMid = 0;
-            double d = 0.0, e = 0.0;
+	        double xMid = 0;
+	        double d = 0.0, e = 0.0;
 
-            // set up
             double fxmin = f(xmin);
             double fxmax = f(xmax);
             double root = xmax;
@@ -30,13 +28,14 @@ namespace MathNet.Numerics.RootFinding
 
             for (int i = 0; i <= maxIterations; i++)
             {
+                // adjust bounds
                 if (Math.Sign(froot) == Math.Sign(fxmax))
                 {
-                    // Rename xMin_, root_, xMax_ and adjust bounds
                     xmax = xmin;
                     fxmax = fxmin;
                     e = d = root - xmin;
                 }
+
                 if (Math.Abs(fxmax) < Math.Abs(froot))
                 {
                     xmin = root;
@@ -46,20 +45,22 @@ namespace MathNet.Numerics.RootFinding
                     froot = fxmax;
                     fxmax = fxmin;
                 }
-                // Convergence check
-                xAcc1 = 2.0 * Precision.DoubleMachinePrecision * Math.Abs(root) + 0.5 * accuracy;
+
+                // convergence check
+                double xAcc1 = 2.0 * Precision.DoubleMachinePrecision * Math.Abs(root) + 0.5 * accuracy;
                 xMid = (xmax - root) / 2.0;
-                if (Math.Abs(xMid) <= xAcc1 || Close(froot, 0.0))
+                if (Math.Abs(xMid) <= xAcc1 || froot.AlmostEqualWithAbsoluteError(0, froot, accuracy))
                 {
                     return root;
                 }
-                if (Math.Abs(e) >= xAcc1 &&
-                    Math.Abs(fxmin) > Math.Abs(froot))
-                {
 
+                if (Math.Abs(e) >= xAcc1 && Math.Abs(fxmin) > Math.Abs(froot))
+                {
                     // Attempt inverse quadratic interpolation
-                    s = froot / fxmin;
-                    if (Close(xmin, xmax))
+                    double s = froot / fxmin;
+                    double p;
+                    double q;
+                    if (xmin.AlmostEqual(xmax))
                     {
                         p = 2.0 * xMid * s;
                         q = 1.0 - s;
@@ -67,22 +68,27 @@ namespace MathNet.Numerics.RootFinding
                     else
                     {
                         q = fxmin / fxmax;
-                        r = froot / 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) q = -q;  // Check whether in bounds
+
+                    if (p > 0.0)
+                    {
+                        // Check whether in bounds
+                        q = -q;
+                    }
                     p = Math.Abs(p);
-                    min1 = 3.0 * xMid * q - Math.Abs(xAcc1 * q);
-                    min2 = Math.Abs(e * q);
-                    if (2.0 * p < Math.Min(min1, min2))
+                    if (2.0 * p < Math.Min(3.0 * xMid * q - Math.Abs(xAcc1 * q), Math.Abs(e * q)))
                     {
-                        e = d;                // Accept interpolation
+                        // Accept interpolation
+                        e = d;
                         d = p / q;
                     }
                     else
                     {
-                        d = xMid;  // Interpolation failed, use bisection
+                        // Interpolation failed, use bisection
+                        d = xMid;
                         e = d;
                     }
                 }
@@ -95,9 +101,13 @@ namespace MathNet.Numerics.RootFinding
                 xmin = root;
                 fxmin = froot;
                 if (Math.Abs(d) > xAcc1)
+                {
                     root += d;
+                }
                 else
+                {
                     root += Sign(xAcc1, xMid);
+                }
                 froot = f(root);
             }
 
@@ -111,10 +121,5 @@ namespace MathNet.Numerics.RootFinding
         {
             return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);
         }
-
-        static bool Close(double d1, double d2)
-        {
-            return Math.Abs(d1 - d2) <= double.Epsilon;
-        }
 	}
 }