From 7ea35bfa1c0b4950feae65d49c0e6f2cbf3691d9 Mon Sep 17 00:00:00 2001
From: Eugene Brevdo <ebrevdo@gmail.com>
Date: Thu, 3 Mar 2016 19:39:41 -0800
Subject: Initial implementation of igamma and igammac.

---
 test/array.cpp | 47 ++++++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 46 insertions(+), 1 deletion(-)

(limited to 'test/array.cpp')

diff --git a/test/array.cpp b/test/array.cpp
index 96aef31c7..a37874cc2 100644
--- a/test/array.cpp
+++ b/test/array.cpp
@@ -295,7 +295,6 @@ template<typename ArrayType> void array_real(const ArrayType& m)
   VERIFY_IS_APPROX(Eigen::pow(m1,2*exponents), m1.square().square());
   VERIFY_IS_APPROX(m1.pow(2*exponents), m1.square().square());
   VERIFY_IS_APPROX(pow(m1(0,0), exponents), ArrayType::Constant(rows,cols,m1(0,0)*m1(0,0)));
-  
 
   VERIFY_IS_APPROX(m3.pow(RealScalar(0.5)), m3.sqrt());
   VERIFY_IS_APPROX(pow(m3,RealScalar(0.5)), m3.sqrt());
@@ -305,6 +304,14 @@ template<typename ArrayType> void array_real(const ArrayType& m)
 
   VERIFY_IS_APPROX(log10(m3), log(m3)/log(10));
 
+  // Smoke test to check any compilation issues
+  ArrayType m1_abs_p1 = m1.abs() + 1;
+  ArrayType m2_abs_p1 = m2.abs() + 1;
+  VERIFY_IS_APPROX(Eigen::igamma(m1_abs_p1, m2_abs_p1), Eigen::igamma(m1_abs_p1, m2_abs_p1));
+  VERIFY_IS_APPROX(Eigen::igammac(m1_abs_p1, m2_abs_p1), Eigen::igammac(m1_abs_p1, m2_abs_p1));
+  VERIFY_IS_APPROX(Eigen::igamma(m2_abs_p1, m1_abs_p1), Eigen::igamma(m2_abs_p1, m1_abs_p1));
+  VERIFY_IS_APPROX(Eigen::igammac(m2_abs_p1, m1_abs_p1), Eigen::igammac(m2_abs_p1, m1_abs_p1));
+
   // scalar by array division
   const RealScalar tiny = sqrt(std::numeric_limits<RealScalar>::epsilon());
   s1 += Scalar(tiny);
@@ -323,6 +330,44 @@ template<typename ArrayType> void array_real(const ArrayType& m)
                     std::numeric_limits<RealScalar>::infinity());
     VERIFY_IS_EQUAL(numext::digamma(Scalar(-1)),
                     std::numeric_limits<RealScalar>::infinity());
+
+    Scalar a_s[] = {Scalar(0), Scalar(1), Scalar(1.5), Scalar(4), Scalar(0.0001), Scalar(10000.5)};
+    Scalar x_s[] = {Scalar(0), Scalar(1), Scalar(1.5), Scalar(4), Scalar(0.0001), Scalar(10000.5)};
+
+    // location i*6+j corresponds to a_s[i], x_s[j].
+    Scalar nan = std::numeric_limits<Scalar>::quiet_NaN();
+    Scalar igamma_s[][6] = {
+      {0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+      {0.0, 0.6321205588285578, 0.7768698398515702, 0.9816843611112658,
+       9.999500016666262e-05, 1.0},
+      {0.0, 0.4275932955291202, 0.608374823728911, 0.9539882943107686,
+       7.522076445089201e-07, 1.0},
+      {0.0, 0.01898815687615381, 0.06564245437845008, 0.5665298796332909,
+       4.166333347221828e-18, 1.0},
+      {0.0, 0.9999780593618628, 0.9999899967080838, 0.9999996219837988,
+       0.9991370418689945, 1.0},
+      {0.0, 0.0, 0.0, 0.0, 0.0, 0.5013297751014064}};
+    Scalar igammac_s[][6] = {
+      {1.0, 1.0, 1.0, 1.0, 1.0, 1.0},
+      {1.0, 0.36787944117144233, 0.22313016014842982,
+       0.018315638888734182, 0.9999000049998333, 0.0},
+      {1.0, 0.5724067044708798, 0.3916251762710878,
+       0.04601170568923136, 0.9999992477923555, 0.0},
+      {1.0, 0.9810118431238462, 0.9343575456215499,
+       0.4334701203667089, 1.0, 0.0},
+      {1.0, 2.1940638138146658e-05, 1.0003291916285e-05,
+       3.7801620118431334e-07, 0.0008629581310054535, 0.0},
+      {1.0, 1.0, 1.0, 1.0, 1.0, 0.49867022490946517}};
+    for (int i = 0; i < 6; ++i) {
+      for (int j = 0; j < 6; ++j) {
+        //std::cout << numext::igamma(a_s[i], x_s[j]) << " vs. " << igamma_s[i][j] << std::endl;
+        //std::cout << numext::igammac(a_s[i], x_s[j]) << " c.vs. " <<
+        //igammac_s[i][j] << std::endl;
+        std::cout << a_s[i] << ", " << x_s[j] << std::endl;
+        VERIFY_IS_APPROX(numext::igamma(a_s[i], x_s[j]), igamma_s[i][j]);
+        VERIFY_IS_APPROX(numext::igammac(a_s[i], x_s[j]), igammac_s[i][j]);
+      }
+    }
   }
 #endif  // EIGEN_HAS_C99_MATH
 
-- 
cgit v1.2.3