diff --git a/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp b/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp
index 3adc86f74f..58c0d6b1ea 100644
--- a/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp
+++ b/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp
@@ -23,6 +23,7 @@
      unsigned set_crossover_locations(const Seq& aj, const Seq& bj, const Real& z, unsigned int* crossover_locations)
      {
         BOOST_MATH_STD_USING
+        using nothrow_policy = typename boost::math::policies::normalise<boost::math::policies::policy<>, boost::math::policies::rounding_error<boost::math::policies::ignore_error>>::type;
         unsigned N_terms = 0;
 
         if(aj.size() == 1 && bj.size() == 1)
@@ -55,13 +56,13 @@
               Real t = (-sqrt(sq) - b + z) / 2;
               if (t >= 0)
               {
-                 crossover_locations[N_terms] = itrunc(t);
+                 crossover_locations[N_terms] = itrunc(t, nothrow_policy());
                  ++N_terms;
               }
               t = (sqrt(sq) - b + z) / 2;
               if (t >= 0)
               {
-                 crossover_locations[N_terms] = itrunc(t);
+                 crossover_locations[N_terms] = itrunc(t, nothrow_policy());
                  ++N_terms;
               }
            }
@@ -71,13 +72,13 @@
               Real t = (-sqrt(sq) - b - z) / 2;
               if (t >= 0)
               {
-                 crossover_locations[N_terms] = itrunc(t);
+                 crossover_locations[N_terms] = itrunc(t, nothrow_policy());
                  ++N_terms;
               }
               t = (sqrt(sq) - b - z) / 2;
               if (t >= 0)
               {
-                 crossover_locations[N_terms] = itrunc(t);
+                 crossover_locations[N_terms] = itrunc(t, nothrow_policy());
                  ++N_terms;
               }
            }
@@ -110,7 +111,7 @@
            unsigned n = 0;
            for (auto bi = bj.begin(); bi != bj.end(); ++bi, ++n)
            {
-              crossover_locations[n] = *bi >= 0 ? 0 : itrunc(-*bi) + 1;
+              crossover_locations[n] = *bi >= 0 ? 0 : itrunc(-*bi, nothrow_policy()) + 1;
            }
            std::sort(crossover_locations, crossover_locations + bj.size(), std::less<Real>());
            N_terms = (unsigned)bj.size();
@@ -122,6 +123,7 @@
      std::pair<Real, Real> hypergeometric_pFq_checked_series_impl(const Seq& aj, const Seq& bj, const Real& z, const Policy& pol, const Terminal& termination, long long& log_scale)
      {
         BOOST_MATH_STD_USING
+        using nothrow_policy = typename boost::math::policies::normalise<Policy, boost::math::policies::rounding_error<boost::math::policies::ignore_error>>::type;
         Real result = 1;
         Real abs_result = 1;
         Real term = 1;
@@ -129,8 +131,15 @@
         Real tol = boost::math::policies::get_epsilon<Real, Policy>();
         std::uintmax_t k = 0;
         Real upper_limit(sqrt(boost::math::tools::max_value<Real>())), diff;
+        if ((tools::max_value<Real>() / fabs(z) < upper_limit))
+            upper_limit = tools::max_value<Real>() / fabs(z);
+        for (auto pa = aj.begin(); pa != aj.end(); ++pa)
+        {
+            if (tools::max_value<Real>() / fabs(*pa) < upper_limit)
+                upper_limit = tools::max_value<Real>() / fabs(*pa);
+        }
         Real lower_limit(1 / upper_limit);
-        long long log_scaling_factor = lltrunc(boost::math::tools::log_max_value<Real>()) - 2;
+        long long log_scaling_factor = lltrunc(boost::math::tools::log_max_value<Real>(), nothrow_policy()) - 2;
         Real scaling_factor = exp(Real(log_scaling_factor));
         Real term_m1;
         long long local_scaling = 0;
@@ -234,7 +243,7 @@
               log_scale += log_scaling_factor;
               local_scaling += log_scaling_factor;
            }
