[apps/shared] New automatic zoom on curves

Initial zoom for displaying curves is computed according to the
following rules :
  - For polar and parametric curves, the algorithm has not changed.
    Vertical and horizontal ranges are chosen in order to display the
    full curve with orthonormal graduations.
  - For cartesian curves, the horizontal range is chosen in order to
    display all points of interest (roots, extrema, asymptotes...).
    The vertical range is fitted to be able to display those points, but
    also cuts out points outside of the function's order of magnitude.

Change-Id: Idf8233fc2e6586b85d34c4da152c83e75513d85c

apps/graph/graph/graph_controller.cpp

@@ -35,116 +35,40 @@ void GraphController::viewWillAppear() {
   selectFunctionWithCursor(indexFunctionSelectedByCursor());
 }
 
-bool GraphController::defautRangeIsNormalized() const {
+bool GraphController::defaultRangeIsNormalized() const {
   return functionStore()->displaysNonCartesianFunctions();
 }
 
-void GraphController::interestingFunctionRange(ExpiringPointer<ContinuousFunction> f, float tMin, float tMax, float step, float * xm, float * xM, float * ym, float * yM) const {
-  Poincare::Context * context = textFieldDelegateApp()->localContext();
-  const int balancedBound = std::floor((tMax-tMin)/2/step);
-  for (int j = -balancedBound; j <= balancedBound ; j++) {
-    float t = (tMin+tMax)/2 + step * j;
-    Coordinate2D<float> xy = f->evaluateXYAtParameter(t, context);
-    float x = xy.x1();
-    float y = xy.x2();
-    if (!std::isnan(x) && !std::isinf(x) && !std::isnan(y) && !std::isinf(y)) {
-      *xm = std::min(*xm, x);
-      *xM = std::max(*xM, x);
-      *ym = std::min(*ym, y);
-      *yM = std::max(*yM, y);
-    }
-  }
+void GraphController::interestingRanges(float * xm, float * xM, float * ym, float * yM) const {
+  privateComputeRanges(true, xm, xM, ym, yM);
 }
 
-void GraphController::interestingRanges(float * xm, float * xM, float * ym, float * yM) const {
-  float resultxMin = FLT_MAX;
-  float resultxMax = -FLT_MAX;
-  float resultyMin = FLT_MAX;
-  float resultyMax = -FLT_MAX;
+Shared::InteractiveCurveViewRangeDelegate::Range GraphController::computeYRange(Shared::InteractiveCurveViewRange * interactiveCurveViewRange) {
+  float xm = interactiveCurveViewRange->xMin(),
+        xM = interactiveCurveViewRange->xMax(),
+        ym = FLT_MAX,
+        yM = -FLT_MAX;
+  privateComputeRanges(false, &xm, &xM, &ym, &yM);
+  return Shared::InteractiveCurveViewRangeDelegate::Range{.min = ym, .max = yM};
+}
+
+void GraphController::privateComputeRanges(bool tuneXRange, float * xm, float * xM, float * ym, float * yM) const {
+  Poincare::Context * context = textFieldDelegateApp()->localContext();
+  float resultXMin = tuneXRange ? FLT_MAX : *xm;
+  float resultXMax = tuneXRange ? -FLT_MAX : *xM;
+  float resultYMin = FLT_MAX;
+  float resultYMax = -FLT_MAX;
   assert(functionStore()->numberOfActiveFunctions() > 0);
-  int functionsCount = 0;
-  if (functionStore()->displaysNonCartesianFunctions(&functionsCount)) {
-    for (int i = 0; i < functionsCount; i++) {
-      ExpiringPointer<ContinuousFunction> f = functionStore()->modelForRecord(functionStore()->activeRecordAtIndex(i));
-      if (f->plotType() == ContinuousFunction::PlotType::Cartesian) {
-        continue;
-      }
-      /* Scan x-range from the middle to the extrema in order to get balanced
-       * y-range for even functions (y = 1/x). */
-      double tMin = f->tMin();
-      double tMax = f->tMax();
-      assert(!std::isnan(tMin));
-      assert(!std::isnan(tMax));
-      assert(!std::isnan(f->rangeStep()));
-      interestingFunctionRange(f, tMin, tMax, f->rangeStep(), &resultxMin, &resultxMax, &resultyMin, &resultyMax);
-    }
-    if (resultxMin > resultxMax) {
-      resultxMin = - Range1D::k_default;
-      resultxMax = Range1D::k_default;
-    }
-  } else {
-    resultxMin = const_cast<GraphController *>(this)->interactiveCurveViewRange()->xMin();
-    resultxMax = const_cast<GraphController *>(this)->interactiveCurveViewRange()->xMax();
-  }
-  /* In practice, a step smaller than a pixel's width is needed for sampling
-   * the values of a function. Otherwise some relevant extremal values may be
-   * missed. */
+  int functionsCount = functionStore()->numberOfActiveFunctions();
   for (int i = 0; i < functionsCount; i++) {
     ExpiringPointer<ContinuousFunction> f = functionStore()->modelForRecord(functionStore()->activeRecordAtIndex(i));
-    if (f->plotType() != ContinuousFunction::PlotType::Cartesian) {
-      continue;
-    }
-    /* Scan x-range from the middle to the extrema in order to get balanced
-     * y-range for even functions (y = 1/x). */
-    assert(!std::isnan(f->tMin()));
-    assert(!std::isnan(f->tMax()));
-    const double tMin = std::max(f->tMin(), resultxMin);
-    const double tMax = std::min(f->tMax(), resultxMax);
-    const double step = (tMax - tMin) / (2.0 * (m_view.bounds().width() - 1.0));
-    interestingFunctionRange(f, tMin, tMax, step, &resultxMin, &resultxMax, &resultyMin, &resultyMax);
-  }
-  if (resultyMin > resultyMax) {
-    resultyMin = - Range1D::k_default;
-    resultyMax = Range1D::k_default;
+    f->rangeForDisplay(&resultXMin, &resultXMax, &resultYMin, &resultYMax, context, tuneXRange);
   }
 
-  *xm = resultxMin;
-  *xM = resultxMax;
-  *ym = resultyMin;
-  *yM = resultyMax;
-}
-
-float GraphController::interestingXHalfRange() const {
-  float characteristicRange = 0.0f;
-  Poincare::Context * context = textFieldDelegateApp()->localContext();
-  ContinuousFunctionStore * store = functionStore();
-  int nbActiveFunctions = store->numberOfActiveFunctions();
-  double tMin = INFINITY;
-  double tMax = -INFINITY;
-  for (int i = 0; i < nbActiveFunctions; i++) {
-    ExpiringPointer<ContinuousFunction> f = store->modelForRecord(store->activeRecordAtIndex(i));
-    float fRange = f->expressionReduced(context).characteristicXRange(context, Poincare::Preferences::sharedPreferences()->angleUnit());
-    if (!std::isnan(fRange) && !std::isinf(fRange)) {
-      characteristicRange = std::max(fRange, characteristicRange);
-    }
-    // Compute the combined range of the functions
-    assert(f->plotType() == ContinuousFunction::PlotType::Cartesian); // So that tMin tMax represents xMin xMax
-    tMin = std::min<double>(tMin, f->tMin());
-    tMax = std::max<double>(tMax, f->tMax());
-  }
-  constexpr float rangeMultiplicator = 1.6f;
-  if (characteristicRange > 0.0f ) {
-    return rangeMultiplicator * characteristicRange;
-  }
-  float defaultXHalfRange = InteractiveCurveViewRangeDelegate::interestingXHalfRange();
-  assert(tMin <= tMax);
-  if (tMin >= -defaultXHalfRange && tMax <= defaultXHalfRange) {
-    /* If the combined Range of the functions is smaller than the default range,
-     * use it. */
-    float f = rangeMultiplicator * (float)std::max(std::fabs(tMin), std::fabs(tMax));
-    return (std::isnan(f) || std::isinf(f)) ? defaultXHalfRange : f;
-  }
-  return defaultXHalfRange;
+  *xm = resultXMin;
+  *xM = resultXMax;
+  *ym = resultYMin;
+  *yM = resultYMax;
 }
 
 void GraphController::selectFunctionWithCursor(int functionIndex) {

apps/graph/graph/graph_controller.h

@@ -20,7 +20,6 @@ public:
   void viewWillAppear() override;
   bool displayDerivativeInBanner() const { return m_displayDerivativeInBanner; }
   void setDisplayDerivativeInBanner(bool displayDerivative) { m_displayDerivativeInBanner = displayDerivative; }
-  float interestingXHalfRange() const override;
   void interestingRanges(float * xm, float * xM, float * ym, float * yM) const override;
 private:
   int estimatedBannerNumberOfLines() const override { return 1 + m_displayDerivativeInBanner; }
@@ -34,10 +33,12 @@ private:
   GraphView * functionGraphView() override { return &m_view; }
   CurveParameterController * curveParameterController() override { return &m_curveParameterController; }
   ContinuousFunctionStore * functionStore() const override { return static_cast<ContinuousFunctionStore *>(Shared::FunctionGraphController::functionStore()); }
-  bool defautRangeIsNormalized() const override;
+  bool defaultRangeIsNormalized() const override;
   void interestingFunctionRange(Shared::ExpiringPointer<Shared::ContinuousFunction> f, float tMin, float tMax, float step, float * xm, float * xM, float * ym, float * yM) const;
   bool shouldSetDefaultOnModelChange() const override;
   void jumpToLeftRightCurve(double t, int direction, int functionsCount, Ion::Storage::Record record) override;
+  Range computeYRange(Shared::InteractiveCurveViewRange * interactiveCurveViewRange) override;
+  void privateComputeRanges(bool tuneXRange, float * xm, float * xM, float * ym, float * yM) const;
 
   Shared::RoundCursorView m_cursorView;
   BannerView m_bannerView;

apps/sequence/graph/graph_controller.cpp

@@ -101,4 +101,52 @@ double GraphController::defaultCursorT(Ion::Storage::Record record) {
   return std::fmax(0.0, std::round(Shared::FunctionGraphController::defaultCursorT(record)));
 }
 
+InteractiveCurveViewRangeDelegate::Range GraphController::computeYRange(InteractiveCurveViewRange * interactiveCurveViewRange) {
+  Poincare::Context * context = textFieldDelegateApp()->localContext();
+  float min = FLT_MAX;
+  float max = -FLT_MAX;
+  float xMin = interactiveCurveViewRange->xMin();
+  float xMax = interactiveCurveViewRange->xMax();
+  assert(functionStore()->numberOfActiveFunctions() > 0);
+  for (int i = 0; i < functionStore()->numberOfActiveFunctions(); i++) {
+    ExpiringPointer<Shared::Function> f = functionStore()->modelForRecord(functionStore()->activeRecordAtIndex(i));
+    /* Scan x-range from the middle to the extrema in order to get balanced
+     * y-range for even functions (y = 1/x). */
+    double tMin = f->tMin();
+    if (std::isnan(tMin)) {
+      tMin = xMin;
+    } else if (f->shouldClipTRangeToXRange()) {
+      tMin = std::max<double>(tMin, xMin);
+    }
+    double tMax = f->tMax();
+    if (std::isnan(tMax)) {
+      tMax = xMax;
+    } else if (f->shouldClipTRangeToXRange()) {
+      tMax = std::min<double>(tMax, xMax);
+    }
+  /* In practice, a step smaller than a pixel's width is needed for sampling
+   * the values of a function. Otherwise some relevant extremal values may be
+   * missed. */
+    float rangeStep = f->rangeStep();
+    const float step = std::isnan(rangeStep) ? curveView()->pixelWidth() / 2.0f : rangeStep;
+    const int balancedBound = std::floor((tMax-tMin)/2/step);
+    for (int j = -balancedBound; j <= balancedBound ; j++) {
+      float t = (tMin+tMax)/2 + step * j;
+      Coordinate2D<float> xy = f->evaluateXYAtParameter(t, context);
+      float x = xy.x1();
+      if (!std::isnan(x) && !std::isinf(x) && x >= xMin && x <= xMax) {
+        float y = xy.x2();
+        if (!std::isnan(y) && !std::isinf(y)) {
+          min = std::min(min, y);
+          max = std::max(max, y);
+        }
+      }
+    }
+  }
+  InteractiveCurveViewRangeDelegate::Range range;
+  range.min = min;
+  range.max = max;
+  return range;
+}
+
 }

apps/sequence/graph/graph_controller.h

@@ -27,6 +27,7 @@ private:
   bool handleEnter() override;
   bool moveCursorHorizontally(int direction, int scrollSpeed = 1) override;
   double defaultCursorT(Ion::Storage::Record record) override;
+  InteractiveCurveViewRangeDelegate::Range computeYRange(Shared::InteractiveCurveViewRange * interactiveCurveViewRange) override;
   CurveViewRange * interactiveCurveViewRange() override { return m_graphRange; }
   SequenceStore * functionStore() const override { return static_cast<SequenceStore *>(Shared::FunctionGraphController::functionStore()); }
   GraphView * functionGraphView() override { return &m_view; }

apps/shared/continuous_function.cpp

@@ -261,6 +261,278 @@ void ContinuousFunction::setTMax(float tMax) {
   setCache(nullptr);
 }
 
+void ContinuousFunction::rangeForDisplay(float * xMin, float * xMax, float * yMin, float * yMax, Poincare::Context * context, bool tuneXRange) const {
+  if (plotType() == PlotType::Cartesian) {
+    interestingXAndYRangesForDisplay(xMin, xMax, yMin, yMax, context, tuneXRange);
+  } else {
+    fullXYRange(xMin, xMax, yMin, yMax, context);
+  }
+}
+
+void ContinuousFunction::fullXYRange(float * xMin, float * xMax, float * yMin, float * yMax, Context * context) const {
+  assert(yMin && yMax);
+  assert(!(std::isinf(tMin()) || std::isinf(tMax()) || std::isnan(rangeStep())));
+
+  float resultXMin = FLT_MAX, resultXMax = - FLT_MAX, resultYMin = FLT_MAX, resultYMax = - FLT_MAX;
+  for (float t = tMin(); t <= tMax(); t += rangeStep()) {
+    Coordinate2D<float> xy = privateEvaluateXYAtParameter(t, context);
+    if (!std::isfinite(xy.x1()) || !std::isfinite(xy.x2())) {
+      continue;
+    }
+    resultXMin = std::min(xy.x1(), resultXMin);
+    resultXMax = std::max(xy.x1(), resultXMax);
+    resultYMin = std::min(xy.x2(), resultYMin);
+    resultYMax = std::max(xy.x2(), resultYMax);
+  }
+  if (xMin) {
+    *xMin = resultXMin;
+  }
+  if (xMax) {
+    *xMax = resultXMax;
+  }
+  *yMin = resultYMin;
+  *yMax = resultYMax;
+}
+
+static float evaluateAndRound(const ContinuousFunction * f, float x, Context * context, float precision = 1e-5) {
+  /* When evaluating sin(x)/x close to zero using the standard sine function,
+   * one can detect small varitions, while the cardinal sine is supposed to be
+   * locally monotonous. To smooth our such variations, we round the result of
+   * the evaluations. As we are not interested in precise results but only in
+   * ordering, this approximation is sufficient. */
+  return precision * std::round(f->evaluateXYAtParameter(x, context).x2() / precision);
+}
+
+
+/* TODO : These three methods perform checks that will also be relevant for the
+ * equation solver. Remember to factorize this code when integrating the new
+ * solver. */
+static bool boundOfIntervalOfDefinitionIsReached(float y1, float y2) {
+  return std::isfinite(y1) && !std::isinf(y2) && std::isnan(y2);
+}
+static bool rootExistsOnInterval(float y1, float y2) {
+  return ((y1 < 0.f && y2 > 0.f) || (y1 > 0.f && y2 < 0.f));
+}
+static bool extremumExistsOnInterval(float y1, float y2, float y3) {
+  return (y1 < y2 && y2 > y3) || (y1 > y2 && y2 < y3);
+}
+
+/* This function checks whether an interval contains an extremum or an
+ * asymptote, by recursively computing the slopes. In case of an extremum, the
+ * slope should taper off toward the center. */
+static bool isExtremum(const ContinuousFunction * f, float x1, float x2, float x3, float y1, float y2, float y3, Context * context, int iterations = 3) {
+  if (iterations <= 0) {
+    return false;
+  }
+  float x[2] = {x1, x3}, y[2] = {y1, y3};
+  float xm, ym;
+  for (int i = 0; i < 2; i++) {
+    xm = (x[i] + x2) / 2.f;
+    ym = evaluateAndRound(f, xm, context);
+    bool res = ((y[i] < ym) != (ym < y2)) ? isExtremum(f, x[i], xm, x2, y[i], ym, y2, context, iterations - 1) : std::fabs(ym - y[i]) >= std::fabs(y2 - ym);
+    if (!res) {
+      return false;
+    }
+  }
+  return true;
+}
+
+enum class PointOfInterest : uint8_t {
+  None,
+  Bound,
+  Extremum,
+  Root
+};
+
+void ContinuousFunction::interestingXAndYRangesForDisplay(float * xMin, float * xMax, float * yMin, float * yMax, Context * context, bool tuneXRange) const {
+  assert(xMin && xMax && yMin && yMax);
+  assert(plotType() == PlotType::Cartesian);
+
+  /* Constants of the algorithm. */
+  constexpr float defaultMaxInterval = 2e5f;
+  constexpr float minDistance = 1e-2f;
+  constexpr float asymptoteThreshold = 2e-1f;
+  constexpr float stepFactor = 1.1f;
+  constexpr int maxNumberOfPoints = 3;
+  constexpr float breathingRoom = 0.3f;
+  constexpr float maxRatioBetweenPoints = 100.f;
+
+  const bool hasIntervalOfDefinition = std::isfinite(tMin()) && std::isfinite(tMax());
+  float center, maxDistance;
+  if (!tuneXRange) {
+    center = (*xMax + *xMin) / 2.f;
+    maxDistance = (*xMax - *xMin) / 2.f;
+  } else if (hasIntervalOfDefinition) {
+    center = (tMax() + tMin()) / 2.f;
+    maxDistance = (tMax() - tMin()) / 2.f;
+  } else {
+    center = 0.f;
+    maxDistance = defaultMaxInterval / 2.f;
+  }
+
+  float resultX[2] = {FLT_MAX, - FLT_MAX};
+  float resultYMin = FLT_MAX, resultYMax = - FLT_MAX;
+  float asymptote[2] = {FLT_MAX, - FLT_MAX};
+  int numberOfPoints;
+  float xFallback, yFallback[2] = {NAN, NAN};
+  float firstResult;
+  float dXOld, dXPrev, dXNext, yOld, yPrev, yNext;
+
+  /* Look for an point of interest at the center. */
+  const float a = center - minDistance - FLT_EPSILON, b = center + FLT_EPSILON, c = center + minDistance + FLT_EPSILON;
+  const float ya = evaluateAndRound(this, a, context), yb = evaluateAndRound(this, b, context), yc = evaluateAndRound(this, c, context);
+  if (boundOfIntervalOfDefinitionIsReached(ya, yc) ||
+      boundOfIntervalOfDefinitionIsReached(yc, ya) ||
+      rootExistsOnInterval(ya, yc) ||
+      extremumExistsOnInterval(ya, yb, yc) || ya == yc)
+  {
+    resultX[0] = resultX[1] = center;
+    if (extremumExistsOnInterval(ya, yb, yc) && isExtremum(this, a, b, c, ya, yb, yc, context)) {
+      resultYMin = resultYMax = yb;
+    }
+  }
+
+  /* We search for points of interest by exploring the function leftward from
+   * the center and then rightward, hence the two iterations. */
+  for (int i = 0; i < 2; i++) {
+    /* Initialize the search parameters. */
+    numberOfPoints = 0;
+    firstResult = NAN;
+    xFallback = NAN;
+    dXPrev = i == 0 ? - minDistance : minDistance;
+    dXNext = dXPrev * stepFactor;
+    yPrev = evaluateAndRound(this, center + dXPrev, context);
+    yNext = evaluateAndRound(this, center + dXNext, context);
+
+    while(std::fabs(dXPrev) < maxDistance) {
+      /* Update the slider. */
+      dXOld = dXPrev;
+      dXPrev = dXNext;
+      dXNext *= stepFactor;
+      yOld = yPrev;
+      yPrev = yNext;
+      yNext = evaluateAndRound(this, center + dXNext, context);
+      if (std::isinf(yNext)) {
+        continue;
+      }
+      /* Check for a change in the profile. */
+      const PointOfInterest variation = boundOfIntervalOfDefinitionIsReached(yPrev, yNext) ? PointOfInterest::Bound :
+                                        rootExistsOnInterval(yPrev, yNext) ? PointOfInterest::Root :
+                                        extremumExistsOnInterval(yOld, yPrev, yNext) ? PointOfInterest::Extremum :
+                                        PointOfInterest::None;
+      switch (static_cast<uint8_t>(variation)) {
+        /* The fall through is intentional, as we only want to update the Y
+         * range when an extremum is detected, but need to update the X range
+         * in all cases. */
+        case static_cast<uint8_t>(PointOfInterest::Extremum):
+        if (isExtremum(this, center + dXOld, center + dXPrev, center + dXNext, yOld, yPrev, yNext, context)) {
+          resultYMin = std::min(resultYMin, yPrev);
+          resultYMax = std::max(resultYMax, yPrev);
+        }
+      case static_cast<uint8_t>(PointOfInterest::Bound):
+        /* We only count extrema / discontinuities for limiting the number
+         * of points. This prevents cos(x) and cos(x)+2 from having different
+         * profiles. */
+        if (++numberOfPoints == maxNumberOfPoints) {
+          /* When too many points are encountered, we prepare their erasure by
+           * setting a restore point. */
+          xFallback = dXNext + center;
+          yFallback[0] = resultYMin;
+          yFallback[1] = resultYMax;
+        }
+      case static_cast<uint8_t>(PointOfInterest::Root):
+        asymptote[i] = i == 0 ? FLT_MAX : - FLT_MAX;
+        resultX[0] = std::min(resultX[0], center + (i == 0 ? dXNext : dXPrev));
+        resultX[1] = std::max(resultX[1], center + (i == 1 ? dXNext : dXPrev));
+        if (std::isnan(firstResult)) {
+          firstResult = dXNext;
+        }
+        break;
+      default:
+        const float slopeNext = (yNext - yPrev) / (dXNext - dXPrev), slopePrev = (yPrev - yOld) / (dXPrev - dXOld);
+        if ((std::fabs(slopeNext) < asymptoteThreshold) && (std::fabs(slopePrev) > asymptoteThreshold)) {
+          // Horizontal asymptote begins
+          asymptote[i] = (i == 0) ? std::min(asymptote[i], center + dXNext) : std::max(asymptote[i], center + dXNext);
+        } else if ((std::fabs(slopeNext) < asymptoteThreshold) && (std::fabs(slopePrev) > asymptoteThreshold)) {
+          // Horizontal asymptote invalidates : it might be an asymptote when
+          // going the other way.
+          asymptote[1 - i] = (i == 1) ? std::min(asymptote[1 - i], center + dXPrev) : std::max(asymptote[1 - i], center + dXPrev);
+        }
+      }
+    }
+    if (std::fabs(resultX[i]) > std::fabs(firstResult) * maxRatioBetweenPoints && !std::isnan(xFallback)) {
+      /* When there are too many points, cut them if their orders are too
+       * different. */
+      resultX[i] = xFallback;
+      resultYMin = yFallback[0];
+      resultYMax = yFallback[1];
+    }
+  }
+
+  if (tuneXRange) {
+    /* Cut after horizontal asymptotes. */
+    resultX[0] = std::min(resultX[0], asymptote[0]);
+    resultX[1] = std::max(resultX[1], asymptote[1]);
+    if (resultX[0] >= resultX[1]) {
+      /* Fallback to default range. */
+      resultX[0] = - Range1D::k_default;
+      resultX[1] = Range1D::k_default;
+    } else {
+      /* Add breathing room around points of interest. */
+      float xRange = resultX[1] - resultX[0];
+      resultX[0] -= breathingRoom * xRange;
+      resultX[1] += breathingRoom * xRange;
+      /* Round to the next integer. */
+      resultX[0] = std::floor(resultX[0]);
+      resultX[1] = std::ceil(resultX[1]);
+    }
+    *xMin = std::min(resultX[0], *xMin);
+    *xMax = std::max(resultX[1], *xMax);
+  }
+  *yMin = std::min(resultYMin, *yMin);
+  *yMax = std::max(resultYMax, *yMax);
+
+  refinedYRangeForDisplay(*xMin, *xMax, yMin, yMax, context);
+}
+
+void ContinuousFunction::refinedYRangeForDisplay(float xMin, float xMax, float * yMin, float * yMax, Context * context) const {
+  /* This methods computes the Y range that will be displayed for the cartesian
+   * function, given an X range (xMin, xMax) and bounds yMin and yMax that must
+   * be inside the Y range.*/
+  assert(plotType() == PlotType::Cartesian);
+  assert(yMin && yMax);
+
+  constexpr int sampleSize = 80;
+
+  float sampleYMin = FLT_MAX, sampleYMax = -FLT_MAX;
+  const float step = (xMax - xMin) / (sampleSize - 1);
+  float x, y;
+  float sum = 0.f;
+  int pop = 0;
+
+  for (int i = 1; i < sampleSize; i++) {
+    x = xMin + i * step;
+    y = privateEvaluateXYAtParameter(x, context).x2();
+    sampleYMin = std::min(sampleYMin, y);
+    sampleYMax = std::max(sampleYMax, y);
+    if (std::isfinite(y) && std::fabs(y) > FLT_EPSILON) {
+      sum += std::log(std::fabs(y));
+      pop++;
+    }
+  }
+  /* sum/n is the log mean value of the function, which can be interpreted as
+   * its average order of magnitude. Then, bound is the value for the next
+   * order of magnitude and is used to cut the Y range. */
+  float bound = (pop > 0) ? std::exp(sum / pop + 1.f) : FLT_MAX;
+  *yMin = std::min(*yMin, std::max(sampleYMin, -bound));
+  *yMax = std::max(*yMax, std::min(sampleYMax, bound));
+  if (*yMin == *yMax) {
+    float d = (*yMin == 0.f) ? 1.f : *yMin * 0.2f;
+    *yMin -= d;
+    *yMax += d;
+  }
+}
+
 void * ContinuousFunction::Model::expressionAddress(const Ion::Storage::Record * record) const {
   return (char *)record->value().buffer+sizeof(RecordDataBuffer);
 }

apps/shared/continuous_function.h

@@ -70,6 +70,8 @@ public:
   void setTMax(float tMax);
   float rangeStep() const override { return plotType() == PlotType::Cartesian ? NAN : (tMax() - tMin())/k_polarParamRangeSearchNumberOfPoints; }
 
+  void rangeForDisplay(float * xMin, float * xMax, float * yMin, float * yMax, Poincare::Context * context, bool tuneXRange = true) const;
+
   // Extremum
   Poincare::Coordinate2D<double> nextMinimumFrom(double start, double step, double max, Poincare::Context * context) const;
   Poincare::Coordinate2D<double> nextMaximumFrom(double start, double step, double max, Poincare::Context * context) const;
@@ -88,6 +90,12 @@ private:
   typedef Poincare::Coordinate2D<double> (*ComputePointOfInterest)(Poincare::Expression e, char * symbol, double start, double step, double max, Poincare::Context * context);
   Poincare::Coordinate2D<double> nextPointOfInterestFrom(double start, double step, double max, Poincare::Context * context, ComputePointOfInterest compute) const;
   template <typename T> Poincare::Coordinate2D<T> privateEvaluateXYAtParameter(T t, Poincare::Context * context) const;
+
+  // Ranges
+  void fullXYRange(float * xMin, float * xMax, float * yMin, float * yMax, Poincare::Context * context) const;
+  void interestingXAndYRangesForDisplay(float * xMin, float * xMax, float * yMin, float * yMax, Poincare::Context * context, bool tuneXRange = true) const;
+  void refinedYRangeForDisplay(float xMin, float xMax, float * yMin, float * yMax, Poincare::Context * context) const;
+
   /* RecordDataBuffer is the layout of the data buffer of Record
    * representing a ContinuousFunction. See comment on
    * Shared::Function::RecordDataBuffer about packing. */

apps/shared/function_graph_controller.cpp

@@ -73,53 +73,6 @@ void FunctionGraphController::reloadBannerView() {
   reloadBannerViewForCursorOnFunction(m_cursor, record, functionStore(), AppsContainer::sharedAppsContainer()->globalContext());
 }
 
-InteractiveCurveViewRangeDelegate::Range FunctionGraphController::computeYRange(InteractiveCurveViewRange * interactiveCurveViewRange) {
-  Poincare::Context * context = textFieldDelegateApp()->localContext();
-  float min = FLT_MAX;
-  float max = -FLT_MAX;
-  float xMin = interactiveCurveViewRange->xMin();
-  float xMax = interactiveCurveViewRange->xMax();
-  assert(functionStore()->numberOfActiveFunctions() > 0);
-  for (int i = 0; i < functionStore()->numberOfActiveFunctions(); i++) {
-    ExpiringPointer<Function> f = functionStore()->modelForRecord(functionStore()->activeRecordAtIndex(i));
-    /* Scan x-range from the middle to the extrema in order to get balanced
-     * y-range for even functions (y = 1/x). */
-    double tMin = f->tMin();
-    if (std::isnan(tMin)) {
-      tMin = xMin;
-    } else if (f->shouldClipTRangeToXRange()) {
-      tMin = std::max<double>(tMin, xMin);
-    }
-    double tMax = f->tMax();
-    if (std::isnan(tMax)) {
-      tMax = xMax;
-    } else if (f->shouldClipTRangeToXRange()) {
-      tMax = std::min<double>(tMax, xMax);
-    }
-  /* In practice, a step smaller than a pixel's width is needed for sampling
-   * the values of a function. Otherwise some relevant extremal values may be
-   * missed. */
-    float rangeStep = f->rangeStep();
-    const float step = std::isnan(rangeStep) ? curveView()->pixelWidth() / 2.0f : rangeStep;
-    const int balancedBound = std::floor((tMax-tMin)/2/step);
-    for (int j = -balancedBound; j <= balancedBound ; j++) {
-      float t = (tMin+tMax)/2 + step * j;
-      Coordinate2D<float> xy = f->evaluateXYAtParameter(t, context);
-      float x = xy.x1();
-      if (!std::isnan(x) && !std::isinf(x) && x >= xMin && x <= xMax) {
-        float y = xy.x2();
-        if (!std::isnan(y) && !std::isinf(y)) {
-          min = std::min(min, y);
-          max = std::max(max, y);
-        }
-      }
-    }
-  }
-  InteractiveCurveViewRangeDelegate::Range range;
-  range.min = min;
-  range.max = max;
-  return range;
-}
 
 double FunctionGraphController::defaultCursorT(Ion::Storage::Record record) {
   return (interactiveCurveViewRange()->xMin()+interactiveCurveViewRange()->xMax())/2.0f;

apps/shared/function_graph_controller.h

@@ -44,9 +44,6 @@ private:
   virtual FunctionGraphView * functionGraphView() = 0;
   virtual FunctionCurveParameterController * curveParameterController() = 0;
 
-  // InteractiveCurveViewRangeDelegate
-  InteractiveCurveViewRangeDelegate::Range computeYRange(InteractiveCurveViewRange * interactiveCurveViewRange) override;
-
   // InteractiveCurveViewController
   bool moveCursorVertically(int direction) override;
   uint32_t modelVersion() override;

apps/shared/interactive_curve_view_range.cpp

@@ -90,72 +90,11 @@ void InteractiveCurveViewRange::roundAbscissa() {
 void InteractiveCurveViewRange::normalize() {
   /* We center the ranges on the current range center, and put each axis so that
    * 1cm = 2 current units. */
-  m_yAuto = false;
-
-  const float unit = std::max(xGridUnit(), yGridUnit());
-
-  // Set x range
-  float newXHalfRange = NormalizedXHalfRange(unit);
-  float newXMin = xCenter() - newXHalfRange;
-  float newXMax = xCenter() + newXHalfRange;
-  if (!std::isnan(newXMin) && !std::isnan(newXMax)) {
-    m_xRange.setMax(newXMax, k_lowerMaxFloat, k_upperMaxFloat);
-    MemoizedCurveViewRange::protectedSetXMin(newXMin, k_lowerMaxFloat, k_upperMaxFloat);
-  }
-  // Set y range
-  float newYHalfRange = NormalizedYHalfRange(unit);
-  float newYMin = yCenter() - newYHalfRange;
-  float newYMax = yCenter() + newYHalfRange;
-  if (!std::isnan(newYMin) && !std::isnan(newYMax)) {
-    m_yRange.setMax(newYMax, k_lowerMaxFloat, k_upperMaxFloat);
-    MemoizedCurveViewRange::protectedSetYMin(newYMin, k_lowerMaxFloat, k_upperMaxFloat);
-  }
-}
-
-void InteractiveCurveViewRange::setTrigonometric() {
-  m_yAuto = false;
-  // Set x range
-  float x = (Preferences::sharedPreferences()->angleUnit() == Preferences::AngleUnit::Degree) ? 600.0f : 10.5f;
-  m_xRange.setMax(x, k_lowerMaxFloat, k_upperMaxFloat);
-  MemoizedCurveViewRange::protectedSetXMin(-x, k_lowerMaxFloat, k_upperMaxFloat);
-  // Set y range
-  float y = 1.6f;
-  m_yRange.setMax(y, k_lowerMaxFloat, k_upperMaxFloat);
-  MemoizedCurveViewRange::protectedSetYMin(-y, k_lowerMaxFloat, k_upperMaxFloat);
-}
-
-void InteractiveCurveViewRange::setDefault() {
-  if (m_delegate == nullptr) {
-    return;
-  }
-  if (!m_delegate->defautRangeIsNormalized()) {
-    m_yAuto = true;
-    m_xRange.setMax(m_delegate->interestingXHalfRange(), k_lowerMaxFloat, k_upperMaxFloat);
-    setXMin(-xMax());
-    return;
-  }
 
   m_yAuto = false;
-  // Compute the interesting ranges
-  float a,b,c,d;
-  m_delegate->interestingRanges(&a, &b, &c, &d);
-  m_xRange.setMin(a, k_lowerMaxFloat, k_upperMaxFloat);
-  m_xRange.setMax(b, k_lowerMaxFloat, k_upperMaxFloat);
-  m_yRange.setMin(c, k_lowerMaxFloat, k_upperMaxFloat);
-  m_yRange.setMax(d, k_lowerMaxFloat, k_upperMaxFloat);
 
-  // Add margins
   float xRange = xMax() - xMin();
   float yRange = yMax() - yMin();
-  m_xRange.setMin(m_delegate->addMargin(xMin(), xRange, false, true), k_lowerMaxFloat, k_upperMaxFloat);
-  // Use MemoizedCurveViewRange::protectedSetXMax to update xGridUnit
-  MemoizedCurveViewRange::protectedSetXMax(m_delegate->addMargin(xMax(), xRange, false, false), k_lowerMaxFloat, k_upperMaxFloat);
-  m_yRange.setMin(m_delegate->addMargin(yMin(), yRange, true, true), k_lowerMaxFloat, k_upperMaxFloat);
-  MemoizedCurveViewRange::protectedSetYMax(m_delegate->addMargin(yMax(), yRange, true, false), k_lowerMaxFloat, k_upperMaxFloat);
-
-  // Normalize the axes, so that a polar circle is displayed as a circle
-  xRange = xMax() - xMin();
-  yRange = yMax() - yMin();
   float xyRatio = xRange/yRange;
 
   const float unit = std::max(xGridUnit(), yGridUnit());
@@ -177,6 +116,49 @@ void InteractiveCurveViewRange::setDefault() {
   }
 }
 
+void InteractiveCurveViewRange::setTrigonometric() {
+  m_yAuto = false;
+  // Set x range
+  float x = (Preferences::sharedPreferences()->angleUnit() == Preferences::AngleUnit::Degree) ? 600.0f : 10.5f;
+  m_xRange.setMax(x, k_lowerMaxFloat, k_upperMaxFloat);
+  MemoizedCurveViewRange::protectedSetXMin(-x, k_lowerMaxFloat, k_upperMaxFloat);
+  // Set y range
+  float y = 1.6f;
+  m_yRange.setMax(y, k_lowerMaxFloat, k_upperMaxFloat);
+  MemoizedCurveViewRange::protectedSetYMin(-y, k_lowerMaxFloat, k_upperMaxFloat);
+}
+
+void InteractiveCurveViewRange::setDefault() {
+  if (m_delegate == nullptr) {
+    return;
+  }
+
+  // Compute the interesting range
+  m_yAuto = false;
+  float xm, xM, ym, yM;
+  m_delegate->interestingRanges(&xm, &xM, &ym, &yM);
+  m_xRange.setMin(xm, k_lowerMaxFloat, k_upperMaxFloat);
+  m_xRange.setMax(xM, k_lowerMaxFloat, k_upperMaxFloat);
+  m_yRange.setMin(ym, k_lowerMaxFloat, k_upperMaxFloat);
+  m_yRange.setMax(yM, k_lowerMaxFloat, k_upperMaxFloat);
+
+  // Add margins
+  float xRange = xMax() - xMin();
+  float yRange = yMax() - yMin();
+  m_xRange.setMin(m_delegate->addMargin(xMin(), xRange, false, true), k_lowerMaxFloat, k_upperMaxFloat);
+  // Use MemoizedCurveViewRange::protectedSetXMax to update xGridUnit
+  MemoizedCurveViewRange::protectedSetXMax(m_delegate->addMargin(xMax(), xRange, false, false), k_lowerMaxFloat, k_upperMaxFloat);
+  m_yRange.setMin(m_delegate->addMargin(yMin(), yRange, true, true), k_lowerMaxFloat, k_upperMaxFloat);
+  MemoizedCurveViewRange::protectedSetYMax(m_delegate->addMargin(yMax(), yRange, true, false), k_lowerMaxFloat, k_upperMaxFloat);
+
+  if (!m_delegate->defaultRangeIsNormalized()) {
+    return;
+  }
+
+  // Normalize the axes, so that a polar circle is displayed as a circle
+  normalize();
+}
+
 void InteractiveCurveViewRange::centerAxisAround(Axis axis, float position) {
   if (std::isnan(position)) {
     return;

apps/shared/interactive_curve_view_range_delegate.h

@@ -12,7 +12,7 @@ public:
   bool didChangeRange(InteractiveCurveViewRange * interactiveCurveViewRange);
   virtual float interestingXMin() const { return -interestingXHalfRange(); }
   virtual float interestingXHalfRange() const { return 10.0f; }
-  virtual bool defautRangeIsNormalized() const { return false; }
+  virtual bool defaultRangeIsNormalized() const { return false; }
   virtual void interestingRanges(float * xm, float * xM, float * ym, float * yM) const { assert(false); }
   virtual float addMargin(float x, float range, bool isVertical, bool isMin) = 0;
 protected: