From 03ebffa09d1bdfd81ff6fce672157e80f789bb4a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C3=89milie=20Feral?= <emilie.feral@numworks.com>
Date: Fri, 12 Jan 2018 10:49:42 +0100
Subject: [PATCH] [apps] Graph: first version of a function minimum finding
 algorithm

---
 apps/graph/cartesian_function.cpp | 96 +++++++++++++++++++++++++++++++
 apps/graph/cartesian_function.h   |  9 +++
 2 files changed, 105 insertions(+)

diff --git a/apps/graph/cartesian_function.cpp b/apps/graph/cartesian_function.cpp
index 15f749e07..ebd958ca5 100644
--- a/apps/graph/cartesian_function.cpp
+++ b/apps/graph/cartesian_function.cpp
@@ -1,4 +1,8 @@
 #include "cartesian_function.h"
+#include <float.h>
+#include <cmath>
+
+using namespace Poincare;
 
 namespace Graph {
 
@@ -37,8 +41,100 @@ double CartesianFunction::sumBetweenBounds(double start, double end, Poincare::C
   return integral.approximateToScalar<double>(*context);
 }
 
+CartesianFunction::Point CartesianFunction::mininimumBetweenBounds(double start, double end, Poincare::Context * context) const {
+  Point p = brentAlgorithm(start, end, context);
+  if (evaluateAtAbscissa(p.abscissa-k_sqrtEps, context) < p.value || evaluateAtAbscissa(p.abscissa+k_sqrtEps, context) < p.value || std::isnan(p.value)) {
+    p.abscissa = NAN;
+    p.value = NAN;
+  }
+  return p;
+}
+
 char CartesianFunction::symbol() const {
   return 'x';
 }
 
+CartesianFunction::Point CartesianFunction::brentAlgorithm(double ax, double bx, Context * context) const {
+  double e = 0.0;
+  double a = ax;
+  double b = bx;
+  double x = a+k_goldenRatio*(b-a);
+  double v = x;
+  double w = x;
+  double fx = evaluateAtAbscissa(x, context);
+  double fw = fx;
+  double fv = fw;
+
+  double d = NAN;
+  double u, fu;
+
+  for (int i = 0; i < 100; i++) {
+    double m = 0.5*(a+b);
+    double tol1 = k_sqrtEps*std::fabs(x)+1E-10;
+    double tol2 = 2.0*tol1;
+    if (std::fabs(x-m) <= tol2-0.5*(b-a)) {
+      Point result = {.abscissa = x, .value = fx};
+      return result;
+    }
+    double p = 0;
+    double q = 0;
+    double r = 0;
+    if (std::fabs(e) > tol1) {
+      r = (x-w)*(fx-fv);
+      q = (x-v)*(fx-fw);
+      p = (x-v)*q -(x-w)*r;
+      q = 2.0*(q-r);
+      if (q>0.0) {
+        p = -p;
+      } else {
+        q = -q;
+      }
+      r = e;
+      e = d;
+    }
+    if (std::fabs(p) < std::fabs(0.5*q*r) && p<q*(a-x) && p<q*(b-x)) {
+      d = p/q;
+      u= x+d;
+      if (u-a < tol2 || b-u < tol2) {
+        d = x < m ? tol1 : -tol1;
+      }
+    } else {
+      e = x<m ? b-x : a-x;
+      d = k_goldenRatio*e;
+    }
+    u = x + (std::fabs(d) >= tol1 ? d : (d>0 ? tol1 : -tol1));
+    fu = evaluateAtAbscissa(u, context);
+    if (fu <= fx) {
+      if (u<x) {
+        b = x;
+      } else {
+        a = x;
+      }
+      v = w;
+      fv = fw;
+      w = x;
+      fw = fx;
+      x = u;
+      fx = fu;
+    } else {
+      if (u<x) {
+        a = u;
+      } else {
+        b = u;
+      }
+      if (fu <= fw || w == x) {
+        v = w;
+        fv = fw;
+        w = u;
+        fw = fu;
+      } else if (fu <= fv || v == x || v == w) {
+        v = u;
+        fv = fu;
+      }
+    }
+  }
+  Point result = {.abscissa = NAN, .value = NAN};
+  return result;
+}
+
 }
diff --git a/apps/graph/cartesian_function.h b/apps/graph/cartesian_function.h
index 9d4ed6f28..9d5a9b91f 100644
--- a/apps/graph/cartesian_function.h
+++ b/apps/graph/cartesian_function.h
@@ -13,8 +13,17 @@ public:
   void setDisplayDerivative(bool display);
   double approximateDerivative(double x, Poincare::Context * context) const;
   double sumBetweenBounds(double start, double end, Poincare::Context * context) const override;
+  struct Point {
+    double abscissa;
+    double value;
+  };
+  Point mininimumBetweenBounds(double start, double end, Poincare::Context * context) const;
   char symbol() const override;
 private:
+  constexpr static double k_sqrtEps = 1.4901161193847656E-8; // sqrt(DBL_EPSILON)
+  constexpr static double k_goldenRatio = 0.381966011250105151795413165634361882279690820194237137864; // (3-sqrt(5))/2
+  double bracketFunction(double x, double z, Poincare::Context * context);
+  Point brentAlgorithm(double ax, double bx, Poincare::Context * context) const;
   bool m_displayDerivative;
 };