diff --git a/src/stan/agrad/rev/gamma_q.hpp b/src/stan/agrad/rev/gamma_q.hpp
index b1e69418781..ba0ca04cdc0 100644
--- a/src/stan/agrad/rev/gamma_q.hpp
+++ b/src/stan/agrad/rev/gamma_q.hpp
@@ -7,6 +7,7 @@
 #include <stan/agrad/rev/op/dv_vari.hpp>
 #include <stan/agrad/rev/op/vd_vari.hpp>
 #include <stan/math/functions/gamma_q.hpp>
+#include <stan/prob/internal_math.hpp>
 #include <boost/math/special_functions/gamma.hpp>
 #include <boost/math/special_functions/digamma.hpp>
 
@@ -21,31 +22,15 @@ namespace stan {
                      avi,bvi) {
         }
         void chain() {
-
-          double u = stan::math::gamma_q(avi_->val_, bvi_->val_);
-      
-          double S = 0;
-          double s = 1;
-          double l = std::log(bvi_->val_);
-          double g = boost::math::tgamma(avi_->val_);
-          double dig = boost::math::digamma(avi_->val_);
-      
-          int k = 0;
-          double delta = s / (avi_->val_ * avi_->val_);
-      
-          while (std::fabs(delta) > 1e-6) {
-            S += delta;
-            ++k;
-            s *= - bvi_->val_ / k;
-            delta = s / ((k + avi_->val_) * (k + avi_->val_));
-          }
-          
-
-          avi_->adj_ += adj_ * ((1.0 - u) * ( dig - l ) + std::exp( avi_->val_ * l ) * S / g);
-          bvi_->adj_ -= adj_ * (std::exp(-bvi_->val_) * std::pow(bvi_->val_, avi_->val_ - 1.0) / g);
+          avi_->adj_ += adj_ 
+            * stan::math::gradRegIncGamma(avi_->val_, bvi_->val_,
+                                          boost::math::tgamma(avi_->val_),
+                                          boost::math::digamma(avi_->val_));
+          bvi_->adj_ -= adj_ 
+            * boost::math::gamma_p_derivative(avi_->val_, bvi_->val_);
         }
       };
-
+      
       class gamma_q_vd_vari : public op_vd_vari {
       public:
         gamma_q_vd_vari(vari* avi, double b) :
@@ -53,26 +38,10 @@ namespace stan {
                      avi,b) {
         }
         void chain() {
-
-          double u = stan::math::gamma_q(avi_->val_, bd_);
-      
-          double S = 0;
-          double s = 1;
-          double l = std::log(bd_);
-          double g = boost::math::tgamma(avi_->val_);
-          double dig = boost::math::digamma(avi_->val_);
-      
-          int k = 0;
-          double delta = s / (avi_->val_ * avi_->val_);
-      
-          while (std::fabs(delta) > 1e-6) {
-            S += delta;
-            ++k;
-            s *= - bd_ / k;
-            delta = s / ((k + avi_->val_) * (k + avi_->val_));
-          }
-
-          avi_->adj_ += adj_ * ((1.0 - u) * ( dig - l ) + std::exp( avi_->val_ * l ) * S / g);
+          avi_->adj_ += adj_ 
+            * stan::math::gradRegIncGamma(avi_->val_, bd_,
+                                          boost::math::tgamma(avi_->val_),
+                                          boost::math::digamma(avi_->val_));
         }
       };
 
@@ -83,7 +52,8 @@ namespace stan {
                      a,bvi) {
         }
         void chain() {
-          bvi_->adj_ -= adj_ * (std::exp(-bvi_->val_) * std::pow(bvi_->val_, ad_ - 1.0) / boost::math::tgamma(ad_));     
+          bvi_->adj_ -= adj_ 
+            * boost::math::gamma_p_derivative(ad_, bvi_->val_);
         }
       };
     }
diff --git a/src/stan/prob/internal_math.hpp b/src/stan/prob/internal_math.hpp
index 5315eaa6a8f..7ccb5e223b5 100644
--- a/src/stan/prob/internal_math.hpp
+++ b/src/stan/prob/internal_math.hpp
@@ -4,6 +4,7 @@
 #include <math.h>
 #include <boost/math/special_functions/gamma.hpp>
 #include <boost/math/special_functions/beta.hpp>
