Merge pull request #878 from kvedala/brents-method

feat: Brent's method to find extrema

DIRECTORY.md

@@ -131,6 +131,7 @@
 
 ## Numerical Methods
   * [Bisection Method](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/bisection_method.cpp)
+  * [Brent Method Extrema](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/brent_method_extrema.cpp)
   * [Durand Kerner Roots](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/durand_kerner_roots.cpp)
   * [False Position](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/false_position.cpp)
   * [Gaussian Elimination](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/gaussian_elimination.cpp)

numerical_methods/brent_method_extrema.cpp

@@ -0,0 +1,215 @@
+/**
+ * \file
+ * \brief Find real extrema of a univariate real function in a given interval
+ * using [Brent's method](https://en.wikipedia.org/wiki/Brent%27s_method).
+ *
+ * Refer the algorithm discoverer's publication
+ * [online](https://maths-people.anu.edu.au/~brent/pd/rpb011i.pdf) and also
+ * associated book:
+ * > R. P. Brent, Algorithms for Minimization without
+ * > Derivatives, Prentice-Hall, Englewood Cliffs, New Jersey, 1973
+ *
+ * \see golden_search_extrema.cpp
+ *
+ * \author [Krishna Vedala](https://github.com/kvedala)
+ */
+#define _USE_MATH_DEFINES  ///< required for MS Visual C++
+#include <cassert>
+#include <cmath>
+#include <functional>
+#include <iostream>
+#include <limits>
+
+#define EPSILON \
+    std::sqrt(  \
+        std::numeric_limits<double>::epsilon())  ///< system accuracy limit
+
+/**
+ * @brief Get the real root of a function in the given interval.
+ *
+ * @param f function to get root for
+ * @param lim_a lower limit of search window
+ * @param lim_b upper limit of search window
+ * @return root found in the interval
+ */
+double get_minima(const std::function<double(double)> &f, double lim_a,
+                  double lim_b) {
+    uint32_t iters = 0;
+
+    if (lim_a > lim_b) {
+        std::swap(lim_a, lim_b);
+    } else if (std::abs(lim_a - lim_b) <= EPSILON) {
+        std::cerr << "Search range must be greater than " << EPSILON << "\n";
+        return lim_a;
+    }
+
+    // golden ratio value
+    const double M_GOLDEN_RATIO = (3.f - std::sqrt(5.f)) / 2.f;
+
+    double v = lim_a + M_GOLDEN_RATIO * (lim_b - lim_a);
+    double u, w = v, x = v;
+    double fu, fv = f(v);
+    double fw = fv, fx = fv;
+
+    double mid_point = (lim_a + lim_b) / 2.f;
+    double p = 0, q = 0, r = 0;
+
+    double d, e = 0;
+    double tolerance, tolerance2;
+
+    do {
+        mid_point = (lim_a + lim_b) / 2.f;
+        tolerance = EPSILON * std::abs(x);
+        tolerance2 = 2 * tolerance;
+
+        if (std::abs(e) > tolerance2) {
+            // fit parabola
+            r = (x - w) * (fx - fv);
+            q = (x - v) * (fx - fw);
+            p = (x - v) * q - (x - w) * r;
+            q = 2.f * (q - r);
+            if (q > 0)
+                p = -p;
+            else
+                q = -q;
+            r = e;
+            e = d;
+        }
+
+        if (std::abs(p) < std::abs(0.5 * q * r) && p < q * (lim_b - x)) {
+            // parabolic interpolation step
+            d = p / q;
+            u = x + d;
+            if (u - lim_a < tolerance2 || lim_b - u < tolerance2)
+                d = x < mid_point ? tolerance : -tolerance;
+        } else {
+            // golden section interpolation step
+            e = (x < mid_point ? lim_b : lim_a) - x;
+            d = M_GOLDEN_RATIO * e;
+        }
+
+        // evaluate not too close to x
+        if (std::abs(d) >= tolerance)
+            u = d;
+        else if (d > 0)
+            u = tolerance;
+        else
+            u = -tolerance;
+        u += x;
+        fu = f(u);
+
+        // update variables
+        if (fu <= fx) {
+            if (u < x)
+                lim_b = x;
+            else
+                lim_a = x;
+            v = w;
+            fv = fw;
+            w = x;
+            fw = fx;
+            x = u;
+            fx = fu;
+        } else {
+            if (u < x)
+                lim_a = u;
+            else
+                lim_b = u;
+            if (fu <= fw || x == w) {
+                v = w;
+                fv = fw;
+                w = u;
+                fw = fu;
+            } else if (fu <= fv || v == x || v == w) {
+                v = u;
+                fv = fu;
+            }
+        }
+
+        iters++;
+    } while (std::abs(x - mid_point) > (tolerance - (lim_b - lim_a) / 2.f));
+
+    std::cout << " (iters: " << iters << ") ";
+
+    return x;
+}
+
+/**
+ * @brief Test function to find root for the function
+ * \f$f(x)= (x-2)^2\f$
+ * in the interval \f$[1,5]\f$
+ * \n Expected result = 2
+ */
+void test1() {
+    // define the function to minimize as a lambda function
+    std::function<double(double)> f1 = [](double x) {
+        return (x - 2) * (x - 2);
+    };
+
+    std::cout << "Test 1.... ";
+
+    double minima = get_minima(f1, -1, 5);
+
+    std::cout << minima << "...";
+
+    assert(std::abs(minima - 2) < EPSILON);
+    std::cout << "passed\n";
+}
+
+/**
+ * @brief Test function to find root for the function
+ * \f$f(x)= x^{\frac{1}{x}}\f$
+ * in the interval \f$[-2,10]\f$
+ * \n Expected result: \f$e\approx 2.71828182845904509\f$
+ */
+void test2() {
+    // define the function to maximize as a lambda function
+    // since we are maximixing, we negated the function return value
+    std::function<double(double)> func = [](double x) {
+        return -std::pow(x, 1.f / x);
+    };
+
+    std::cout << "Test 2.... ";
+
+    double minima = get_minima(func, -2, 5);
+
+    std::cout << minima << " (" << M_E << ")...";
+
+    assert(std::abs(minima - M_E) < EPSILON);
+    std::cout << "passed\n";
+}
+
+/**
+ * @brief Test function to find *maxima* for the function
+ * \f$f(x)= \cos x\f$
+ * in the interval \f$[0,12]\f$
+ * \n Expected result: \f$\pi\approx 3.14159265358979312\f$
+ */
+void test3() {
+    // define the function to maximize as a lambda function
+    // since we are maximixing, we negated the function return value
+    std::function<double(double)> func = [](double x) { return std::cos(x); };
+
+    std::cout << "Test 3.... ";
+
+    double minima = get_minima(func, -4, 12);
+
+    std::cout << minima << " (" << M_PI << ")...";
+
+    assert(std::abs(minima - M_PI) < EPSILON);
+    std::cout << "passed\n";
+}
+
+/** Main function */
+int main() {
+    std::cout.precision(18);
+
+    std::cout << "Computations performed with machine epsilon: " << EPSILON
+              << "\n";
+
+    test1();
+    test2();
+    test3();
+
+    return 0;
+}