Improve robustness of igamma and igammac to bad inputs.

Check for nan inputs and propagate them immediately. Limit the number of internal iterations to 2000 (same number as used by scipy.special.gammainc). This prevents an infinite loop when the function is called with nan or very large arguments.

Original change by mfirgunov@google.com

1 files changed, 10 insertions, 2 deletions

diff --git a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h
index 5d1b8fcc2..7deca3e12 100644
--- a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h
+++ b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h
@@ -535,6 +535,10 @@ struct igammac_impl {
       return nan;
     }
 
+    if (numext::isnan(a) || numext::isnan(x)) { // propagate nans
+      return nan;
+    }
+
     if ((x < one) || (x < a)) {
       /* The checks above ensure that we meet the preconditions for
        * igamma_impl::Impl(), so call it, rather than igamma_impl::Run().
@@ -592,7 +596,7 @@ struct igammac_impl {
     qkm1 = z * x;
     ans = pkm1 / qkm1;
 
-    while (true) {
+    for (int i = 0; i < 2000; i++) {
       c += one;
       y += one;
       z += two;
@@ -724,6 +728,10 @@ struct igamma_impl {
       return nan;
     }
 
+    if (numext::isnan(a) || numext::isnan(x)) { // propagate nans
+      return nan;
+    }
+
     if ((x > one) && (x > a)) {
       /* The checks above ensure that we meet the preconditions for
        * igammac_impl::Impl(), so call it, rather than igammac_impl::Run().
@@ -770,7 +778,7 @@ struct igamma_impl {
     c = one;
     ans = one;
 
-    while (true) {
+    for (int i = 0; i < 2000; i++) {
       r += one;
       c *= x/r;
       ans += c;