diff --git a/trng/gamma_dist.hpp b/trng/gamma_dist.hpp
index 057b1ee..91decf5 100644
--- a/trng/gamma_dist.hpp
+++ b/trng/gamma_dist.hpp
@@ -118,23 +118,7 @@ namespace trng {
     // inverse cumulative density function
     TRNG_CUDA_ENABLE
     result_type icdf_(result_type x) const {
-      if (x <= math::numeric_limits<result_type>::epsilon())
-        return 0;
-      if (P.kappa() == 1)  // special case of exponential distribution
-        return -math::ln(1 - x) * P.theta();
-      const result_type ln_Gamma_kappa{math::ln_Gamma(P.kappa())};
-      result_type y{P.kappa()}, y_old;
-      int num_iterations{0};
-      do {
-        ++num_iterations;
-        y_old = y;
-        const result_type f0{math::GammaP(P.kappa(), y) - x};
-        const result_type f1{math::exp((P.kappa() - 1) * math::ln(y) - y - ln_Gamma_kappa)};
-        const result_type f2{f1 * (P.kappa() - 1 - y) / y};
-        y -= f0 / f1 * (1 + f0 * f2 / (2 * f1 * f1));
-      } while (num_iterations < 16 &&
-               math::abs((y - y_old) / y) > 16 * math::numeric_limits<result_type>::epsilon());
-      return y * P.theta();
+      return math::inv_GammaP(P.kappa(), x) * P.theta();
     }
 
   public:
@@ -202,7 +186,7 @@ namespace trng {
         return 0;
       if (x == 1)
         return math::numeric_limits<result_type>::infinity();
-      return icdf_(x);
+      return math::inv_GammaP(P.kappa(), x) * P.theta();
     }
   };
 
diff --git a/trng/special_functions.hpp b/trng/special_functions.hpp
index 9258cbe..2587f01 100644
--- a/trng/special_functions.hpp
+++ b/trng/special_functions.hpp
@@ -586,10 +586,10 @@ namespace trng {
       template<typename T>
       TRNG_CUDA_ENABLE T inv_GammaP(T a, T p) {
         const T eps{sqrt(numeric_limits<T>::epsilon())};
-        T a1{a - 1};
-        T glna{ln_Gamma(a)};
-        T lna1{ln(a1)};
-        T afac{exp(a1 * (lna1 - 1) - glna)};
+        const T a1{a - 1};
+        const T glna{ln_Gamma(a)};
+        const T lna1{ln(a1)};
+        const T afac{exp(a1 * (lna1 - 1) - glna)};
         T x;
         // initial guess
         if (a > T{1}) {
@@ -603,7 +603,7 @@ namespace trng {
           x = p < t ? pow(p / t, 1 / a) : 1 - ln1p(-(p - t) / (1 - t));
         }
         // refinement by Halley's method
-        for (int i{0}; i < 32; ++i) {
+        for (int i{0}; i < numeric_limits<T>::digits; ++i) {
           if (x <= 0) {
             x = 0;
             break;