diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp
index ed6383991..2cadf7f66 100644
--- a/include/boost/math/special_functions/gamma.hpp
+++ b/include/boost/math/special_functions/gamma.hpp
@@ -1207,6 +1207,68 @@ T tgamma_delta_ratio_imp(T z, T delta, const Policy& pol)
    return tgamma_delta_ratio_imp_lanczos(z, delta, pol, lanczos_type());
 }
 
+template <class T, class Policy>
+T tgamma_ratio_imp(T x, T y, const Policy& pol)
+{
+   BOOST_MATH_STD_USING
+
+   if(x <= 0)
+      policies::raise_domain_error<T>("boost::math::tgamma_ratio<%1%>(%1%, %1%)", "Gamma function ratios only implemented for positive arguments (got a=%1%).", x, pol);
+   if(y <= 0)
+      policies::raise_domain_error<T>("boost::math::tgamma_ratio<%1%>(%1%, %1%)", "Gamma function ratios only implemented for positive arguments (got b=%1%).", y, pol);
+
+   if((x < max_factorial<T>::value) && (y < max_factorial<T>::value))
+   {
+      // Rather than subtracting values, lets just call the gamma functions directly:
+      return boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
+   }
+   T prefix = 1;
+   if(x < 1)
+   {
+      if(y < 2 * max_factorial<T>::value)
+      {
+         // We need to sidestep on x as well, otherwise we'll underflow
+         // before we get to factor in the prefix term:
+         prefix /= x;
+         x += 1;
+         while(y >=  max_factorial<T>::value)
+         {
+            y -= 1;
+            prefix /= y;
+         }
+         return prefix * boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
+      }
+      //
+      // result is almost certainly going to underflow to zero, try logs just in case:
+      //
+      return exp(boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol));
+   }
+   if(y < 1)
+   {
+      if(x < 2 * max_factorial<T>::value)
+      {
+         // We need to sidestep on y as well, otherwise we'll overflow
+         // before we get to factor in the prefix term:
+         prefix *= y;
+         y += 1;
+         while(x >= max_factorial<T>::value)
+         {
+            x -= 1;
+            prefix *= x;
+         }
+         return prefix * boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
+      }
+      //
+      // Result will almost certainly overflow, try logs just in case:
+      //
+      return exp(boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol));
+   }
+   //
+   // Regular case, x and y both large and similar in magnitude:
+   //
+   return boost::math::tgamma_delta_ratio(x, y - x, pol);
+}
+
 template <class T, class Policy>
 T gamma_p_derivative_imp(T a, T x, const Policy& pol)
 {
@@ -1609,7 +1671,7 @@ inline typename tools::promote_args<T1, T2>::type
       policies::discrete_quantile<>,
       policies::assert_undefined<> >::type forwarding_policy;
 
-   return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::tgamma_delta_ratio_imp(static_cast<value_type>(a), static_cast<value_type>(static_cast<value_type>(b) - static_cast<value_type>(a)), forwarding_policy()), "boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)");
+   return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::tgamma_ratio_imp(static_cast<value_type>(a), static_cast<value_type>(b), forwarding_policy()), "boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)");
 }
 template <class T1, class T2>
 inline typename tools::promote_args<T1, T2>::type 
diff --git a/test/test_tgamma_ratio.hpp b/test/test_tgamma_ratio.hpp
index 48403f143..d8d9101a0 100644
--- a/test/test_tgamma_ratio.hpp
+++ b/test/test_tgamma_ratio.hpp
@@ -116,5 +116,18 @@ void test_tgamma_ratio(T, const char* name)
 
    do_test_tgamma_ratio<T>(tgamma_ratio_data, name, "tgamma ratios");
 
+   //
+   // A few special spot tests:
+   //
+   BOOST_MATH_STD_USING
+   T tol = boost::math::tools::epsilon<T>() * 20;
+   if(std::numeric_limits<T>::max_exponent > 200)
+   {
+      BOOST_CHECK_CLOSE_FRACTION(boost::math::tgamma_ratio(ldexp(T(1), -500), T(180.25)), 8.0113754557649679470816892372669519037339812035512e-178L, tol);
+      BOOST_CHECK_CLOSE_FRACTION(boost::math::tgamma_ratio(ldexp(T(1), -525), T(192.25)), 1.5966560279353205461166489184101261541784867035063e-197L, tol);
+      BOOST_CHECK_CLOSE_FRACTION(boost::math::tgamma_ratio(T(182.25), ldexp(T(1), -500)), 4.077990437521002194346763299159975185747917450788e+181L, tol);
+      BOOST_CHECK_CLOSE_FRACTION(boost::math::tgamma_ratio(T(193.25), ldexp(T(1), -525)), 1.2040790040958522422697601672703926839178050326148e+199L, tol);
+      BOOST_CHECK_CLOSE_FRACTION(boost::math::tgamma_ratio(T(193.25), T(194.75)), 0.00037151765099653237632823607820104961270831942138159L, tol);
+   }
 }