[apps/probability] Helper for contonued fractions and infinite series

apps/probability/Makefile

@@ -4,6 +4,8 @@ app_headers += apps/probability/app.h
 app_probability_test_src = $(addprefix apps/probability/,\
   law/erf_inv.cpp \
   law/law.cpp \
+  law/regularized_gamma.cpp \
+  law/helper.cpp \
 )
 
 app_probability_src = $(addprefix apps/probability/,\
@@ -45,9 +47,6 @@ i18n_files += $(addprefix apps/probability/,\
   base.universal.i18n\
 )
 
-tests_src += $(addprefix apps/probability/law/,\
-  regularized_gamma.cpp \
-)
 tests_src += $(addprefix apps/probability/test/,\
   erf_inv.cpp\
   regularized_gamma.cpp\

apps/probability/law/chi_squared_law.cpp

@@ -1,6 +1,5 @@
 #include "chi_squared_law.h"
 #include "regularized_gamma.h"
-#include <poincare/solver.h>
 #include <cmath>
 
 namespace Probability {

apps/probability/law/helper.cpp

@@ -0,0 +1,79 @@
+#include "helper.h"
+#include "law.h"
+#include <cmath>
+#include <float.h>
+#include <assert.h>
+
+void increaseIfTooSmall(double * d) {
+  assert(!std::isinf(*d) && !std::isnan(*d));
+  const double small = 1e-50;
+  if (std::fabs(*d) < small) {
+    *d = small;
+  }
+}
+
+bool Helper::ContinuedFractionEvaluation(ValueForIndex a, ValueForIndex b, int maxNumberOfIterations, double * result, double context1, double context2) {
+  // Lentz's method with Thompson and Barnett's modification
+  double hnminus1 = b(0, context1, context2);
+  increaseIfTooSmall(&hnminus1);
+  double cnminus1 = hnminus1;
+  double dnminus1 = 0.0;
+  double hn = hnminus1;
+
+  int iterationCount = 1;
+  while (iterationCount < maxNumberOfIterations) {
+    // Compute a(n) and b(n)
+    const double an = a(iterationCount, context1, context2);
+    const double bn = b(iterationCount, context1, context2);
+
+    // Compute c(n) and d(n)
+    double cn = bn + an / cnminus1;
+    double dn = bn + an * dnminus1;
+    increaseIfTooSmall(&cn);
+    increaseIfTooSmall(&dn);
+    dn = 1 / dn;
+
+    // Compute delta(n) and h(n)
+    const double deltaN = cn * dn;
+    hn = hnminus1 * deltaN;
+
+    if (std::isinf(hn) || std::isnan(hn)) {
+      return false;
+    }
+    if (std::fabs(deltaN - 1.0) < 10e-15) {
+      break;
+    }
+
+    // Update c, d, h
+    dnminus1 = dn;
+    cnminus1 = cn;
+    hnminus1 = hn;
+    iterationCount++;
+  }
+
+  if (iterationCount >= maxNumberOfIterations) {
+    return false;
+  }
+
+  *result = hn;
+  return true;
+}
+
+bool Helper::InfiniteSeriesEvaluation(double firstTerm, TermUpdater termUpdater, double termLimit, int maxNumberOfIterations, double * result, double context1, double context2) {
+  double iterationCount = 0.0;
+  double currentTerm = firstTerm;
+  double sum = currentTerm;
+  while (iterationCount < maxNumberOfIterations
+      && !std::isinf(sum)
+      && std::fabs(currentTerm/sum) > termLimit)
+  {
+    iterationCount+= 1.0;
+    currentTerm = termUpdater(currentTerm, iterationCount, context1, context2);
+    sum+= currentTerm;
+  }
+  if (iterationCount >= maxNumberOfIterations) {
+    return false;
+  }
+  *result = sum;
+  return true;
+}

apps/probability/law/helper.h

@@ -0,0 +1,17 @@
+#ifndef PROBABILITY_LAW_HELPER_H
+#define PROBABILITY_LAW_HELPER_H
+
+class Helper {
+public:
+  // Continued fractions
+  typedef double (*ValueForIndex)(double index, double s, double x);
+  static bool ContinuedFractionEvaluation(ValueForIndex a, ValueForIndex b, int maxNumberOfIterations, double * result, double context1, double context2);
+
+  // Infinite series
+  typedef double (*TermUpdater)(double previousTerm, double index, double s, double x);
+  static bool InfiniteSeriesEvaluation(double firstTerm, TermUpdater termUpdater, double termLimit, int maxNumberOfIterations, double * result, double context1, double context2);
+
+};
+
+#endif
+

apps/probability/law/regularized_gamma.cpp

@@ -1,87 +1,10 @@
 #include "regularized_gamma.h"
 #include "law.h"
+#include "helper.h"
 #include <cmath>
 #include <float.h>
 #include <assert.h>
 
-void increaseIfTooSmall(double * d) {
-  assert(!std::isinf(*d) && !std::isnan(*d));
-  const double small = 1e-50;
-  if (std::fabs(*d) < small) {
-    *d = small;
-  }
-}
-
-typedef double (*ValueForIndex)(double index, double s, double x);
-
-static inline bool ContinuedFractionEvaluation(ValueForIndex a, ValueForIndex b, int maxNumberOfIterations, double * result, double context1, double context2) {
-  // Lentz's method with Thompson and Barnett's modification
-  double hnminus1 = b(0, context1, context2);
-  increaseIfTooSmall(&hnminus1);
-  double cnminus1 = hnminus1;
-  double dnminus1 = 0.0;
-  double hn = hnminus1;
-
-  int iterationCount = 1;
-  while (iterationCount < maxNumberOfIterations) {
-    // Compute a(n) and b(n)
-    const double an = a(iterationCount, context1, context2);
-    const double bn = b(iterationCount, context1, context2);
-
-    // Compute c(n) and d(n)
-    double cn = bn + an / cnminus1;
-    double dn = bn + an * dnminus1;
-    increaseIfTooSmall(&cn);
-    increaseIfTooSmall(&dn);
-    dn = 1 / dn;
-
-    // Compute delta(n) and h(n)
-    const double deltaN = cn * dn;
-    hn = hnminus1 * deltaN;
-
-    if (std::isinf(hn) || std::isnan(hn)) {
-      return false;
-    }
-    if (std::fabs(deltaN - 1.0) < 10e-15) {
-      break;
-    }
-
-    // Update c, d, h
-    dnminus1 = dn;
-    cnminus1 = cn;
-    hnminus1 = hn;
-    iterationCount++;
-  }
-
-  if (iterationCount >= maxNumberOfIterations) {
-    return false;
-  }
-
-  *result = hn;
-  return true;
-}
-
-typedef double (*TermUpdater)(double previousTerm, double index, double s, double x);
-static inline bool InfiniteSeriesEvaluation(double firstTerm, TermUpdater termUpdater, double termLimit, int maxNumberOfIterations, double * result, double context1, double context2) {
-  double iterationCount = 0.0;
-  double currentTerm = firstTerm;
-  double sum = currentTerm;
-  while (iterationCount < maxNumberOfIterations
-      && !std::isinf(sum)
-      && std::fabs(currentTerm/sum) > termLimit)
-  {
-    iterationCount+= 1.0;
-    currentTerm = termUpdater(currentTerm, iterationCount, context1, context2);
-    sum+= currentTerm;
-  }
-  if (iterationCount >= maxNumberOfIterations) {
-    return false;
-  }
-  *result = sum;
-  return true;
-}
-
-
 bool regularizedGamma(double s, double x, double epsilon, int maxNumberOfIterations, double * result) {
   // TODO Put interruption instead of maxNumberOfIterations
 
@@ -125,7 +48,7 @@ bool regularizedGamma(double s, double x, double epsilon, int maxNumberOfIterati
      * Other method : Lentz's method with Thompson and Barnett's modification */
 
     double continuedFractionValue = 0.0;
-    if (!ContinuedFractionEvaluation(
+    if (!Helper::ContinuedFractionEvaluation(
           [](double index, double s, double x) { return index * (s - index); },
           [](double index, double s, double x) { return (2.0 * index + 1.0) + x - s; },
           maxNumberOfIterations,
@@ -142,7 +65,7 @@ bool regularizedGamma(double s, double x, double epsilon, int maxNumberOfIterati
   /* For x < s + 1: Compute using the infinite series representation
    * incompleteGamma(s,x) = 1/gamma(a) * e^-x * x^s * sum(x^n/((s+n)(s+n-1)...) */
   double infiniteSeriesValue = 0.0;
-  if (!InfiniteSeriesEvaluation(
+  if (!Helper::InfiniteSeriesEvaluation(
         1.0/s,
         [](double previousTerm, double index, double s, double x) { return previousTerm * x / (s + index); },
         epsilon,