Implemented large a series expansion

Author: JacobHass8

include/boost/math/special_functions/gamma.hpp

@@ -1290,6 +1290,33 @@ BOOST_MATH_GPU_ENABLED T incomplete_tgamma_large_x(const T& a, const T& x, const
    return result;
 }
 
+template <class T>
+struct incomplete_tgamma_lower_large_a_series
+{
+   typedef T result_type;
+   BOOST_MATH_GPU_ENABLED incomplete_tgamma_lower_large_a_series(const T& a, const T& x)
+      : a_poch(a + 1), z(x), term(1 / a) {}
+   BOOST_MATH_GPU_ENABLED T operator()()
+   {
+      T result = term;
+      term *= z / a_poch;
+      a_poch += 1;
+      return result;
+   }
+   T a_poch, z, term;
+};
+
+template <class T, class Policy>
+T incomplete_tgamma_lower_large_a(const T& a, const T&x, const Policy & pol)
+{
+   BOOST_MATH_STD_USING
+   incomplete_tgamma_lower_large_a_series<T> s(a, x);
+   boost::math::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
+   T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
+   boost::math::policies::check_series_iterations<T>("boost::math::tgamma_p<%1%>(%1%,%1%)", max_iter, pol);
+   return result;
+}
+
 
 //
 // Main incomplete gamma entry point, handles all four incomplete gamma's:
@@ -1813,6 +1840,45 @@ BOOST_MATH_GPU_ENABLED T lgamma_incomplete_imp(T a, T x, const Policy& pol)
    return log(gamma_q(a, x, pol));
 }
 
+// Calculate log of incomplete gamma function
+template <class T, class Policy>
+BOOST_MATH_GPU_ENABLED T lgamma_incomplete_lower_imp(T a, T x, const Policy& pol)
+{
+   using namespace boost::math;  // temporary until we're in the right namespace
+
+   BOOST_MATH_STD_USING_CORE
+
+   // Check for invalid inputs (a < 0 or x < 0)
+   constexpr auto function = "boost::math::lgamma_p<%1%>(%1%, %1%)";
+   if(a <= 0)
+      return policies::raise_domain_error<T>(function, "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a, pol);
+   if(x < 0)
+      return policies::raise_domain_error<T>(function, "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x, pol);
+
+   // Need to change this to be more appropriate! 
+   if (a > 100){
+      return log(detail::incomplete_tgamma_lower_large_a(a, x, pol)) + a * log(x) - x - lgamma(a, pol);
+   }
+
+   //
+   // Can't do better than taking the log of P, but...
+   //
+   // Figure out whether we need P or Q, since if we calculate P and it's too close to unity
+   // we will lose precision in the result, selection logic here is extracted from gamma_incomplete_imp_final:
+   //
+   bool need_p = false;
+   if ((x < 0.5) && (T(-0.4) / log(x) < a))
+      need_p = true;
+   else if ((x < 1.1) && (x >= 0.5) && (x * 0.75f < a))
+      need_p = true;
+   else if ((x < a) && (x >= 1.1))
+      need_p = true;
+
+   if (need_p)
+      return log(gamma_p(a, x, pol));
+   return log1p(-gamma_q(a, x, pol), pol);
+}
+
 //
 // Ratios of two gamma functions:
 //