diff --git a/apps/shared/continuous_function.cpp b/apps/shared/continuous_function.cpp index 0d316a5e1..36d6b4390 100644 --- a/apps/shared/continuous_function.cpp +++ b/apps/shared/continuous_function.cpp @@ -263,7 +263,7 @@ void ContinuousFunction::setTMax(float tMax) { void ContinuousFunction::rangeForDisplay(float * xMin, float * xMax, float * yMin, float * yMax, Poincare::Context * context, bool tuneXRange) const { if (plotType() == PlotType::Cartesian) { - Function::rangeForDisplay(xMin, xMax, yMin, yMax, context, tuneXRange); + protectedRangeForDisplay(xMin, xMax, yMin, yMax, context, tuneXRange, true); } else { fullXYRange(xMin, xMax, yMin, yMax, context); } diff --git a/apps/shared/continuous_function.h b/apps/shared/continuous_function.h index 7d1580c01..c66be00f7 100644 --- a/apps/shared/continuous_function.h +++ b/apps/shared/continuous_function.h @@ -92,7 +92,6 @@ private: template Poincare::Coordinate2D privateEvaluateXYAtParameter(T t, Poincare::Context * context) const; void fullXYRange(float * xMin, float * xMax, float * yMin, float * yMax, Poincare::Context * context) const; - void refinedYRangeForDisplay(float xMin, float xMax, float * yMin, float * yMax, Poincare::Context * context) const override { protectedRefinedYRangeForDisplay(xMin, xMax, yMin, yMax, context, true); } /* RecordDataBuffer is the layout of the data buffer of Record * representing a ContinuousFunction. See comment on diff --git a/apps/shared/function.cpp b/apps/shared/function.cpp index 9491c3956..6cfc69963 100644 --- a/apps/shared/function.cpp +++ b/apps/shared/function.cpp @@ -2,6 +2,7 @@ #include "range_1D.h" #include "poincare_helpers.h" #include "poincare/src/parsing/parser.h" +#include #include #include #include @@ -78,257 +79,21 @@ Function::RecordDataBuffer * Function::recordData() const { return reinterpret_cast(const_cast(d.buffer)); } -// Ranges -static float evaluateAndRound(const Function * f, float x, Context * context, float precision = 1e-5) { +void Function::protectedRangeForDisplay(float * xMin, float * xMax, float * yMin, float * yMax, Poincare::Context * context, bool tuneXRange, bool boundByMagnitude) const { + Zoom::ValueAtAbscissa evaluation = [](float x, Context * context, const void * auxiliary) { /* 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 + * one can detect small variations, 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); + constexpr float precision = 1e-5; + return precision * std::round(static_cast(auxiliary)->evaluateXYAtParameter(x, context).x2() / precision); + }; + Zoom::InterestingRangesForDisplay(evaluation, xMin, xMax, yMin, yMax, tMin(), tMax(), context, this, tuneXRange); + + evaluation = [](float x, Context * context, const void * auxiliary) { + return static_cast(auxiliary)->evaluateXYAtParameter(x, context).x2(); + }; + Zoom::RefinedYRangeForDisplay(evaluation, *xMin, *xMax, yMin, yMax, context, this, boundByMagnitude); } - -/* 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 Function * 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 Function::rangeForDisplay(float * xMin, float * xMax, float * yMin, float * yMax, Context * context, bool tuneXRange) const { - assert(xMin && xMax && yMin && yMax); - - /* 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 a 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(variation)) { - /* The fallthrough 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(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(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(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 Function::protectedRefinedYRangeForDisplay(float xMin, float xMax, float * yMin, float * yMax, Context * context, bool boundByMagnitude) const { - /* This methods computes the Y range that will be displayed for cartesian - * functions and sequences, given an X range (xMin, xMax) and bounds yMin and - * yMax that must be inside the Y range.*/ - assert(yMin && yMax); - - constexpr int sampleSize = Ion::Display::Width / 4; - constexpr float boundNegligigbleThreshold = 0.2f; - - 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 = evaluateXYAtParameter(x, context).x2(); - if (!std::isfinite(y)) { - continue; - } - sampleYMin = std::min(sampleYMin, y); - sampleYMax = std::max(sampleYMax, y); - if (std::fabs(y) > FLT_EPSILON) { - sum += std::log(std::fabs(y)); - pop++; - } - } - /* sum/pop 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. */ - if (boundByMagnitude) { - float bound = (pop > 0) ? std::exp(sum / pop + 1.f) : FLT_MAX; - sampleYMin = std::max(sampleYMin, - bound); - sampleYMax = std::min(sampleYMax, bound); - } - *yMin = std::min(*yMin, sampleYMin); - *yMax = std::max(*yMax, sampleYMax); - if (*yMin == *yMax) { - float d = (*yMin == 0.f) ? 1.f : *yMin * 0.2f; - *yMin -= d; - *yMax += d; - } - /* Round out the smallest bound to 0 if it is negligible compare to the - * other one. This way, we can display the X axis for positive functions such - * as sqrt(x) even if we do not sample close to 0. */ - if (*yMin > 0.f && *yMin / *yMax < boundNegligigbleThreshold) { - *yMin = 0.f; - } else if (*yMax < 0.f && *yMax / *yMin < boundNegligigbleThreshold) { - *yMax = 0.f; - } -} - } diff --git a/apps/shared/function.h b/apps/shared/function.h index 453940025..49442d1b3 100644 --- a/apps/shared/function.h +++ b/apps/shared/function.h @@ -55,7 +55,7 @@ public: virtual Poincare::Expression sumBetweenBounds(double start, double end, Poincare::Context * context) const = 0; // Range - virtual void rangeForDisplay(float * xMin, float * xMax, float * yMin, float * yMax, Poincare::Context * context, bool tuneXRange = true) const; + virtual void rangeForDisplay(float * xMin, float * xMax, float * yMin, float * yMax, Poincare::Context * context, bool tuneXRange = true) const = 0; protected: /* RecordDataBuffer is the layout of the data buffer of Record @@ -93,11 +93,9 @@ protected: bool m_active; }; - void protectedRefinedYRangeForDisplay(float xMin, float xMax, float * yMin, float * yMax, Poincare::Context * context, bool boundByMagnitude) const; + void protectedRangeForDisplay(float * xMin, float * xMax, float * yMin, float * yMax, Poincare::Context * context, bool tuneXRange, bool boundByMagnitude) const; private: - virtual void refinedYRangeForDisplay(float xMin, float xMax, float * yMin, float * yMax, Poincare::Context * context) const = 0; - RecordDataBuffer * recordData() const; }; diff --git a/apps/shared/range_1D.h b/apps/shared/range_1D.h index 58ef2537f..b89b316d8 100644 --- a/apps/shared/range_1D.h +++ b/apps/shared/range_1D.h @@ -3,6 +3,7 @@ #include #include +#include #if __EMSCRIPTEN__ #include @@ -19,7 +20,7 @@ public: /* If m_min and m_max are too close, we cannot divide properly the range by * the number of pixels, which creates a drawing problem. */ constexpr static float k_minFloat = 1E-4f; - constexpr static float k_default = 10.0f; + constexpr static float k_default = Poincare::Zoom::k_defaultHalfRange; Range1D(float min = -k_default, float max = k_default) : m_min(min), m_max(max) diff --git a/apps/shared/sequence.h b/apps/shared/sequence.h index e3c6b92ad..58ac6e7ad 100644 --- a/apps/shared/sequence.h +++ b/apps/shared/sequence.h @@ -74,6 +74,10 @@ public: Poincare::Expression sumBetweenBounds(double start, double end, Poincare::Context * context) const override; constexpr static int k_initialRankNumberOfDigits = 3; // m_initialRank is capped by 999 + + //Range + void rangeForDisplay(float * xMin, float * xMax, float * yMin, float * yMax, Poincare::Context * context, bool tuneXRange = true) const override { protectedRangeForDisplay(xMin, xMax, yMin, yMax, context, tuneXRange, false); }; + private: constexpr static const KDFont * k_layoutFont = KDFont::LargeFont; @@ -153,7 +157,6 @@ private: }; template T templatedApproximateAtAbscissa(T x, SequenceContext * sqctx) const; - void refinedYRangeForDisplay(float xMin, float xMax, float * yMin, float * yMax, Poincare::Context * context) const override { protectedRefinedYRangeForDisplay(xMin, xMax, yMin, yMax, context, false); } size_t metaDataSize() const override { return sizeof(RecordDataBuffer); } const Shared::ExpressionModel * model() const override { return &m_definition; } RecordDataBuffer * recordData() const; diff --git a/poincare/Makefile b/poincare/Makefile index 790c7fa6c..61ee32fe3 100644 --- a/poincare/Makefile +++ b/poincare/Makefile @@ -156,6 +156,7 @@ poincare_src += $(addprefix poincare/src/,\ vector_cross.cpp \ vector_dot.cpp \ vector_norm.cpp \ + zoom.cpp \ ) poincare_src += $(addprefix poincare/src/parsing/,\ diff --git a/poincare/include/poincare/zoom.h b/poincare/include/poincare/zoom.h new file mode 100644 index 000000000..69a33d5f7 --- /dev/null +++ b/poincare/include/poincare/zoom.h @@ -0,0 +1,50 @@ +#ifndef POINCARE_ZOOM_H +#define POINCARE_ZOOM_H + +#include +#include + +namespace Poincare { + +class Zoom { +public: + static constexpr float k_defaultHalfRange = 10.f; + + typedef float (*ValueAtAbscissa)(float abscissa, Context * context, const void * auxiliary); + + static void InterestingRangesForDisplay(ValueAtAbscissa evaluation, float * xMin, float * xMax, float * yMin, float * yMax, float tMin, float tMax, Context * context, const void * auxiliary, bool tuneXRange = false); + static void RefinedYRangeForDisplay(ValueAtAbscissa evaluation, float xMin, float xMax, float * yMin, float * yMax, Context * context, const void * auxiliary, bool boundByMagnitude = false); + +private: + static constexpr int k_peakNumberOfPointsOfInterest = 3; + static constexpr int k_sampleSize = Ion::Display::Width / 4; + static constexpr float k_defaultMaxInterval = 2e5f; + static constexpr float k_minimalDistance = 1e-2f; + static constexpr float k_asymptoteThreshold = 2e-1f; + static constexpr float k_stepFactor = 1.1f; + static constexpr float k_breathingRoom = 0.3f; + static constexpr float k_forceXAxisThreshold = 0.2f; + static constexpr float k_maxRatioBetweenPointsOfInterest = 100.f; + + enum class PointOfInterest : uint8_t { + None, + Bound, + Extremum, + Root + }; + + /* TODO : These 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); } + /* IsConvexAroundExtremum 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 IsConvexAroundExtremum(ValueAtAbscissa evaluation, float x1, float x2, float x3, float y1, float y2, float y3, Context * context, const void * auxiliary, int iterations = 3); +}; + +} + +#endif diff --git a/poincare/src/zoom.cpp b/poincare/src/zoom.cpp new file mode 100644 index 000000000..315f255bf --- /dev/null +++ b/poincare/src/zoom.cpp @@ -0,0 +1,226 @@ +#include +#include +#include + +namespace Poincare { + +constexpr int + Zoom::k_peakNumberOfPointsOfInterest, + Zoom::k_sampleSize; +constexpr float + Zoom::k_defaultMaxInterval, + Zoom::k_minimalDistance, + Zoom::k_asymptoteThreshold, + Zoom::k_stepFactor, + Zoom::k_breathingRoom, + Zoom::k_forceXAxisThreshold, + Zoom::k_defaultHalfRange, + Zoom::k_maxRatioBetweenPointsOfInterest; + +void Zoom::InterestingRangesForDisplay(ValueAtAbscissa evaluation, float * xMin, float * xMax, float * yMin, float * yMax, float tMin, float tMax, Context * context, const void * auxiliary, bool tuneXRange) { + assert(xMin && xMax && yMin && yMax); + + 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 = k_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 a point of interest at the center. */ + const float a = center - k_minimalDistance - FLT_EPSILON, b = center + FLT_EPSILON, c = center + k_minimalDistance + FLT_EPSILON; + const float ya = evaluation(a, context, auxiliary), yb = evaluation(b, context, auxiliary), yc = evaluation(c, context, auxiliary); + 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) && IsConvexAroundExtremum(evaluation, a, b, c, ya, yb, yc, context, auxiliary)) { + 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 ? - k_minimalDistance : k_minimalDistance; + dXNext = dXPrev * k_stepFactor; + yPrev = evaluation(center + dXPrev, context, auxiliary); + yNext = evaluation(center + dXNext, context, auxiliary); + + while(std::fabs(dXPrev) < maxDistance) { + /* Update the slider. */ + dXOld = dXPrev; + dXPrev = dXNext; + dXNext *= k_stepFactor; + yOld = yPrev; + yPrev = yNext; + yNext = evaluation(center + dXNext, context, auxiliary); + 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(variation)) { + /* The fallthrough 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(PointOfInterest::Extremum): + if (IsConvexAroundExtremum(evaluation, center + dXOld, center + dXPrev, center + dXNext, yOld, yPrev, yNext, context, auxiliary)) { + resultYMin = std::min(resultYMin, yPrev); + resultYMax = std::max(resultYMax, yPrev); + } + case static_cast(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 == k_peakNumberOfPointsOfInterest) { + /* 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(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) < k_asymptoteThreshold) && (std::fabs(slopePrev) > k_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) < k_asymptoteThreshold) && (std::fabs(slopePrev) > k_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) * k_maxRatioBetweenPointsOfInterest && !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] = - k_defaultHalfRange; + resultX[1] = k_defaultHalfRange; + } else { + /* Add breathing room around points of interest. */ + float xRange = resultX[1] - resultX[0]; + resultX[0] -= k_breathingRoom * xRange; + resultX[1] += k_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); +} + +void Zoom::RefinedYRangeForDisplay(ValueAtAbscissa evaluation, float xMin, float xMax, float * yMin, float * yMax, Context * context, const void * auxiliary, bool boundByMagnitude) { + /* This methods computes the Y range that will be displayed for cartesian + * functions and sequences, given an X range (xMin, xMax) and bounds yMin and + * yMax that must be inside the Y range.*/ + assert(yMin && yMax); + + float sampleYMin = FLT_MAX, sampleYMax = -FLT_MAX; + const float step = (xMax - xMin) / (k_sampleSize - 1); + float x, y; + float sum = 0.f; + int pop = 0; + + for (int i = 1; i < k_sampleSize; i++) { + x = xMin + i * step; + y = evaluation(x, context, auxiliary); + if (!std::isfinite(y)) { + continue; + } + sampleYMin = std::min(sampleYMin, y); + sampleYMax = std::max(sampleYMax, y); + if (std::fabs(y) > FLT_EPSILON) { + sum += std::log(std::fabs(y)); + pop++; + } + } + /* sum/pop 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. */ + if (boundByMagnitude) { + float bound = (pop > 0) ? std::exp(sum / pop + 1.f) : FLT_MAX; + sampleYMin = std::max(sampleYMin, - bound); + sampleYMax = std::min(sampleYMax, bound); + } + *yMin = std::min(*yMin, sampleYMin); + *yMax = std::max(*yMax, sampleYMax); + if (*yMin == *yMax) { + float d = (*yMin == 0.f) ? 1.f : *yMin * 0.2f; + *yMin -= d; + *yMax += d; + } + /* Round out the smallest bound to 0 if it is negligible compare to the + * other one. This way, we can display the X axis for positive functions such + * as sqrt(x) even if we do not sample close to 0. */ + if (*yMin > 0.f && *yMin / *yMax < k_forceXAxisThreshold) { + *yMin = 0.f; + } else if (*yMax < 0.f && *yMax / *yMin < k_forceXAxisThreshold) { + *yMax = 0.f; + } +} + +bool Zoom::IsConvexAroundExtremum(ValueAtAbscissa evaluation, float x1, float x2, float x3, float y1, float y2, float y3, Context * context, const void * auxiliary, int iterations) { + 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 = evaluation(xm, context, auxiliary); + bool res = ((y[i] < ym) != (ym < y2)) ? IsConvexAroundExtremum(evaluation, x[i], xm, x2, y[i], ym, y2, context, auxiliary, iterations - 1) : std::fabs(ym - y[i]) >= std::fabs(y2 - ym); + if (!res) { + return false; + } + } + return true; +} + +}