+#include <boost/math/special_functions/fpclassify.hpp>
 
 namespace stan {
     
@@ -173,33 +174,31 @@ namespace stan {
           
       }
       
-      // Gradient of the regularized incomplete gamma functions igamma(a, g)
-      double gradRegIncGamma(double a, double z, double g, double dig, double precision = 1e-6) 
-      {
-          
-          using boost::math::gamma_p;
-          
-          double S = 0;
-          double s = 1;
-          double l = std::log(z);
-          
-          int k = 0;
-          double delta = s / (a * a);
-          
-          while (fabs(delta) > precision)
-          {
-              S += delta;
-              ++k;
-              s *= - z / k;
-              delta = s / ((k + a) * (k + a));
-          }
-          
-          // Precomputed values
-          // dig -> digamma(a)
-          // g   -> g(a)
-          return gamma_p(a, z) * ( dig - l ) + std::exp( a * l ) * S / g;
-          
+    // Gradient of the regularized incomplete gamma functions igamma(a, g)
+    // Precomputed values
+    // g   = boost::math::tgamma(a)
+    // dig = boost::math::digamma(a)
+    double gradRegIncGamma(double a, double z, double g, double dig, 
+                           double precision = 1e-6) {
+      using boost::math::gamma_p;
+          
+      double S = 0;
+      double s = 1;
+      double l = std::log(z);
+          
+      int k = 0;
+      double delta = s / (a * a);
+          
+      while (fabs(delta) > precision) {
+        S += delta;
+        ++k;
+        s *= - z / k;
+        delta = s / ((k + a) * (k + a));
+        if (boost::math::isinf(delta))
+          throw std::domain_error("stan::math::gradRegIncGamma not converging");
       }
+      return gamma_p(a, z) * ( dig - l ) + std::exp( a * l ) * S / g;
+    }
       
   }
 
diff --git a/src/test/agrad/rev/gamma_q_test.cpp b/src/test/agrad/rev/gamma_q_test.cpp
index a0a63dda69c..09e43daea7e 100644
--- a/src/test/agrad/rev/gamma_q_test.cpp
+++ b/src/test/agrad/rev/gamma_q_test.cpp
@@ -21,6 +21,26 @@ TEST(AgradRev,gamma_q_var_var) {
   b = -1.0;
   EXPECT_THROW(gamma_q(a,b), std::domain_error);
 }
+TEST(AgradRevGammaQ,infLoopInVersion2_0_1_var_var) {
+  // FIXME: causes infinite loop in 2.0.1 gradient calcs
+  AVAR a = 8.01006;
+  AVAR b = 2.47579e+215;
+  AVEC x = createAVEC(a,b);
+
+  AVAR f = gamma_q(a,b);
+  VEC g;
+  EXPECT_THROW(f.grad(x,g), std::domain_error);
+}
+TEST(AgradRevGammaQ,infLoopInVersion2_0_1_var_double) {
+  // FIXME: causes infinite loop in 2.0.1 gradient calcs
+  AVAR a = 8.01006;
+  double b = 2.47579e+215;
+  AVEC x = createAVEC(a);
+
+  AVAR f = gamma_q(a,b);
+  VEC g;
+  EXPECT_THROW(f.grad(x,g), std::domain_error);
+}
 TEST(AgradRev,gamma_q_double_var) {
   double a = 0.5;
   AVAR b = 1.0;
diff --git a/src/test/prob/internal_math_test.cpp b/src/test/prob/internal_math_test.cpp
new file mode 100644
index 00000000000..2067b584759
--- /dev/null
+++ b/src/test/prob/internal_math_test.cpp
@@ -0,0 +1,21 @@
+#include <gtest/gtest.h>
+#include <stan/prob/internal_math.hpp>
+
+TEST(ProbInternalMath, gradRegIncGamma_typical) {
+  double a = 0.5;
+  double b = 1.0;
+  double g = 1.77245;
+  double dig = -1.96351;
+  
+  EXPECT_FLOAT_EQ(0.38984156, stan::math::gradRegIncGamma(a, b, g, dig));
+}
+
+TEST(ProbInternalMath, gradRegIncGamma_infLoopInVersion2_0_1) {
+  double a = 8.01006;
+  double b = 2.47579e+215;
+  double g = 5143.28;
+  double dig = 2.01698;
+  
+  EXPECT_THROW(stan::math::gradRegIncGamma(a, b, g, dig),
+               std::domain_error);
+}