diff --git a/doc/sf_and_dist/igamma.qbk b/doc/sf_and_dist/igamma.qbk
index ba5c69a6d..d246307c1 100644
--- a/doc/sf_and_dist/igamma.qbk
+++ b/doc/sf_and_dist/igamma.qbk
@@ -251,14 +251,14 @@ The lower incomplete gamma is computed from its series representation:
 4) [equation igamma8]
 
 Or by subtraction of the upper integral from either [Gamma](a) or 1
-when /x > a and x > 1.1/.
+when /x - (1/(3x)) > a and x > 1.1/.
 
 The upper integral is computed from Legendre's continued fraction representation:
 
 5) [equation igamma9]
 
-When /x > 1.1/ or by subtraction of the lower integral from either [Gamma](a) or 1
-when /x < a/.
+When /(x > 1.1)/ or by subtraction of the lower integral from either [Gamma](a) or 1
+when /x - (1/(3x))  < a/.
 
 For /x < 1.1/ computation of the upper integral is more complex as the continued 
 fraction representation is unstable in this area.  However there is another 
@@ -279,9 +279,9 @@ That therefore imposes a similar limit on the precision of this
 function in this region.
 
 For /x < 1.1/ the crossover point where the result is ~0.5 no longer
-occurs for /x ~ y/.  Using /x * 1.1 < a/ as the crossover criterion
+occurs for /x ~ y/.  Using /x * 0.75 < a/ as the crossover criterion
 for /0.5 < x <= 1.1/ keeps the maximum value computed (whether
-it's the upper or lower interval) to around 0.6.   Likewise for
+it's the upper or lower interval) to around 0.75.   Likewise for
 /x <= 0.5/ then using /-0.4 / log(x) < a/ as the crossover criterion
 keeps the maximum value computed to around 0.7
 (whether it's the upper or lower interval).
diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp
index 44232936d..661165bf7 100644
--- a/include/boost/math/special_functions/gamma.hpp
+++ b/include/boost/math/special_functions/gamma.hpp
@@ -733,22 +733,26 @@ T regularised_gamma_prefix(T a, T z, const Policy& pol, const lanczos::undefined
 // Upper gamma fraction for very small a:
 //
 template <class T, class Policy>
-inline T tgamma_small_upper_part(T a, T x, const Policy& pol)
+inline T tgamma_small_upper_part(T a, T x, const Policy& pol, T* pgam = 0)
 {
    BOOST_MATH_STD_USING  // ADL of std functions.
    //
    // Compute the full upper fraction (Q) when a is very small:
    //
    T result;
-   result = boost::math::tgamma1pm1(a, pol) - boost::math::powm1(x, a, pol);
+   result = boost::math::tgamma1pm1(a, pol);
+   if(pgam)
+      *pgam = (result + 1) / a;
+   T p = boost::math::powm1(x, a, pol);
+   result -= p;
    result /= a;
    detail::small_gamma2_series<T> s(a, x);
    boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
 #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
    T zero = 0;
-   result -= pow(x, a) * tools::sum_series(s, boost::math::policies::digits<T, Policy>(), max_iter, zero);
+   result -= (p + 1) * tools::sum_series(s, boost::math::policies::digits<T, Policy>(), max_iter, zero);
 #else
-   result -= pow(x, a) * tools::sum_series(s, boost::math::policies::digits<T, Policy>(), max_iter);
+   result -= (p + 1) * tools::sum_series(s, boost::math::policies::digits<T, Policy>(), max_iter);
 #endif
    policies::check_series_iterations("boost::math::tgamma_small_upper_part<%1%>(%1%, %1%)", max_iter, pol);
    return result;
@@ -877,9 +881,10 @@ T gamma_incomplete_imp(T a, T x, bool normalised, bool invert,
       {
          // Compute Q:
          invert = !invert;
-         result = tgamma_small_upper_part(a, x, pol);
+         T g;
+         result = tgamma_small_upper_part(a, x, pol, &g);
          if(normalised)
-            result /= boost::math::tgamma(a, pol);
+            result /= g;
          if(p_derivative)
             *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
       }
@@ -887,9 +892,9 @@ T gamma_incomplete_imp(T a, T x, bool normalised, bool invert,
    else if(x < 1.1)
    {
       //
-      // Changover here occurs when P ~ 0.6 or Q ~ 0.4:
+      // Changover here occurs when P ~ 0.75 or Q ~ 0.25:
       //
-      if(x * 1.1f < a)
+      if(x * 0.75f < a)
       {
          // Compute P:
          result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol);
@@ -902,9 +907,10 @@ T gamma_incomplete_imp(T a, T x, bool normalised, bool invert,
       {
          // Compute Q:
          invert = !invert;
-         result = tgamma_small_upper_part(a, x, pol);
+         T g;
+         result = tgamma_small_upper_part(a, x, pol, &g);
          if(normalised)
-            result /= boost::math::tgamma(a, pol);
+            result /= g;
          if(p_derivative)
             *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
       }
@@ -981,11 +987,14 @@ T gamma_incomplete_imp(T a, T x, bool normalised, bool invert,
          // Regular case where the result will not be too close to 0.5.
          //
          // Changeover here occurs at P ~ Q ~ 0.5
+         // Note that series computation of P is about x2 faster than continued fraction
+         // calculation of Q, so try and use the CF only when really necessary, especially
+         // for small x.
          //
          result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol);
          if(p_derivative)
             *p_derivative = result;
-         if(x < a)
+         if(x - (1 / (3 * x)) < a)
          {
             // Compute P:
             if(result != 0)
@@ -1001,6 +1010,8 @@ T gamma_incomplete_imp(T a, T x, bool normalised, bool invert,
       }
    }
 
+   if(normalised && (result > 1))
+      result = 1;
    if(invert)
    {
       T gam = normalised ? 1 : boost::math::tgamma(a, pol);
diff --git a/minimax/f.cpp b/minimax/f.cpp
index 6850d27ab..f35c7e40d 100644
--- a/minimax/f.cpp
+++ b/minimax/f.cpp
@@ -4,8 +4,8 @@
 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
 
 #define L22
-#include "../tools/ntl_rr_lanczos.hpp"
-#include "../tools/ntl_rr_digamma.hpp"
+//#include "../tools/ntl_rr_lanczos.hpp"
+//#include "../tools/ntl_rr_digamma.hpp"
 #include <boost/math/bindings/rr.hpp>
 #include <boost/math/tools/polynomial.hpp>
 #include <boost/math/special_functions.hpp>
diff --git a/minimax/main.cpp b/minimax/main.cpp
index 99df43976..ca7a5071e 100644
--- a/minimax/main.cpp
+++ b/minimax/main.cpp
@@ -409,11 +409,11 @@ void do_test_n(T, const char* name, unsigned count)
       std::string msg = "Max Error found at ";
       msg += name;
       msg += " precision = ";
-      msg.append(62 - 17 - msg.size(), ' ');
+      //msg.append(62 - 17 - msg.size(), ' ');
       std::cout << msg << "Poly: " << std::setprecision(6) 
-         << std::setw(15) << std::left 
+         //<< std::setw(15) << std::left 
          << boost::math::tools::real_cast<T>(max_error) 
-         << "Cheb: " << boost::math::tools::real_cast<T>(max_cheb_error) << std::endl;
+         << " Cheb: " << boost::math::tools::real_cast<T>(max_cheb_error) << std::endl;
    }
    else
    {
diff --git a/performance/test_igamma.cpp b/performance/test_igamma.cpp
index b81de2bfa..64b989e30 100644
--- a/performance/test_igamma.cpp
+++ b/performance/test_igamma.cpp
@@ -188,3 +188,124 @@ BOOST_MATH_PERFORMANCE_TEST(igamma_test, "igamma-gsl")
 }
 
 #endif
+
+#ifdef TEST_DCDFLIB
+#include <dcdflib.h>
+namespace dcd{
+
+inline double gamma_q(double x, double y)
+{ 
+   double ans, qans;
+   int i = 0;
+   gamma_inc (&x, &y, &ans, &qans, &i); 
+   return qans;
+}
+
+inline double gamma_p(double x, double y)
+{ 
+   double ans, qans;
+   int i = 0;
+   gamma_inc (&x, &y, &ans, &qans, &i); 
+   return ans;
+}
+
+inline double gamma_q_inv(double x, double y)
+{ 
+   double ans, p, nul;
+   int i = 0;
+   p = 1 - y;
+   nul = 0;
+   gamma_inc_inv (&x, &ans, &nul, &p, &y, &i); 
+   return ans;
+}
+
+inline double gamma_p_inv(double x, double y)
+{ 
+   double ans, p, nul;
+   int i = 0;
+   p = 1 - y;
+   nul = 0;
+   gamma_inc_inv (&x, &ans, &nul, &y, &p, &i); 
+   return ans;
+}
+
+}
+
+template <std::size_t N>
+double igamma_evaluate2_dcd(const boost::array<boost::array<T, 6>, N>& data)
+{
+   double result = 0;
+   for(unsigned i = 0; i < N; ++i)
+   {
+      result += dcd::gamma_p(data[i][0], data[i][1]);
+      result += dcd::gamma_q(data[i][0], data[i][1]);
+   }
+   return result;
+}
+
+BOOST_MATH_PERFORMANCE_TEST(igamma_test, "igamma-dcd")
+{
+   double result = igamma_evaluate2_dcd(igamma_big_data);
+   result += igamma_evaluate2_dcd(igamma_int_data);
+   result += igamma_evaluate2_dcd(igamma_med_data);
+   result += igamma_evaluate2_dcd(igamma_small_data);
+
+   consume_result(result);
+   set_call_count(
+      2 * (sizeof(igamma_big_data) 
+      + sizeof(igamma_int_data) 
+      + sizeof(igamma_med_data)
+      + sizeof(igamma_small_data)) / sizeof(igamma_big_data[0]));
+}
+
+template <std::size_t N>
+double igamma_inv_evaluate2_dcd(const boost::array<boost::array<T, 6>, N>& data)
+{
+   double result = 0;
+   for(unsigned i = 0; i < N; ++i)
+   {
+      result += dcd::gamma_p_inv(data[i][0], data[i][5]);
+      result += dcd::gamma_q_inv(data[i][0], data[i][3]);
+   }
+   return result;
+}
+
+BOOST_MATH_PERFORMANCE_TEST(igamma_inv_test, "igamma_inv-dcd")
+{
+   double result = igamma_inv_evaluate2_dcd(igamma_big_data);
+   result += igamma_inv_evaluate2_dcd(igamma_int_data);
+   result += igamma_inv_evaluate2_dcd(igamma_med_data);
+   result += igamma_inv_evaluate2_dcd(igamma_small_data);
+
+   consume_result(result);
+   set_call_count(
+      2 * (sizeof(igamma_big_data) 
+      + sizeof(igamma_int_data) 
+      + sizeof(igamma_med_data)
+      + sizeof(igamma_small_data)) / sizeof(igamma_big_data[0]));
+}
+
+#endif
+
+double x = 0.165048161769598689119220580323599278926849365234375e-11;
+double y = 0.165048164480104120332981665342231281101703643798828125e-13;
+
+BOOST_MATH_PERFORMANCE_TEST(igamma_scrap, "igamma_scrap")
+{
+   //double result = boost::math::gamma_q(x, y);
+   double result = igamma_evaluate2(igamma_small_data);
+   //result += dcd::gamma_p(x, y);
+
+   consume_result(result);
+   set_call_count(1);
+}
+
+BOOST_MATH_PERFORMANCE_TEST(igamma_scrap, "igamma_scrap-dcd")
+{
+   //double result = dcd::gamma_q(x, y);
+   double result = igamma_evaluate2_dcd(igamma_small_data);
+
+   consume_result(result);
+   set_call_count(1);
+}
+
diff --git a/test/test_gamma_hooks.hpp b/test/test_gamma_hooks.hpp
index f33fd9a2f..b07c5b83a 100644
--- a/test/test_gamma_hooks.hpp
+++ b/test/test_gamma_hooks.hpp
@@ -162,12 +162,117 @@ inline long double gamma_p(long double x, long double y)
 }
 #endif
 
+#ifdef TEST_DCDFLIB
+#define TEST_OTHER
+#include <dcdflib.h>
+
+namespace other{
+float tgamma(float z)
+{
+   double v = z;
+   return (float)gamma_x(&v);
+}
+double tgamma(double z)
+{
+   return gamma_x(&z);
+}
+long double tgamma(long double z)
+{
+   double v = z;
+   return gamma_x(&v);
+}
+float lgamma(float z)
+{
+   double v = z;
+   return (float)gamma_log(&v);
+}
+double lgamma(double z)
+{
+   double v = z;
+   return gamma_log(&v);
+}
+long double lgamma(long double z)
+{
+   double v = z;
+   return gamma_log(&v);
+}
+inline double gamma_q(double x, double y)
+{ 
+   double ans, qans;
+   int i = 0;
+   gamma_inc (&x, &y, &ans, &qans, &i); 
+   return qans;
+}
+inline float gamma_q(float x, float y)
+{ 
+   return (float)gamma_q((double)x, (double)y); 
+}
+inline long double gamma_q(long double x, long double y)
+{ 
+   return gamma_q((double)x, (double)y); 
+}
+inline double gamma_p(double x, double y)
+{ 
+   double ans, qans;
+   int i = 0;
+   gamma_inc (&x, &y, &ans, &qans, &i); 
+   return ans;
+}
+inline float gamma_p(float x, float y)
+{ 
+   return (float)gamma_p((double)x, (double)y); 
+}
+inline long double gamma_p(long double x, long double y)
+{ 
+   return gamma_p((double)x, (double)y); 
+}
+
+inline double gamma_q_inv(double x, double y)
+{ 
+   double ans, p, nul;
+   int i = 0;
+   p = 1 - y;
+   nul = 0;
+   gamma_inc_inv (&x, &ans, &nul, &p, &y, &i); 
+   return ans;
+}
+inline float gamma_q_inv(float x, float y)
+{ 
+   return (float)gamma_q_inv((double)x, (double)y); 
+}
+inline long double gamma_q_inv(long double x, long double y)
+{ 
+   return gamma_q_inv((double)x, (double)y); 
+}
+inline double gamma_p_inv(double x, double y)
+{ 
+   double ans, p, nul;
+   int i = 0;
+   p = 1 - y;
+   nul = 0;
+   gamma_inc_inv (&x, &ans, &nul, &y, &p, &i); 
+   return ans;
+}
+inline float gamma_p_inv(float x, float y)
+{ 
+   return (float)gamma_p_inv((double)x, (double)y); 
+}
+inline long double gamma_p_inv(long double x, long double y)
+{ 
+   return gamma_p_inv((double)x, (double)y); 
+}
+
+}
+#endif
+
 #ifdef TEST_OTHER
 namespace other{
    boost::math::concepts::real_concept tgamma(boost::math::concepts::real_concept){ return 0; }
    boost::math::concepts::real_concept lgamma(boost::math::concepts::real_concept){ return 0; }
    boost::math::concepts::real_concept gamma_q(boost::math::concepts::real_concept, boost::math::concepts::real_concept){ return 0; }
    boost::math::concepts::real_concept gamma_p(boost::math::concepts::real_concept, boost::math::concepts::real_concept){ return 0; }
+   boost::math::concepts::real_concept gamma_p_inv(boost::math::concepts::real_concept x, boost::math::concepts::real_concept y){ return 0; }
+   boost::math::concepts::real_concept gamma_q_inv(boost::math::concepts::real_concept x, boost::math::concepts::real_concept y){ return 0; }
 }
 #endif
 
diff --git a/test/test_igamma.cpp b/test/test_igamma.cpp
index a47e5d1d7..9646ef154 100644
--- a/test/test_igamma.cpp
+++ b/test/test_igamma.cpp
@@ -361,7 +361,7 @@ void do_test_gamma_2(const T& data, const char* type_name, const char* test_name
       bind_func(funcp, 0, 1),
       extract_result(3));
    handle_test_result(result, data[result.worst()], result.worst(), type_name, "boost::math::gamma_q", test_name);
-#if defined(TEST_CEPHES) || defined(TEST_GSL)
+#if defined(TEST_OTHER)
    //
    // test other gamma_q(T, T) against data:
    //
@@ -388,7 +388,7 @@ void do_test_gamma_2(const T& data, const char* type_name, const char* test_name
       bind_func(funcp, 0, 1),
       extract_result(5));
    handle_test_result(result, data[result.worst()], result.worst(), type_name, "boost::math::gamma_p", test_name);
-#if defined(TEST_CEPHES) || defined(TEST_GSL)
+#if defined(TEST_OTHER)
    //
    // test other gamma_p(T, T) against data:
    //
diff --git a/test/test_igamma_inv.cpp b/test/test_igamma_inv.cpp
index 812369479..79ece855b 100644
--- a/test/test_igamma_inv.cpp
+++ b/test/test_igamma_inv.cpp
@@ -315,6 +315,29 @@ void do_test_gamma_inv(const T& data, const char* type_name, const char* test_na
       bind_func(funcp, 0, 1),
       extract_result(3));
    handle_test_result(result, data[result.worst()], result.worst(), type_name, "boost::math::gamma_q_inv", test_name);
+#ifdef TEST_OTHER
+   if(boost::is_floating_point<value_type>::value)
+   {
+      funcp = other::gamma_p_inv;
+      //
+      // test gamma_p_inv(T, T) against data:
+      //
+      result = boost::math::tools::test(
+         data,
+         bind_func(funcp, 0, 1),
+         extract_result(2));
+      handle_test_result(result, data[result.worst()], result.worst(), type_name, "other::gamma_p_inv", test_name);
+      //
+      // test gamma_q_inv(T, T) against data:
+      //
+      funcp = other::gamma_q_inv;
+      result = boost::math::tools::test(
+         data,
+         bind_func(funcp, 0, 1),
+         extract_result(3));
+      handle_test_result(result, data[result.worst()], result.worst(), type_name, "other::gamma_q_inv", test_name);
+   }
+#endif
 }
 
 template <class T>