From 4c342cd93305ae36f41d774b799fd9b7807905be Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9a=20Saviot?= Date: Tue, 13 Aug 2019 17:59:00 +0200 Subject: [PATCH] [apps/probability] Helper for contonued fractions and infinite series --- apps/probability/Makefile | 5 +- apps/probability/law/chi_squared_law.cpp | 1 - apps/probability/law/helper.cpp | 79 ++++++++++++++++++++ apps/probability/law/helper.h | 17 +++++ apps/probability/law/regularized_gamma.cpp | 83 +--------------------- 5 files changed, 101 insertions(+), 84 deletions(-) create mode 100644 apps/probability/law/helper.cpp create mode 100644 apps/probability/law/helper.h diff --git a/apps/probability/Makefile b/apps/probability/Makefile index 705d1aea7..17573d725 100644 --- a/apps/probability/Makefile +++ b/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\ diff --git a/apps/probability/law/chi_squared_law.cpp b/apps/probability/law/chi_squared_law.cpp index 3ad7db661..c212fa77d 100644 --- a/apps/probability/law/chi_squared_law.cpp +++ b/apps/probability/law/chi_squared_law.cpp @@ -1,6 +1,5 @@ #include "chi_squared_law.h" #include "regularized_gamma.h" -#include #include namespace Probability { diff --git a/apps/probability/law/helper.cpp b/apps/probability/law/helper.cpp new file mode 100644 index 000000000..8f6f565c1 --- /dev/null +++ b/apps/probability/law/helper.cpp @@ -0,0 +1,79 @@ +#include "helper.h" +#include "law.h" +#include +#include +#include + +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; +} diff --git a/apps/probability/law/helper.h b/apps/probability/law/helper.h new file mode 100644 index 000000000..83f8924f8 --- /dev/null +++ b/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 + diff --git a/apps/probability/law/regularized_gamma.cpp b/apps/probability/law/regularized_gamma.cpp index 208632c8b..da210fd94 100644 --- a/apps/probability/law/regularized_gamma.cpp +++ b/apps/probability/law/regularized_gamma.cpp @@ -1,87 +1,10 @@ #include "regularized_gamma.h" #include "law.h" +#include "helper.h" #include #include #include -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,