.

utils/xalgorithm.h

@@ -2,9 +2,14 @@
 
 #include <algorithm>
 #include <cmath>
+#include <functional>
 
 template <typename T>
-double golden_section_search(const T& f, double a, double b, double tol)
+double golden_section_search(
+    const T& f, 
+    double a, 
+    double b, 
+    const double time_accuracy)
 {
     const double invphi = (std::sqrt(5.0) - 1) / 2.0;
     const double invphi2 = (3 - std::sqrt(5.0)) / 2.0;
@@ -14,11 +19,11 @@ double golden_section_search(const T& f, double a, double b, double tol)
 
     double h = b - a;
 
-    if (h <= tol) {
+    if (h <= time_accuracy) {
         return (a + b) / 2.0;
     }
 
-    const int n = (int)ceil(std::log(tol / h) / std::log(invphi));
+    const int n = static_cast<int>(ceil(std::log(time_accuracy / h) / std::log(invphi)));
 
     double c = a + invphi2 * h;
     double d = a + invphi * h;
@@ -51,3 +56,57 @@ double golden_section_search(const T& f, double a, double b, double tol)
         return (c + b) / 2.0;
     }
 }
+
+inline double golden_section_search_stdf(
+    double a, 
+    double b, 
+    const double time_accuracy,
+    const std::function<double(const double&)>& f
+    )
+{
+    const double invphi = (std::sqrt(5.0) - 1) / 2.0;
+    const double invphi2 = (3 - std::sqrt(5.0)) / 2.0;
+
+    a = std::min(a, b);
+    b = std::max(a, b);
+
+    double h = b - a;
+
+    if (h <= time_accuracy) {
+        return (a + b) / 2.0;
+    }
+
+    const int n = static_cast<int>(ceil(std::log(time_accuracy / h) / std::log(invphi)));
+
+    double c = a + invphi2 * h;
+    double d = a + invphi * h;
+    double yc = f(c);
+    double yd = f(d);
+
+    for (int k = 0; k < n - 1; k++) {
+        if (yc < yd) {
+            b = d;
+            d = c;
+            yd = yc;
+            h = invphi * h;
+            c = a + invphi2 * h;
+            yc = f(c);
+        }
+        else {
+            a = c;
+            c = d;
+            yc = yd;
+            h = invphi * h;
+            d = a + invphi * h;
+            yd = f(d);
+        }
+    }
+
+    if (yc < yd) {
+        return (a + d) / 2.0;
+    }
+    else {
+        return (c + b) / 2.0;
+    }
+}
+