Initial implementation of igamma and igammac.
diff --git a/test/array.cpp b/test/array.cpp
index 96aef31..a37874c 100644
--- a/test/array.cpp
+++ b/test/array.cpp
@@ -295,7 +295,6 @@
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 @@
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 @@
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