From f6481533626a35b71ddc0a2b0f0a20c2ac3bdf28 Mon Sep 17 00:00:00 2001
From: Heiko Bauke <heiko.bauke@mail.de>
Date: Tue, 27 Feb 2024 00:38:11 +0100
Subject: [PATCH] speed-up unit tests by reducing number of test values

while increasing the covered numerical range, replace implementation of inverse cdf for better numerical stability
---
 tests/test_distributions.cc | 45 ++++++++++++++++++++++++++-----------
 trng/chi_square_dist.hpp    | 22 +-----------------
 2 files changed, 33 insertions(+), 34 deletions(-)

diff --git a/tests/test_distributions.cc b/tests/test_distributions.cc
index 019da47..4312c8d 100644
--- a/tests/test_distributions.cc
+++ b/tests/test_distributions.cc
@@ -37,6 +37,7 @@
 #include <numeric>
 #include <ciso646>
 #include <sstream>
+#include <algorithm>
 
 #include <catch2/catch.hpp>
 
@@ -133,7 +134,6 @@ double chi_percentil(const std::vector<double> &p, const std::vector<int> &count
 template<typename dist>
 void continuous_dist_test_integrate_pdf(const dist &d) {
   using result_type = typename dist::result_type;
-  // const int samples{1024 * 16 + 1};
   const int samples(static_cast<int>(std::min(
       1024ll * 1024, static_cast<long long>(std::round(
                          1 / std::sqrt(std::numeric_limits<result_type>::epsilon()))))));
@@ -152,14 +152,19 @@ void continuous_dist_test_integrate_pdf(const dist &d) {
 template<typename dist>
 void continuous_dist_test_icdf(const dist &d) {
   using result_type = typename dist::result_type;
-  const int bins{1024 * 1024};
-  const result_type dp{result_type(1) / result_type(bins)};
   const auto eps{256 * std::numeric_limits<result_type>::epsilon()};
-  for (int i{1}; i < bins; ++i) {
-    const result_type p{i * dp};
+  std::vector<result_type> p_values;
+  p_values.push_back(result_type(1) / result_type(2));
+  for (int i{2}; i < std::numeric_limits<result_type>::digits; ++i) {
+    p_values.push_back(result_type(0) + std::pow(p_values[0], i));
+    p_values.push_back(result_type(1) - std::pow(p_values[0], i));
+  }
+  std::sort(p_values.begin(), p_values.end());
+  for (const auto p : p_values) {
     const result_type x{d.icdf(p)};
     const result_type y{d.cdf(x)};
-    if (not(std::abs(y - p) < eps) or i == 1)
+    if (not(std::abs(y - p) <
+            eps))  // switch into REQUIRE macro in failure case only, for performance reasons
       REQUIRE(std::abs(y - p) < eps);
   }
 }
@@ -203,10 +208,18 @@ void continuous_dist_test_streamable(dist &d) {
 
 template<typename T>
 void continuous_dist_test(T &d) {
-  SECTION("integrate pdf") { continuous_dist_test_integrate_pdf(d); }
-  SECTION("icdf") { continuous_dist_test_icdf(d); }
-  SECTION("chi2 test") { continuous_dist_test_chi2_test(d); }
-  SECTION("streamable") { continuous_dist_test_streamable(d); }
+  SECTION("integrate pdf") {
+    continuous_dist_test_integrate_pdf(d);
+  }
+  SECTION("icdf") {
+    continuous_dist_test_icdf(d);
+  }
+  SECTION("chi2 test") {
+    continuous_dist_test_chi2_test(d);
+  }
+  SECTION("streamable") {
+    continuous_dist_test_streamable(d);
+  }
 }
 
 
@@ -285,9 +298,15 @@ void discrete_dist_test_streamable(dist &d) {
 
 template<typename T>
 void discrete_dist_test(T &d) {
-  SECTION("test pdf & cdf") { discrete_dist_test_pdf(d); }
-  SECTION("chi2_test") { discrete_dist_test_chi2_test(d); }
-  SECTION("streamable") { discrete_dist_test_streamable(d); }
+  SECTION("test pdf & cdf") {
+    discrete_dist_test_pdf(d);
+  }
+  SECTION("chi2_test") {
+    discrete_dist_test_chi2_test(d);
+  }
+  SECTION("streamable") {
+    discrete_dist_test_streamable(d);
+  }
 }
 
 
diff --git a/trng/chi_square_dist.hpp b/trng/chi_square_dist.hpp
index 9e0b358..43d5ffc 100644
--- a/trng/chi_square_dist.hpp
+++ b/trng/chi_square_dist.hpp
@@ -112,27 +112,7 @@ namespace trng {
     // inverse cumulative density function
     TRNG_CUDA_ENABLE
     result_type icdf_(result_type x) const {
-      if (x <= math::numeric_limits<result_type>::epsilon())
-        return 0;
-      const result_type kappa{P.nu() / result_type(2)};
-      const result_type theta{2};
-      if (kappa == 1)  // special case of exponential distribution
-        return -math::ln(1 - x) * theta;
-      const result_type ln_Gamma_kappa{math::ln_Gamma(kappa)};
-      result_type y{kappa}, y_old;
-      if (kappa < 1 and x < result_type(1) / result_type(2))
-        y = x * x;
-      int num_iterations{0};
-      do {
-        ++num_iterations;
-        y_old = y;
-        const result_type f0{math::GammaP(kappa, y) - x};
-        const result_type f1{math::pow(y, kappa - 1) * math::exp(-y - ln_Gamma_kappa)};
-        const result_type f2{f1 * (kappa - 1 - y) / y};
-        y -= f0 / f1 * (1 + f0 * f2 / (2 * f1 * f1));
-      } while (num_iterations < 16 and
-               math::abs((y - y_old) / y) > 16 * math::numeric_limits<result_type>::epsilon());
-      return y * theta;
+      return math::inv_GammaP(P.nu() / result_type(2), x) * 2;
     }
 
   public: