diff --git a/include/boost/math/special_functions/math_fwd.hpp b/include/boost/math/special_functions/math_fwd.hpp
index 0e0815790..ca16fc8d9 100644
--- a/include/boost/math/special_functions/math_fwd.hpp
+++ b/include/boost/math/special_functions/math_fwd.hpp
@@ -561,6 +561,12 @@ namespace boost
    template <class RT1, class RT2, class Policy>
    BOOST_MATH_GPU_ENABLED tools::promote_args_t<RT1, RT2> gamma_q(RT1 a, RT2 z, const Policy&);
 
+   template <class RT1, class RT2>
+   BOOST_MATH_GPU_ENABLED tools::promote_args_t<RT1, RT2> lgamma_q(RT1 a, RT2 z);
+
+   template <class RT1, class RT2, class Policy>
+   BOOST_MATH_GPU_ENABLED tools::promote_args_t<RT1, RT2> lgamma_q(RT1 a, RT2 z, const Policy&);
+
    template <class RT1, class RT2>
    BOOST_MATH_GPU_ENABLED tools::promote_args_t<RT1, RT2> gamma_p(RT1 a, RT2 z);
 
@@ -1516,6 +1522,9 @@ namespace boost
 \
    template <class RT1, class RT2>\
    BOOST_MATH_GPU_ENABLED inline boost::math::tools::promote_args_t<RT1, RT2> gamma_q(RT1 a, RT2 z){ return boost::math::gamma_q(a, z, Policy()); }\
+\
+   template <class RT1, class RT2>\
+   BOOST_MATH_GPU_ENABLED inline boost::math::tools::promote_args_t<RT1, RT2> lgamma_q(RT1 a, RT2 z){ return boost::math::lgamma_q(a, z, Policy()); }\
 \
    template <class RT1, class RT2>\
    BOOST_MATH_GPU_ENABLED inline boost::math::tools::promote_args_t<RT1, RT2> gamma_p(RT1 a, RT2 z){ return boost::math::gamma_p(a, z, Policy()); }\
diff --git a/test/compile_test/instantiate.hpp b/test/compile_test/instantiate.hpp
index 26fa53d7b..1d186aff7 100644
--- a/test/compile_test/instantiate.hpp
+++ b/test/compile_test/instantiate.hpp
@@ -263,6 +263,7 @@ void instantiate(RealType)
    boost::math::tgamma_lower(v1, v2);
    boost::math::gamma_p(v1, v2);
    boost::math::gamma_q(v1, v2);
+   boost::math::lgamma_q(v1, v2);
    boost::math::gamma_p_inv(v1, v2);
    boost::math::gamma_q_inv(v1, v2);
    boost::math::gamma_p_inva(v1, v2);
@@ -542,6 +543,7 @@ void instantiate(RealType)
    boost::math::tgamma_lower(v1 * 1, v2 - 0);
    boost::math::gamma_p(v1 * 1, v2 + 0);
    boost::math::gamma_q(v1 * 1, v2 + 0);
+   boost::math::lgamma_q(v1 * 1, v2 + 0);
    boost::math::gamma_p_inv(v1 * 1, v2 + 0);
    boost::math::gamma_q_inv(v1 * 1, v2 + 0);
    boost::math::gamma_p_inva(v1 * 1, v2 + 0);
@@ -793,6 +795,7 @@ void instantiate(RealType)
    boost::math::tgamma_lower(v1, v2, pol);
    boost::math::gamma_p(v1, v2, pol);
    boost::math::gamma_q(v1, v2, pol);
+   boost::math::lgamma_q(v1, v2, pol);
    boost::math::gamma_p_inv(v1, v2, pol);
    boost::math::gamma_q_inv(v1, v2, pol);
    boost::math::gamma_p_inva(v1, v2, pol);
@@ -1070,6 +1073,7 @@ void instantiate(RealType)
    test::tgamma_lower(v1, v2);
    test::gamma_p(v1, v2);
    test::gamma_q(v1, v2);
+   test::lgamma_q(v1, v2);
    test::gamma_p_inv(v1, v2);
    test::gamma_q_inv(v1, v2);
    test::gamma_p_inva(v1, v2);
@@ -1351,6 +1355,7 @@ void instantiate_mixed(RealType)
    boost::math::gamma_p(i, s);
    boost::math::gamma_p(fr, lr);
    boost::math::gamma_q(i, s);
+   boost::math::lgamma_q(i, s);
    boost::math::gamma_q(fr, lr);
    boost::math::gamma_p_inv(i, fr);
    boost::math::gamma_q_inv(s, fr);
@@ -1566,6 +1571,7 @@ void instantiate_mixed(RealType)
    boost::math::gamma_p(i, s, pol);
    boost::math::gamma_p(fr, lr, pol);
    boost::math::gamma_q(i, s, pol);
+   boost::math::lgamma_q(i, s, pol);
    boost::math::gamma_q(fr, lr, pol);
    boost::math::gamma_p_inv(i, fr, pol);
    boost::math::gamma_q_inv(s, fr, pol);
@@ -1777,7 +1783,9 @@ void instantiate_mixed(RealType)
    test::gamma_p(i, s);
    test::gamma_p(fr, lr);
    test::gamma_q(i, s);
+   test::lgamma_q(i, s);
    test::gamma_q(fr, lr);
+   test::lgamma_q(fr, lr);
    test::gamma_p_inv(i, fr);
    test::gamma_q_inv(s, fr);
    test::gamma_p_inva(i, lr);
diff --git a/test/compile_test/sf_gamma_incl_test.cpp b/test/compile_test/sf_gamma_incl_test.cpp
index f7ad18e9c..74ab85b2d 100644
--- a/test/compile_test/sf_gamma_incl_test.cpp
+++ b/test/compile_test/sf_gamma_incl_test.cpp
@@ -39,6 +39,12 @@ void compile_and_link_test()
    check_result<long double>(boost::math::gamma_q<long double>(l, l));
 #endif
 
+   check_result<float>(boost::math::lgamma_q<float>(f, f));
+   check_result<double>(boost::math::lgamma_q<double>(d, d));
+#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
+   check_result<long double>(boost::math::lgamma_q<long double>(l, l));
+#endif
+
    check_result<float>(boost::math::gamma_p_inv<float>(f, f));
    check_result<double>(boost::math::gamma_p_inv<double>(d, d));
 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
diff --git a/test/test_igamma.hpp b/test/test_igamma.hpp
index c2b6820bc..85b86f421 100644
--- a/test/test_igamma.hpp
+++ b/test/test_igamma.hpp
@@ -280,7 +280,7 @@ void test_spots(T, const char* name = nullptr)
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(800), static_cast<T>(1000.25)), static_cast<T>(-24.43514173634027477093666725985191846106997808357863808910970142L), tolerance * (boost::math::tools::digits<T>() > 54 ? 20 : 1));
    // Once we get large a,x then error start to accumulate no matter what we do:
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(1200), static_cast<T>(1249.75)), static_cast<T>(-2.565496161584661216769813239648606145255794643472303927896044375L), tolerance * (std::is_floating_point<T>::value ? 1 : 4));
-   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(1200), static_cast<T>(1250.25)), static_cast<T>(-2.591934862117586205519309712218581885256650074210410262843591453L), tolerance * ((std::numeric_limits<T>::max_digits10 >= 36) ? 500 : (std::is_same<T, float>::value ? 1 : 50)));
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(1200), static_cast<T>(1250.25)), static_cast<T>(-2.591934862117586205519309712218581885256650074210410262843591453L), tolerance * ((std::numeric_limits<T>::max_digits10 >= 36) ? 750 : (std::is_same<T, float>::value ? 1 : 50)));
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(2200), static_cast<T>(2249.75)), static_cast<T>(-1.933779894897391651410597618307863427927461116308937004149240320L), tolerance * (std::is_floating_point<T>::value ? 1 : 10));
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(2200), static_cast<T>(2250.25)), static_cast<T>(-1.950346484067948344620463026377077515919992808640737320057812268L), tolerance * (std::is_same<T, float>::value ? 1 : (std::is_floating_point<T>::value ? 100 : 200)));
    //