diff --git a/apps/probability/distribution/binomial_distribution.cpp b/apps/probability/distribution/binomial_distribution.cpp
index 76f3e22f4..8a9b386a3 100644
--- a/apps/probability/distribution/binomial_distribution.cpp
+++ b/apps/probability/distribution/binomial_distribution.cpp
@@ -1,4 +1,5 @@
 #include "binomial_distribution.h"
+#include <poincare/binomial_distribution.h>
 #include <assert.h>
 #include <cmath>
 
@@ -22,6 +23,10 @@ I18n::Message BinomialDistribution::parameterDefinitionAtIndex(int index) {
   }
 }
 
+float BinomialDistribution::evaluateAtAbscissa(float x) const {
+  return Poincare::BinomialDistribution::EvaluateAtAbscissa<float>(x, m_parameter1, m_parameter2);
+}
+
 float BinomialDistribution::xMin() const {
   float min = 0.0f;
   float max = m_parameter1 > 0.0f ? m_parameter1 : 1.0f;
@@ -68,13 +73,7 @@ bool BinomialDistribution::authorizedValueAtIndex(float x, int index) const {
 }
 
 double BinomialDistribution::cumulativeDistributiveInverseForProbability(double * probability) {
-  if (m_parameter1 == 0.0 && (m_parameter2 == 0.0 || m_parameter2 == 1.0)) {
-    return NAN;
-  }
-  if (*probability >= 1.0) {
-    return m_parameter1;
-  }
-  return Distribution::cumulativeDistributiveInverseForProbability(probability);
+  return Poincare::BinomialDistribution::CumulativeDistributiveInverseForProbability<double>(*probability, m_parameter1, m_parameter2);
 }
 
 double BinomialDistribution::rightIntegralInverseForProbability(double * probability) {
@@ -87,38 +86,8 @@ double BinomialDistribution::rightIntegralInverseForProbability(double * probabi
   return Distribution::rightIntegralInverseForProbability(probability);
 }
 
-template<typename T>
-T BinomialDistribution::templatedApproximateAtAbscissa(T x) const {
-  if (m_parameter1 == 0) {
-    if (m_parameter2 == 0 || m_parameter2 == 1) {
-      return NAN;
-    }
-    if (floor(x) == 0) {
-      return 1;
-    }
-    return 0;
-  }
-  if (m_parameter2 == 1) {
-    if (floor(x) == m_parameter1) {
-      return 1;
-    }
-    return 0;
-  }
-  if (m_parameter2 == 0) {
-    if (floor(x) == 0) {
-      return 1;
-    }
-    return 0;
-  }
-  if (x > m_parameter1) {
-    return 0;
-  }
-  T lResult = std::lgamma((T)(m_parameter1+1.0)) - std::lgamma(std::floor(x)+(T)1.0) - std::lgamma((T)m_parameter1 - std::floor(x)+(T)1.0)+
-    std::floor(x)*std::log((T)m_parameter2) + ((T)m_parameter1-std::floor(x))*std::log((T)(1.0-m_parameter2));
-  return std::exp(lResult);
+double BinomialDistribution::evaluateAtDiscreteAbscissa(int k) const {
+  return Poincare::BinomialDistribution::EvaluateAtAbscissa<double>((double) k, (double)m_parameter1, (double)m_parameter2);
 }
 
 }
-
-template float Probability::BinomialDistribution::templatedApproximateAtAbscissa(float x) const;
-template double Probability::BinomialDistribution::templatedApproximateAtAbscissa(double x) const;
diff --git a/apps/probability/distribution/binomial_distribution.h b/apps/probability/distribution/binomial_distribution.h
index 7b8811d4e..c55711fa7 100644
--- a/apps/probability/distribution/binomial_distribution.h
+++ b/apps/probability/distribution/binomial_distribution.h
@@ -16,17 +16,12 @@ public:
   float yMax() const override;
   I18n::Message parameterNameAtIndex(int index) override;
   I18n::Message parameterDefinitionAtIndex(int index) override;
-  float evaluateAtAbscissa(float x) const override {
-    return templatedApproximateAtAbscissa(x);
-  }
+  float evaluateAtAbscissa(float x) const override;
   bool authorizedValueAtIndex(float x, int index) const override;
   double cumulativeDistributiveInverseForProbability(double * probability) override;
   double rightIntegralInverseForProbability(double * probability) override;
 protected:
-  double evaluateAtDiscreteAbscissa(int k) const override {
-    return templatedApproximateAtAbscissa((double)k);
-  }
-  template<typename T> T templatedApproximateAtAbscissa(T x) const;
+  double evaluateAtDiscreteAbscissa(int k) const override;
 };
 
 }
