From e900b010c899b04b1710c5e99486b52494a40ea1 Mon Sep 17 00:00:00 2001
From: Rasmus Munk Larsen <rmlarsen@google.com>
Date: Mon, 19 Mar 2018 09:04:54 -0700
Subject: [PATCH] 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
---
 .../src/SpecialFunctions/SpecialFunctionsImpl.h      | 12 ++++++++++--
 1 file 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;