From c7b758f53654ee3c499d4e6140ff7b2bd8cadcf9 Mon Sep 17 00:00:00 2001
From: Gabriel Ozouf <gabriel.ozouf@numworks.com>
Date: Wed, 9 Dec 2020 15:50:01 +0100
Subject: [PATCH] [poincare/expression] Fix solutions to e^x=0

Because of the limitations of the floating-point representation, e^x is
null for x <= 710, causing the nextRoot functions to find roots for it.
To prevent this, we looking for places where a function changes sign, we
actually require the function to take two non-null values of different
signs.
---
 apps/solver/test/equation_store.cpp |  2 ++
 poincare/src/expression.cpp         | 24 +++++++++++++++++++++---
 poincare/test/function_solver.cpp   |  8 +++++++-
 3 files changed, 30 insertions(+), 4 deletions(-)

diff --git a/apps/solver/test/equation_store.cpp b/apps/solver/test/equation_store.cpp
index 608f81f6f..8ffbd0c54 100644
--- a/apps/solver/test/equation_store.cpp
+++ b/apps/solver/test/equation_store.cpp
@@ -82,6 +82,8 @@ QUIZ_CASE(equation_solve) {
   assert_solves_numerically_to("cos(x)=0", -100, 100, {-90.0, 90.0});
   assert_solves_numerically_to("cos(x)=0", -900, 1000, {-810.0, -630.0, -450.0, -270.0, -90.0, 90.0, 270.0, 450.0, 630.0, 810.0});
   assert_solves_numerically_to("√(y)=0", -900, 1000, {0}, "y");
+  assert_solves_numerically_to("ℯ^x=0", -1000, 1000, {});
+  assert_solves_numerically_to("ℯ^x/1000=0", -1000, 1000, {});
 
   // Long variable names
   assert_solves_to("2abcde+3=4", "abcde=1/2");
diff --git a/poincare/src/expression.cpp b/poincare/src/expression.cpp
index 685242180..498e3da38 100644
--- a/poincare/src/expression.cpp
+++ b/poincare/src/expression.cpp
@@ -1005,6 +1005,9 @@ Coordinate2D<double> Expression::nextMaximum(const char * symbol, double start,
 }
 
 double Expression::nextRoot(const char * symbol, double start, double step, double max, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const {
+  if (nullStatus(context) == ExpressionNode::NullStatus::Null) {
+    return start + step;
+  }
   return nextIntersectionWithExpression(symbol, start, step, max,
       [](double x, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, const void * context1, const void * context2, const void * context3) {
         const Expression * expression0 = reinterpret_cast<const Expression *>(context1);
@@ -1139,17 +1142,32 @@ double Expression::nextIntersectionWithExpression(const char * symbol, double st
 
 void Expression::bracketRoot(const char * symbol, double start, double step, double max, double result[2], Solver::ValueAtAbscissa evaluation, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, const Expression expression) const {
   double a = start;
+  double c = a;
   double b = start+step;
+  double fa = evaluation(a, context, complexFormat, angleUnit, this, symbol, &expression);
+  double fc = fa;
+  double fb = evaluation(b, context, complexFormat, angleUnit, this, symbol, &expression);
   while (step > 0.0 ? b <= max : b >= max) {
-    double fa = evaluation(a, context, complexFormat, angleUnit, this, symbol, &expression);
-    double fb = evaluation(b, context, complexFormat, angleUnit, this, symbol, &expression);
-    if (fa*fb <= 0) {
+    if (fa == 0. && ((fc < 0. && fb > 0.) || (fc > 0. && fb < 0.))) {
+      /* If fa is null, we still check that the function changes sign on ]c,b[,
+       * and that fc and fb are not null. Otherwise, it's more likely those
+       * zeroes are caused by approximation errors. */
+      result[0] = a;
+      result[1] = b;
+      return;
+    } else if (fb != 0. && ((fa < 0.) != (fb < 0.))) {
+      /* The function changes sign.
+       * The case fb = 0 is handled in the next pass with fa = 0. */
       result[0] = a;
       result[1] = b;
       return;
     }
+    c = a;
+    fc = fa;
     a = b;
+    fa = fb;
     b = b+step;
+    fb = evaluation(b, context, complexFormat, angleUnit, this, symbol, &expression);
   }
   result[0] = NAN;
   result[1] = NAN;
diff --git a/poincare/test/function_solver.cpp b/poincare/test/function_solver.cpp
index 4dee3c65f..11274a2b0 100644
--- a/poincare/test/function_solver.cpp
+++ b/poincare/test/function_solver.cpp
@@ -236,9 +236,15 @@ QUIZ_CASE(poincare_function_root) {
   {
     constexpr int numberOfRoots = 1;
     Coordinate2D<double> roots[numberOfRoots] = {
-      Coordinate2D<double>(99.8, 0.0)};
+      Coordinate2D<double>(99.9, 0.0)};
     assert_points_of_interest_are(PointOfInterestType::Root, numberOfRoots, roots, "0", nullptr, "a", 100.0, -0.1, -1.0);
   }
+  {
+    constexpr int numberOfRoots = 1;
+    Coordinate2D<double> roots[numberOfRoots] = {
+      Coordinate2D<double>(NAN, 0.0)};
+    assert_points_of_interest_are(PointOfInterestType::Root, numberOfRoots, roots, "ℯ^x", nullptr, "a", -1000.0, 0.1, -800);
+  }
 }
 
 QUIZ_CASE(poincare_function_intersection) {