diff --git a/apps/probability/distribution/distribution.cpp b/apps/probability/distribution/distribution.cpp
index 789a6a88d..f31ec66e4 100644
--- a/apps/probability/distribution/distribution.cpp
+++ b/apps/probability/distribution/distribution.cpp
@@ -7,21 +7,12 @@ namespace Probability {
 
 double Distribution::cumulativeDistributiveFunctionAtAbscissa(double x) const {
   if (!isContinuous()) {
-    int end = std::round(x);
-    double result = 0.0;
-    for (int k = 0; k <=end; k++) {
-      result += evaluateAtDiscreteAbscissa(k);
-      /* Avoid too long loop */
-      if (k > k_maxNumberOfOperations) {
-        break;
-      }
-      if (result >= k_maxProbability) {
-        result = 1.0;
-        break;
-      }
-
-    }
-    return result;
+    return Poincare::Solver::CumulativeDistributiveFunctionForNDefinedFunction<double>(x,
+        [](double k, Poincare::Context * context, Poincare::Preferences::ComplexFormat complexFormat, Poincare::Preferences::AngleUnit angleUnit, const void * context1, const void * context2, const void * context3) {
+        const Distribution * distribution = reinterpret_cast<const Distribution *>(context1);
+        return distribution->evaluateAtDiscreteAbscissa(k);
+      }, nullptr, Poincare::Preferences::ComplexFormat::Real, Poincare::Preferences::AngleUnit::Degree, this);
+    // Context, complex format and angle unit are dummy values
   }
   return 0.0;
 }
diff --git a/apps/probability/distribution/normal_distribution.cpp b/apps/probability/distribution/normal_distribution.cpp
index 691b7d0bd..932d97254 100644
--- a/apps/probability/distribution/normal_distribution.cpp
+++ b/apps/probability/distribution/normal_distribution.cpp
@@ -1,9 +1,7 @@
 #include "normal_distribution.h"
 #include <poincare/normal_distribution.h>
-#include <assert.h>
 #include <cmath>
 #include <float.h>
-#include <ion.h>
 
 namespace Probability {
 
diff --git a/poincare/Makefile b/poincare/Makefile
index 7e103af34..0e8756511 100644
--- a/poincare/Makefile
+++ b/poincare/Makefile
@@ -30,6 +30,7 @@ poincare_src += $(addprefix poincare/src/,\
 
 poincare_src += $(addprefix poincare/src/,\
   init.cpp \
+  binomial_distribution.cpp \
   erf_inv.cpp \
   exception_checkpoint.cpp \
   helpers.cpp \
@@ -44,7 +45,9 @@ poincare_src += $(addprefix poincare/src/,\
   arc_sine.cpp \
   arc_tangent.cpp \
   arithmetic.cpp \
+  binom_cdf.cpp \
   binomial_coefficient.cpp \
+  binomial_distribution_function.cpp \
   ceiling.cpp \
   complex.cpp \
   complex_argument.cpp \
diff --git a/poincare/include/poincare/binom_cdf.h b/poincare/include/poincare/binom_cdf.h
new file mode 100644
index 000000000..48f491256
--- /dev/null
+++ b/poincare/include/poincare/binom_cdf.h
@@ -0,0 +1,46 @@
+#ifndef POINCARE_BINOM_CDF_H
+#define POINCARE_BINOM_CDF_H
+
+#include <poincare/approximation_helper.h>
+#include <poincare/binomial_distribution_function.h>
+
+namespace Poincare {
+
+class BinomCDFNode final : public BinomialDistributionFunctionNode  {
+public:
+
+  // TreeNode
+  size_t size() const override { return sizeof(BinomCDFNode); }
+  int numberOfChildren() const override;
+#if POINCARE_TREE_LOG
+  virtual void logNodeName(std::ostream & stream) const override {
+    stream << "BinomCDF";
+  }
+#endif
+
+  // Properties
+  Type type() const override { return Type::BinomCDF; }
+  Sign sign(Context * context) const override { return Sign::Positive; }
+  Expression setSign(Sign s, ReductionContext reductionContext) override;
+
+private:
+  // Layout
+  Layout createLayout(Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const override;
+  int serialize(char * buffer, int bufferSize, Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const override;
+
+  // Evaluation
+  Evaluation<float> approximate(SinglePrecision p, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const override { return templatedApproximate<float>(context, complexFormat, angleUnit); }
+  Evaluation<double> approximate(DoublePrecision p, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const override { return templatedApproximate<double>(context, complexFormat, angleUnit); }
+  template<typename T> Evaluation<T> templatedApproximate(Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const;
+};
+
+class BinomCDF final : public BinomialDistributionFunction {
+public:
+  BinomCDF(const BinomCDFNode * n) : BinomialDistributionFunction(n) {}
+  static BinomCDF Builder(Expression child0, Expression child1, Expression child2) { return TreeHandle::FixedArityBuilder<BinomCDF, BinomCDFNode>(ArrayBuilder<TreeHandle>(child0, child1, child2).array(), 3); }
+  static constexpr Expression::FunctionHelper s_functionHelper = Expression::FunctionHelper("binomcdf", 3, &UntypedBuilderThreeChildren<BinomCDF>);
+};
+
+}
+
+#endif
diff --git a/poincare/include/poincare/binomial_distribution.h b/poincare/include/poincare/binomial_distribution.h
new file mode 100644
index 000000000..0169fbb84
--- /dev/null
+++ b/poincare/include/poincare/binomial_distribution.h
@@ -0,0 +1,22 @@
+#ifndef POINCARE_BINOMIAL_DISTRIBUTION_H
+#define POINCARE_BINOMIAL_DISTRIBUTION_H
+
+#include <poincare/expression.h>
+#include <poincare/preferences.h>
+
+namespace Poincare {
+
+class BinomialDistribution final {
+public:
+  template<typename T> static T EvaluateAtAbscissa(T x, T n, T p);
+  template<typename T> static T CumulativeDistributiveFunctionAtAbscissa(T x, T n, T p);
+  template<typename T> static T CumulativeDistributiveInverseForProbability(T probability, T n, T p);
+  template<typename T> static bool ParametersAreOK(T n, T p);
+  /* ExpressionParametersAreOK returns true if the expression could be verified.
+   * The result of the verification is *result. */
+  static bool ExpressionParametersAreOK(bool * result, const Expression & n, const Expression & p, Context * context);
+};
+
+}
+
+#endif
diff --git a/poincare/include/poincare/binomial_distribution_function.h b/poincare/include/poincare/binomial_distribution_function.h
new file mode 100644
index 000000000..c99c1d5c4
--- /dev/null
+++ b/poincare/include/poincare/binomial_distribution_function.h
@@ -0,0 +1,26 @@
+#ifndef POINCARE_BINOMIAL_DISTRIBUTION_FUNCTION_H
+#define POINCARE_BINOMIAL_DISTRIBUTION_FUNCTION_H
+
+#include <poincare/expression.h>
+
+namespace Poincare {
+
+// BinomialDistributionFunctions are BinomCDF, InvBinom and BinomPDF
+
+class BinomialDistributionFunctionNode : public ExpressionNode {
+private:
+  // Simplication
+  Expression shallowReduce(ReductionContext reductionContext) override;
+  LayoutShape leftLayoutShape() const override { return LayoutShape::MoreLetters; };
+  LayoutShape rightLayoutShape() const override { return LayoutShape::BoundaryPunctuation; }
+};
+
+class BinomialDistributionFunction : public Expression {
+public:
+  BinomialDistributionFunction(const BinomialDistributionFunctionNode * n) : Expression(n) {}
+  Expression shallowReduce(Context * context, bool * stopReduction = nullptr);
+};
+
+}
+
+#endif
diff --git a/poincare/include/poincare/expression.h b/poincare/include/poincare/expression.h
index 10ada5808..604fc5c98 100644
--- a/poincare/include/poincare/expression.h
+++ b/poincare/include/poincare/expression.h
@@ -27,6 +27,7 @@ class Expression : public TreeHandle {
   friend class ArcTangent;
   friend class Arithmetic;
   friend class BinomialCoefficient;
+  friend class BinomialDistributionFunction;
   friend class Ceiling;
   friend class CommonLogarithm;
   template<typename T>
diff --git a/poincare/include/poincare/expression_node.h b/poincare/include/poincare/expression_node.h
index 75355f414..b579969d4 100644
--- a/poincare/include/poincare/expression_node.h
+++ b/poincare/include/poincare/expression_node.h
@@ -47,6 +47,7 @@ public:
     ArcCosine,
     ArcSine,
     ArcTangent,
+    BinomCDF,
     BinomialCoefficient,
     Ceiling,
     ComplexArgument,
diff --git a/poincare/include/poincare/solver.h b/poincare/include/poincare/solver.h
index 7fde84e69..414a6ea9b 100644
--- a/poincare/include/poincare/solver.h
+++ b/poincare/include/poincare/solver.h
@@ -9,11 +9,25 @@ namespace Poincare {
 
 class Solver {
 public:
+  // Minimum
   typedef double (*ValueAtAbscissa)(double abscissa, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, const void * context1, const void * context2, const void * context3);
   static Coordinate2D BrentMinimum(double ax, double bx, ValueAtAbscissa evaluation, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, const void * context1 = nullptr, const void * context2 = nullptr, const void * context3 = nullptr);
+
+  // Root
   static double BrentRoot(double ax, double bx, double precision, ValueAtAbscissa evaluation, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, const void * context1 = nullptr, const void * context2 = nullptr, const void * context3 = nullptr);
   static Coordinate2D IncreasingFunctionRoot(double ax, double bx, double precision, ValueAtAbscissa evaluation, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, const void * context1 = nullptr, const void * context2 = nullptr, const void * context3 = nullptr);
+
+  // Proba
+
+  // Cumulative distributive inverse for function defined on N (positive integers)
+  template<typename T> static T CumulativeDistributiveInverseForNDefinedFunction(T probability, ValueAtAbscissa evaluation, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, const void * context1 = nullptr, const void * context2 = nullptr, const void * context3 = nullptr);
+
+  // Cumulative distributive function for function defined on N (positive integers)
+  template<typename T> static T CumulativeDistributiveFunctionForNDefinedFunction(T x, ValueAtAbscissa evaluation, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, const void * context1 = nullptr, const void * context2 = nullptr, const void * context3 = nullptr);
+
 private:
+  constexpr static int k_maxNumberOfOperations = 1000000;
+  constexpr static double k_maxProbability = 0.9999995;
   constexpr static double k_sqrtEps = 1.4901161193847656E-8; // sqrt(DBL_EPSILON)
   constexpr static double k_goldenRatio = 0.381966011250105151795413165634361882279690820194237137864; // (3-sqrt(5))/2
 };
diff --git a/poincare/include/poincare_nodes.h b/poincare/include/poincare_nodes.h
index 308baf921..b5b19ea2a 100644
--- a/poincare/include/poincare_nodes.h
+++ b/poincare/include/poincare_nodes.h
@@ -6,6 +6,7 @@
 #include <poincare/arc_cosine.h>
 #include <poincare/arc_sine.h>
 #include <poincare/arc_tangent.h>
+#include <poincare/binom_cdf.h>
 #include <poincare/binomial_coefficient.h>
 #include <poincare/complex_argument.h>
 #include <poincare/confidence_interval.h>
diff --git a/poincare/src/binom_cdf.cpp b/poincare/src/binom_cdf.cpp
new file mode 100644
index 000000000..43e8d87f9
--- /dev/null
+++ b/poincare/src/binom_cdf.cpp
@@ -0,0 +1,40 @@
+#include <poincare/binom_cdf.h>
+#include <poincare/layout_helper.h>
+#include <poincare/binomial_distribution.h>
+#include <poincare/serialization_helper.h>
+#include <assert.h>
+
+namespace Poincare {
+
+constexpr Expression::FunctionHelper BinomCDF::s_functionHelper;
+
+int BinomCDFNode::numberOfChildren() const { return BinomCDF::s_functionHelper.numberOfChildren(); }
+
+Expression BinomCDFNode::setSign(Sign s, ReductionContext reductionContext) {
+  assert(s == Sign::Positive);
+  return BinomCDF(this);
+}
+
+Layout BinomCDFNode::createLayout(Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const {
+  return LayoutHelper::Prefix(BinomCDF(this), floatDisplayMode, numberOfSignificantDigits, BinomCDF::s_functionHelper.name());
+}
+
+int BinomCDFNode::serialize(char * buffer, int bufferSize, Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const {
+  return SerializationHelper::Prefix(this, buffer, bufferSize, floatDisplayMode, numberOfSignificantDigits, BinomCDF::s_functionHelper.name());
+}
+
+template<typename T>
+Evaluation<T> BinomCDFNode::templatedApproximate(Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const {
+  Evaluation<T> xEvaluation = childAtIndex(0)->approximate(T(), context, complexFormat, angleUnit);
+  Evaluation<T> nEvaluation = childAtIndex(1)->approximate(T(), context, complexFormat, angleUnit);
+  Evaluation<T> pEvaluation = childAtIndex(2)->approximate(T(), context, complexFormat, angleUnit);
+
+  const T x = xEvaluation.toScalar();
+  const T n = nEvaluation.toScalar();
+  const T p = pEvaluation.toScalar();
+
+  // CumulativeDistributiveFunctionAtAbscissa handles bad mu and var values
+  return Complex<T>::Builder(BinomialDistribution::CumulativeDistributiveFunctionAtAbscissa(x, n, p));
+}
+
+}
diff --git a/poincare/src/binomial_distribution.cpp b/poincare/src/binomial_distribution.cpp
new file mode 100644
index 000000000..660958beb
--- /dev/null
+++ b/poincare/src/binomial_distribution.cpp
@@ -0,0 +1,175 @@
+#include <poincare/binomial_distribution.h>
+#include <poincare/rational.h>
+#include <cmath>
+#include <float.h>
+#include <assert.h>
+
+namespace Poincare {
+
+template<typename T>
+T BinomialDistribution::EvaluateAtAbscissa(T x, T n, T p) {
+  if (std::isnan(x) || std::isinf(x) || !ParametersAreOK(n, p)){
+    return NAN;
+  }
+  T precision = sizeof(T) == sizeof(double) ? DBL_EPSILON : FLT_EPSILON;
+  bool nIsZero = std::abs(n) < precision;
+  bool pIsZero = std::abs(p) < precision;
+  bool pIsOne = !pIsZero && std::abs(p - (T)1.0) < precision;
+  if (nIsZero) {
+    if (pIsZero || pIsOne) {
+      return NAN;
+    }
+    if (std::floor(x) == 0) {
+      return (T)1;
+    }
+    return (T)0;
+  }
+  if (pIsOne) {
+    return (T)(floor(x) == n ? 1 : 0);
+  }
+  if (pIsZero) {
+    return (T)(floor(x) == 0.0 ? 1 : 0);
+  }
+  if (x > n) {
+    return(T)0;
+  }
+  T lResult = std::lgamma(n+(T)1.0) - std::lgamma(std::floor(x)+(T)1.0) - std::lgamma(n - std::floor(x)+(T)1.0)+
+    std::floor(x)*std::log(p) + (n-std::floor(x))*std::log((T)(1.0)-p);
+  return std::exp(lResult);
+}
+
+template<typename T>
+T BinomialDistribution::CumulativeDistributiveFunctionAtAbscissa(T x, T n, T p) {
+  if (!ParametersAreOK(n, p) || std::isnan(x) || std::isinf(x)) {
+    return NAN;
+  }
+  if (x < (T)0.0) {
+    return (T)0.0;
+  }
+  if (x >= n) {
+    return (T)1.0;
+  }
+  bool isDouble = sizeof(T) == sizeof(double);
+  return Solver::CumulativeDistributiveFunctionForNDefinedFunction<T>(
+      x,
+      [](double abscissa, Context * context, Poincare::Preferences::ComplexFormat complexFormat, Poincare::Preferences::AngleUnit angleUnit, const void * n, const void * p, const void * isDouble) {
+        if (*(bool *)isDouble) {
+          return (double)BinomialDistribution::EvaluateAtAbscissa<T>(abscissa, *(reinterpret_cast<const double *>(n)), *(reinterpret_cast<const double *>(p)));
+        }
+        return (double)BinomialDistribution::EvaluateAtAbscissa<T>(abscissa, *(reinterpret_cast<const float *>(n)), *(reinterpret_cast<const float *>(p)));
+      }, (Context *)nullptr, Preferences::ComplexFormat::Real, Preferences::AngleUnit::Degree, &n, &p, &isDouble);
+    // Context, complex format and angle unit are dummy values
+}
+
+template<typename T>
+T BinomialDistribution::CumulativeDistributiveInverseForProbability(T probability, T n, T p) {
+  if (!ParametersAreOK(n, p) || std::isnan(probability) || std::isinf(probability) || probability < (T)0.0 || probability > (T)1.0) {
+    return NAN;
+  }
+  if (!ParametersAreOK(n, p)) {
+    return NAN;
+  }
+  bool isDouble = sizeof(T) == sizeof(double);
+  T precision = isDouble ? DBL_EPSILON : FLT_EPSILON;
+  bool nIsZero = std::abs(n) < precision;
+  bool pIsZero = std::abs(p) < precision;
+  bool pIsOne = !pIsZero && std::abs(p - (T)1.0) < precision;
+  if (nIsZero && (pIsZero || pIsOne)) {
+    return NAN;
+  }
+
+  if (std::abs(probability - (T)1.0) < precision) {
+    return n;
+  }
+
+  return Solver::CumulativeDistributiveInverseForNDefinedFunction<T>(
+      probability,
+      [](double x, Context * context, Poincare::Preferences::ComplexFormat complexFormat, Poincare::Preferences::AngleUnit angleUnit, const void * n, const void * p, const void * isDouble) {
+        if (*(bool *)isDouble) {
+          return (double)BinomialDistribution::EvaluateAtAbscissa<T>(x, *(reinterpret_cast<const double *>(n)), *(reinterpret_cast<const double *>(p)));
+        }
+        return (double)BinomialDistribution::EvaluateAtAbscissa<T>(x, *(reinterpret_cast<const float *>(n)), *(reinterpret_cast<const float *>(p)));
+      }, (Context *)nullptr, Preferences::ComplexFormat::Real, Preferences::AngleUnit::Degree, &n, &p, &isDouble);
+    // Context, complex format and angle unit are dummy values
+}
+
+template<typename T>
+bool BinomialDistribution::ParametersAreOK(T n, T p) {
+  /* TODO As the cumulative probability are computed by looping over all discrete
+   * abscissa within the interesting range, the complexity of the cumulative
+   * probability is linear with the size of the range. Here we cap the maximal
+   * size of the range to 10000. If one day we want to increase or get rid of
+   *  this cap, we should implement the explicit formula of the cumulative
+   *  probability (which depends on an incomplete beta function) to make the
+   *  comlexity O(1). */
+  //TODO LEA n doit être <= 99999.0f) {
+  return !std::isnan(n) && !std::isnan(p)
+    && !std::isinf(n) && !std::isinf(p)
+    && (n == (int)n) && (n >= (T)0)
+    && (p >= (T)0) && (p <= (T)1);
+}
+
+bool BinomialDistribution::ExpressionParametersAreOK(bool * result, const Expression & n, const Expression & p, Context * context) {
+  assert(result != nullptr);
+  if (n.deepIsMatrix(context) || p.deepIsMatrix(context)) {
+    *result = false;
+    return true;
+  }
+
+  if (n.isUndefined() || p.isUndefined() || Expression::IsInfinity(n, context) || Expression::IsInfinity(p,context)) {
+    *result = false;
+    return true;
+  }
+  if (!n.isReal(context) || !p.isReal(context)) {
+    // We cannot check that n and p are real
+    return false;
+  }
+
+  if (n.type() != ExpressionNode::Type::Rational) {
+    // We cannot check that that n is an integer
+    return false;
+  }
+
+  {
+    const Rational rationalN = static_cast<const Rational &>(n);
+    if (rationalN.isNegative() || !rationalN.isInteger()) {
+      // n is negative or not an integer
+      *result = false;
+      return true;
+    }
+  }
+
+  if (p.type() != ExpressionNode::Type::Rational) {
+    // We cannot check that that p is >= 0 and <= 1
+    return false;
+  }
+
+  const Rational rationalP = static_cast<const Rational &>(p);
+  if (rationalP.isNegative()) {
+    // p is negative
+    *result = false;
+    return true;
+  }
+  Integer num = rationalP.unsignedIntegerNumerator();
+  Integer den = rationalP.integerDenominator();
+  if (den.isLowerThan(num)) {
+    // p > 1
+    *result = false;
+    return true;
+  }
+
+  *result = true;
+  return true;
+}
+
+
+template float BinomialDistribution::EvaluateAtAbscissa<float>(float, float, float);
+template double BinomialDistribution::EvaluateAtAbscissa<double>(double, double, double);
+template float BinomialDistribution::CumulativeDistributiveFunctionAtAbscissa<float>(float, float, float);
+template double BinomialDistribution::CumulativeDistributiveFunctionAtAbscissa<double>(double, double, double);
+template float BinomialDistribution::CumulativeDistributiveInverseForProbability<float>(float, float, float);
+template double BinomialDistribution::CumulativeDistributiveInverseForProbability<double>(double, double, double);
+template bool BinomialDistribution::ParametersAreOK(float, float);
+template bool BinomialDistribution::ParametersAreOK(double, double);
+
+}
diff --git a/poincare/src/binomial_distribution_function.cpp b/poincare/src/binomial_distribution_function.cpp
new file mode 100644
index 000000000..2dcd2b96e
--- /dev/null
+++ b/poincare/src/binomial_distribution_function.cpp
@@ -0,0 +1,41 @@
+#include <poincare/binomial_distribution_function.h>
+#include <poincare/binomial_distribution.h>
+#include <assert.h>
+
+namespace Poincare {
+
+Expression BinomialDistributionFunctionNode::shallowReduce(ReductionContext reductionContext) {
+  return BinomialDistributionFunction(this).shallowReduce(reductionContext.context());
+}
+
+Expression BinomialDistributionFunction::shallowReduce(Context * context, bool * stopReduction) {
+  if (stopReduction != nullptr) {
+    *stopReduction = true;
+  }
+  {
+    Expression e = Expression::defaultShallowReduce();
+    if (e.isUndefined()) {
+      return e;
+    }
+  }
+
+  Expression n = childAtIndex(1);
+  Expression p = childAtIndex(2);
+
+  // Check mu and var
+  bool nAndPOK = false;
+  bool couldCheckNAndP = BinomialDistribution::ExpressionParametersAreOK(&nAndPOK, n, p, context);
+  if (!couldCheckNAndP) {
+    return *this;
+  }
+  if (!nAndPOK) {
+    return replaceWithUndefinedInPlace();
+  }
+
+  if (stopReduction != nullptr) {
+    *stopReduction = false;
+  }
+  return *this;
+}
+
+}
diff --git a/poincare/src/parsing/parser.h b/poincare/src/parsing/parser.h
index 3a3cdba30..cae31c9fd 100644
--- a/poincare/src/parsing/parser.h
+++ b/poincare/src/parsing/parser.h
@@ -99,6 +99,7 @@ private:
     &HyperbolicArcSine::s_functionHelper,
     &ArcTangent::s_functionHelper,
     &HyperbolicArcTangent::s_functionHelper,
+    &BinomCDF::s_functionHelper,
     &BinomialCoefficient::s_functionHelper,
     &Ceiling::s_functionHelper,
     &ConfidenceInterval::s_functionHelper,
diff --git a/poincare/src/solver.cpp b/poincare/src/solver.cpp
index aa3f76619..d9dd78f7c 100644
--- a/poincare/src/solver.cpp
+++ b/poincare/src/solver.cpp
@@ -205,4 +205,51 @@ Coordinate2D Solver::IncreasingFunctionRoot(double ax, double bx, double precisi
   return Coordinate2D(NAN, NAN);
 }
 
+template<typename T>
+T Solver::CumulativeDistributiveInverseForNDefinedFunction(T probability, ValueAtAbscissa evaluation, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, const void * context1, const void * context2, const void * context3) {
+  T precision = sizeof(T) == sizeof(double) ? DBL_EPSILON : FLT_EPSILON;
+  assert(probability <= (((T)1.0) - precision) && probability > precision);
+  T p = 0.0;
+  int k = 0;
+  T delta = 0.0;
+  do {
+    delta = std::fabs(probability-p);
+    p += evaluation(k++, context, complexFormat, angleUnit, context1, context2, context3);
+    if (p >= k_maxProbability && std::fabs(probability-1.0) <= delta) {
+      return (T)(k-1);
+    }
+  } while (std::fabs(probability-p) <= delta && k < k_maxNumberOfOperations && p < 1.0);
+  p -= evaluation(--k, context, complexFormat, angleUnit, context1, context2, context3);
+  if (k == k_maxNumberOfOperations) {
+    return INFINITY;
+  }
+  if (std::isnan(p)) {
+    return NAN;
+  }
+  return k-1;
+}
+
+template<typename T>
+T Solver::CumulativeDistributiveFunctionForNDefinedFunction(T x, ValueAtAbscissa evaluation, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, const void * context1, const void * context2, const void * context3) {
+  int end = std::round(x);
+  double result = 0.0;
+  for (int k = 0; k <=end; k++) {
+    result += evaluation(k, context, complexFormat, angleUnit, context1, context2, context3);
+    /* Avoid too long loop */
+    if (k > k_maxNumberOfOperations) {
+      break;
+    }
+    if (result >= k_maxProbability) {
+      result = 1.0;
+      break;
+    }
+  }
+  return result;
+}
+
+template float Solver::CumulativeDistributiveInverseForNDefinedFunction(float, ValueAtAbscissa, Context *, Preferences::ComplexFormat, Preferences::AngleUnit, const void *, const void *, const void *);
+template double Solver::CumulativeDistributiveInverseForNDefinedFunction(double, ValueAtAbscissa, Context *, Preferences::ComplexFormat, Preferences::AngleUnit, const void *, const void *, const void *);
+template float Solver::CumulativeDistributiveFunctionForNDefinedFunction(float, ValueAtAbscissa, Context *, Preferences::ComplexFormat, Preferences::AngleUnit, const void *, const void *, const void *);
+template double Solver::CumulativeDistributiveFunctionForNDefinedFunction(double, ValueAtAbscissa, Context *, Preferences::ComplexFormat, Preferences::AngleUnit, const void *, const void *, const void *);
+
 }
diff --git a/poincare/src/tree_handle.cpp b/poincare/src/tree_handle.cpp
index 185ddd231..1a06be7a3 100644
--- a/poincare/src/tree_handle.cpp
+++ b/poincare/src/tree_handle.cpp
@@ -266,6 +266,7 @@ template Addition TreeHandle::NAryBuilder<Addition, AdditionNode>(TreeHandle*, s
 template ArcCosine TreeHandle::FixedArityBuilder<ArcCosine, ArcCosineNode>(TreeHandle*, size_t);
 template ArcSine TreeHandle::FixedArityBuilder<ArcSine, ArcSineNode>(TreeHandle*, size_t);
 template ArcTangent TreeHandle::FixedArityBuilder<ArcTangent, ArcTangentNode>(TreeHandle*, size_t);
+template BinomCDF TreeHandle::FixedArityBuilder<BinomCDF, BinomCDFNode>(TreeHandle*, size_t);
 template BinomialCoefficient TreeHandle::FixedArityBuilder<BinomialCoefficient, BinomialCoefficientNode>(TreeHandle*, size_t);
 template BinomialCoefficientLayout TreeHandle::FixedArityBuilder<BinomialCoefficientLayout, BinomialCoefficientLayoutNode>(TreeHandle*, size_t);
 template Ceiling TreeHandle::FixedArityBuilder<Ceiling, CeilingNode>(TreeHandle*, size_t);
diff --git a/poincare/test/approximation.cpp b/poincare/test/approximation.cpp
index 66d56da94..7431010db 100644
--- a/poincare/test/approximation.cpp
+++ b/poincare/test/approximation.cpp
@@ -234,9 +234,15 @@ QUIZ_CASE(poincare_approximation_function) {
   assert_expression_approximates_to<float>("abs([[3+2𝐢,3+4𝐢][5+2𝐢,3+2𝐢]])", "[[3.605551,5][5.385165,3.605551]]");
   assert_expression_approximates_to<double>("abs([[3+2𝐢,3+4𝐢][5+2𝐢,3+2𝐢]])", "[[3.605551275464,5][5.3851648071345,3.605551275464]]");
 
+  assert_expression_approximates_to<float>("binomcdf(5.3, 9, 0.7)", "0.270341", Degree, Cartesian, 6); // FIXME: precision problem
+  assert_expression_approximates_to<double>("binomcdf(5.3, 9, 0.7)", "0.270340902");
+
   assert_expression_approximates_to<float>("binomial(10, 4)", "210");
   assert_expression_approximates_to<double>("binomial(10, 4)", "210");
 
+  assert_expression_approximates_to<float>("binompdf(4, 9, 0.7)", "0.0735138");
+  assert_expression_approximates_to<double>("binompdf(5.3, 9, 0.7)", "0735138");
+
   assert_expression_approximates_to<float>("ceil(0.2)", "1");
   assert_expression_approximates_to<double>("ceil(0.2)", "1");
 
@@ -270,6 +276,9 @@ QUIZ_CASE(poincare_approximation_function) {
   assert_expression_approximates_to<float>("int(x,x, 1, 2)", "1.5");
   assert_expression_approximates_to<double>("int(x,x, 1, 2)", "1.5");
 
+  assert_expression_approximates_to<float>("invbinom(0.9647324002, 15, 0.7)", "13");
+  assert_expression_approximates_to<double>("invbinom(0.9647324002, 15, 0.7)", "13");
+
   assert_expression_approximates_to<float>("invnorm(0.56, 1.3, 2.4)", "1.662326");
   //assert_expression_approximates_to<double>("invnorm(0.56, 1.3, 2.4)", "1.6623258450088"); FIXME precision error