-           if (fabs(abs_result) < lower_limit)
+           else if ((fabs(abs_result) < lower_limit) && (fabs(abs_result) * scaling_factor < upper_limit))
            {
               abs_result *= scaling_factor;
               result *= scaling_factor;
@@ -321,7 +330,7 @@
                  else
                  {
                     int ls = 1;
-                    Real p = log_pochhammer(*ai, s, pol, &ls);
+                    Real p = log_pochhammer(static_cast<Real>(*ai), s, pol, &ls);
                     s1 *= ls;
                     term += p;
                     loop_error_scale = (std::max)(p, loop_error_scale);
@@ -334,7 +343,7 @@
               for (auto bi = bj.begin(); bi != bj.end(); ++bi)
               {
                  int ls = 1;
-                 Real p = log_pochhammer(*bi, s, pol, &ls);
+                 Real p = log_pochhammer(static_cast<Real>(*bi), s, pol, &ls);
                  s2 *= ls;
                  term -= p;
                  loop_error_scale = (std::max)(p, loop_error_scale);
@@ -371,13 +380,13 @@
               if (term <= tools::log_min_value<Real>())
               {
                  // rescale if we can:
-                 long long scale = lltrunc(floor(term - tools::log_min_value<Real>()) - 2);
+                 long long scale = lltrunc(floor(term - tools::log_min_value<Real>()) - 2, nothrow_policy());
                  term -= scale;
                  loop_scale += scale;
               }
                if (term > 10)
                {
-                  int scale = itrunc(floor(term));
+                  long long scale = lltrunc(floor(term), nothrow_policy());
                   term -= scale;
                   loop_scale += scale;
                }
@@ -402,7 +411,7 @@
                      term /= scaling_factor;
                      loop_scale += log_scaling_factor;
                   }
-                  if (fabs(loop_result) < lower_limit)
+                  else if ((fabs(loop_result) < lower_limit) && (fabs(loop_result) * scaling_factor < upper_limit))
                   {
                      loop_result *= scaling_factor;
                      loop_abs_result *= scaling_factor;
diff --git a/include/boost/math/special_functions/hypergeometric_pFq.hpp b/include/boost/math/special_functions/hypergeometric_pFq.hpp
index 070b852dd2..96f68e510e 100644
--- a/include/boost/math/special_functions/hypergeometric_pFq.hpp
+++ b/include/boost/math/special_functions/hypergeometric_pFq.hpp
@@ -58,12 +58,22 @@ namespace boost {
          BOOST_MATH_STD_USING
 
          long long scale = 0;
+         static const char* function = "boost::math::hypergeometric_pFq<%1%>(%1%,%1%,%1%)";
          std::pair<value_type, value_type> r = boost::math::detail::hypergeometric_pFq_checked_series_impl(aj, bj, value_type(z), pol, boost::math::detail::iteration_terminator(boost::math::policies::get_max_series_iterations<forwarding_policy>()), scale);
-         r.first *= exp(Real(scale));
-         r.second *= exp(Real(scale));
+         //
+         // Overflow check:
+         //
+         if (static_cast<Real>(scale) > tools::log_max_value<Real>())
+             return (r.first < 0 ? -1 : 1) * policies::raise_overflow_error<Real, Policy>(function, nullptr, pol);
+         Real mul = exp(Real(scale));
+         if(fabs(r.first) > 1)
+             if(tools::max_value<Real>() / fabs(r.first) < mul)
+                 return (r.first < 0 ? -1 : 1) * policies::raise_overflow_error<Real, Policy>(function, nullptr, pol);
+         r.first *= mul;
+         r.second *= mul;
          if (p_abs_error)
             *p_abs_error = static_cast<Real>(r.second) * boost::math::tools::epsilon<Real>();
-         return policies::checked_narrowing_cast<result_type, Policy>(r.first, "boost::math::hypergeometric_pFq<%1%>(%1%,%1%,%1%)");
+         return policies::checked_narrowing_cast<result_type, Policy>(r.first, function);
       }
 
       template <class Seq, class Real>
diff --git a/test/test_pFq.hpp b/test/test_pFq.hpp
index 0e7f3a0c54..9cfbac30ba 100644
--- a/test/test_pFq.hpp
+++ b/test/test_pFq.hpp
@@ -298,6 +298,31 @@ void test_spots_2F2(T, const char*)
    }
 }
 
+template <class T>
+void test_special_cases(T, const char*)
+{
+    typedef boost::math::policies::policy<boost::math::policies::overflow_error<boost::math::policies::ignore_error>> nothrow_policy;
+    typedef boost::math::policies::policy<boost::math::policies::overflow_error<boost::math::policies::throw_on_error>> throw_policy;
+    T error_rate;
+    BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({ 2 }, { 3 }, static_cast<T>(1e10), &error_rate, throw_policy()), std::overflow_error);
+    BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({ static_cast<T>(2) }, { static_cast<T>(3) }, static_cast<T>(1e10), &error_rate, throw_policy()), std::overflow_error);
+    BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({ 2 }, { 3 }, boost::math::tools::max_value<T>(), &error_rate, throw_policy()), std::overflow_error);
+    BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({ static_cast<T>(2) }, { static_cast<T>(3) }, boost::math::tools::max_value<T>(), &error_rate, throw_policy()), std::overflow_error);
+    BOOST_MATH_IF_CONSTEXPR(std::numeric_limits<T>::is_specialized)
+    {
+        // Can't figure out how to make this work with real_concept (we get an evaluation_error):
+        BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({ boost::math::tools::max_value<T>() }, { static_cast<T>(3) }, static_cast<T>(0.5f), &error_rate, throw_policy()), std::overflow_error);
+    }
+    BOOST_MATH_IF_CONSTEXPR(std::numeric_limits<T>::has_infinity)
+    {
+        BOOST_CHECK_EQUAL(boost::math::hypergeometric_pFq({ 2 }, { 3 }, static_cast<T>(1e10), &error_rate, nothrow_policy()), std::numeric_limits<T>::infinity());
+        BOOST_CHECK_EQUAL(boost::math::hypergeometric_pFq({ static_cast<T>(2) }, { static_cast<T>(3) }, static_cast<T>(1e10), &error_rate, nothrow_policy()), std::numeric_limits<T>::infinity());
+        BOOST_CHECK_EQUAL(boost::math::hypergeometric_pFq({ 2 }, { 3 }, boost::math::tools::max_value<T>(), &error_rate, nothrow_policy()), std::numeric_limits<T>::infinity());
+        BOOST_CHECK_EQUAL(boost::math::hypergeometric_pFq({ static_cast<T>(2) }, { static_cast<T>(3) }, boost::math::tools::max_value<T>(), &error_rate, nothrow_policy()), std::numeric_limits<T>::infinity());
+        BOOST_CHECK_EQUAL(boost::math::hypergeometric_pFq({ boost::math::tools::max_value<T>() }, { static_cast<T>(3) }, static_cast<T>(0.5f), &error_rate, nothrow_policy()), std::numeric_limits<T>::infinity());
+    }
+}
+
 template <class T>
 void test_spots(T z, const char* type_name)
 {
@@ -309,4 +334,5 @@ void test_spots(T z, const char* type_name)
    test_spots_1F2(z, type_name);
    test_spots_2F2(z, type_name);
    test_spots_2F1(z, type_name);
+   test_special_cases(z, type_name);
 }