diff --git a/apps/probability/calculation/left_integral_calculation.cpp b/apps/probability/calculation/left_integral_calculation.cpp
index d1de80e16..eec5e2b4e 100644
--- a/apps/probability/calculation/left_integral_calculation.cpp
+++ b/apps/probability/calculation/left_integral_calculation.cpp
@@ -52,6 +52,9 @@ void LeftIntegralCalculation::compute(int indexKnownElement) {
       return;
     }
     m_upperBound = m_distribution->cumulativeDistributiveInverseForProbability(&m_result);
+    if (std::isnan(m_upperBound)) {
+      m_result = NAN;
+    }
   }
 }
 
diff --git a/apps/probability/calculation/right_integral_calculation.cpp b/apps/probability/calculation/right_integral_calculation.cpp
index 04f79e7ff..a1d760392 100644
--- a/apps/probability/calculation/right_integral_calculation.cpp
+++ b/apps/probability/calculation/right_integral_calculation.cpp
@@ -52,6 +52,9 @@ void RightIntegralCalculation::compute(int indexKnownElement) {
       return;
     }
     m_lowerBound = m_distribution->rightIntegralInverseForProbability(&m_result);
+    if (std::isnan(m_lowerBound)) {
+      m_result = NAN;
+    }
   }
 }
 
diff --git a/apps/probability/distribution/distribution.cpp b/apps/probability/distribution/distribution.cpp
index 3949b974b..777519a99 100644
--- a/apps/probability/distribution/distribution.cpp
+++ b/apps/probability/distribution/distribution.cpp
@@ -150,7 +150,7 @@ double Distribution::cumulativeDistributiveInverseForProbabilityUsingIncreasingF
       this,
       probability,
       nullptr);
-  assert(std::fabs(result.value()) < FLT_EPSILON*100); // TODO FLT_EPSILON is too strict
+  assert(std::isnan(result.value()) || std::fabs(result.value()) < FLT_EPSILON*100); // TODO FLT_EPSILON is too strict
   return result.abscissa();
 }
 
diff --git a/poincare/src/solver.cpp b/poincare/src/solver.cpp
index 94c3e1172..aa3f76619 100644
--- a/poincare/src/solver.cpp
+++ b/poincare/src/solver.cpp
@@ -186,7 +186,7 @@ Coordinate2D Solver::IncreasingFunctionRoot(double ax, double bx, double precisi
       return Coordinate2D(currentAbscissa, eval);
     }
     // The minimal value is already bigger than 0, return NAN.
-    return Coordinate2D(currentAbscissa, NAN);
+    return Coordinate2D(NAN, NAN);
   }
   while (max - min > precision) {
     currentAbscissa = (min + max) / 2.0;
@@ -196,10 +196,13 @@ Coordinate2D Solver::IncreasingFunctionRoot(double ax, double bx, double precisi
     } else if (eval < -DBL_EPSILON) {
       min = currentAbscissa;
     } else {
-      return Coordinate2D(currentAbscissa, eval);
+      break;
     }
   }
-  return Coordinate2D(currentAbscissa, std::fabs(eval) < precision ? eval : NAN);
+  if (std::fabs(eval) < precision) {
+    return Coordinate2D(currentAbscissa, eval);
+  }
+  return Coordinate2D(NAN, NAN);
 }
 
 }