Implement special function log incomplete gamma function

include/boost/math/special_functions/gamma.hpp

@@ -1772,6 +1772,51 @@ BOOST_MATH_GPU_ENABLED T gamma_incomplete_imp(T a, T x, bool normalised, bool in
    return gamma_incomplete_imp_final(T(a), T(x), normalised, invert, pol, p_derivative);
 }
 
+// Calculate log of incomplete gamma function
+template <class T, class Policy>
+T lgamma_incomplete_imp_final(T a, T x, const Policy& pol)
+{
+   using namespace boost::math;  // temporary until we're in the right namespace
+
+   BOOST_MATH_STD_USING_CORE
+
+   if (((x > 1000) && ((a < x) || (fabs(a - 50) / x < 1))) || ((x > tools::log_max_value<T>() - 10) && (x > a)))
+   {
+      //
+      // Take the logarithmic version of the asymtotic expansion:
+      //
+      return log(detail::incomplete_tgamma_large_x(a, x, pol)) + a * log(x) - x - lgamma(a, pol) - log(x);
+   }
+   //
+   // Can't do better than taking the log of Q, but...
+   //
+   // Figure out whether we need P or Q, since if we calculate Q and it's too close to unity
+   // we will loose 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 log1p(-gamma_p(a, x, pol), pol);
+   return log(gamma_q(a, x, pol));
+}
+
+template <class T, class Policy> 
+T lgamma_incomplete_imp(T a, T x, const Policy& pol){
+   constexpr auto function = "boost::math::lgamma_q<%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);
+
+   // If input is valid proceed as normal
+   return lgamma_incomplete_imp_final(T(a), T(x), pol);
+}
 //
 // Ratios of two gamma functions:
 //
@@ -2390,6 +2435,29 @@ BOOST_MATH_GPU_ENABLED inline tools::promote_args_t<T1, T2>
 {
    return gamma_q(a, z, policies::policy<>());
 }
+
+template <class T1, class T2, class Policy>
+inline tools::promote_args_t<T1, T2> lgamma_q(T1 a, T2 z, const Policy& /* pol */)
+{
+   typedef tools::promote_args_t<T1, T2> result_type;
+   typedef typename policies::evaluation<result_type, Policy>::type value_type;
+   typedef typename policies::normalise<
+      Policy,
+      policies::promote_float<false>,
+      policies::promote_double<false>,
+      policies::discrete_quantile<>,
+      policies::assert_undefined<> >::type forwarding_policy;
+
+   return policies::checked_narrowing_cast<result_type, forwarding_policy>(
+      detail::lgamma_incomplete_imp(static_cast<value_type>(a),
+      static_cast<value_type>(z), forwarding_policy()), "lgamma_q<%1%>(%1%, %1%)");
+}
+
+template <class T1, class T2>
+inline tools::promote_args_t<T1, T2> lgamma_q(T1 a, T2 z)
+{
+   return lgamma_q(a, z, policies::policy<>());
+}
 //
 // Regularised lower incomplete gamma:
 //