diff --git a/include/boost/math/special_functions/beta.hpp b/include/boost/math/special_functions/beta.hpp
index fb78aefdd..ae6ee30ff 100644
--- a/include/boost/math/special_functions/beta.hpp
+++ b/include/boost/math/special_functions/beta.hpp
@@ -12,6 +12,7 @@
 #include <boost/math/special_functions/log1p.hpp>
 #include <boost/math/special_functions/expm1.hpp>
 #include <boost/math/tools/roots.hpp>
+#include <cmath>
 
 namespace boost{ namespace math{
 
@@ -548,6 +549,34 @@ T rising_factorial_ratio(T a, T b, int k)
 //
 // Routine for a > 15, b < 1
 //
+// Begin by figuring out how large our table of Pn's should be,
+// quoted accuracies are "guestimates" based on empiracal observation.
+// Note that the table size should never exceed the size of our 
+// tables of factorials.
+//
+template <class T>
+struct Pn_size
+{
+   // This is likely to be enough for ~100 digit accuracy
+   // but it's hard to quantify exactly:
+   BOOST_STATIC_CONSTANT(unsigned, value = 100);
+};
+template <>
+struct Pn_size<float>
+{
+   BOOST_STATIC_CONSTANT(unsigned, value = 15); // ~8-15 digit accuracy
+};
+template <>
+struct Pn_size<double>
+{
+   BOOST_STATIC_CONSTANT(unsigned, value = 30); // 16-20 digit accuracy
+};
+template <>
+struct Pn_size<long double>
+{
+   BOOST_STATIC_CONSTANT(unsigned, value = 50); // ~35-50 digit accuracy
+};
+
 template <class T, class L>
 T beta_small_b_large_a_series(T a, T b, T x, T y, T s0, T mult, const L& l, bool normalised)
 {
@@ -585,7 +614,7 @@ T beta_small_b_large_a_series(T a, T b, T x, T y, T s0, T mult, const L& l, bool
    // recursively, and requires a full history of all the previous values
    // so no choice but to declare a big table and hope it's big enough...
    //
-   T p[30] = { 1 };  // see 9.3.
+   T p[ ::boost::math::detail::Pn_size<T>::value ] = { 1 };  // see 9.3.
    //
    // Now an initial value for J, see 9.6:
    //
@@ -602,8 +631,19 @@ T beta_small_b_large_a_series(T a, T b, T x, T y, T s0, T mult, const L& l, bool
    T t4 = 4 * t * t;
    T b2n = b; 
 
-   for(unsigned n = 1; n < 30; ++n)
+   for(unsigned n = 1; n < sizeof(p)/sizeof(p[0]); ++n)
    {
+      /*
+      // debugging code, enable this if you want to determine whether
+      // the table of Pn's is large enough...
+      //
+      static int max_count = 2;
+      if(n > max_count)
+      {
+         max_count = n;
+         std::cerr << "Max iterations in BGRAT was " << n << std::endl;
+      }
+      */
       //
       // begin by evaluating the next Pn from Eq 9.4:
       //
@@ -630,8 +670,16 @@ T beta_small_b_large_a_series(T a, T b, T x, T y, T s0, T mult, const L& l, bool
       //
       T r = prefix * p[n] * j;
       sum += r;
-      if(fabs(r) < tools::epsilon(r) * sum)
-         break;
+      if(r > 1)
+      {
+         if(fabs(r) < fabs(tools::epsilon(r) * sum))
+            break;
+      }
+      else
+      {
+         if(fabs(r / tools::epsilon(r)) < fabs(sum))
+            break;
+      }
    }
    return sum;
 }
diff --git a/include/boost/math/special_functions/factorials.hpp b/include/boost/math/special_functions/factorials.hpp
index ebfc76dfd..bb6742a21 100644
--- a/include/boost/math/special_functions/factorials.hpp
+++ b/include/boost/math/special_functions/factorials.hpp
@@ -8,6 +8,7 @@
 
 #include <boost/math/special_functions/gamma.hpp>
 #include <boost/array.hpp>
+#include <cmath>
 
 namespace boost{ namespace math{
 
@@ -274,14 +275,14 @@ inline unsigned max_factorial()
 template <class T>
 T factorial(unsigned i)
 {
-   using namspace std;
+   using namespace std;
 
    if(i <= max_factorial<T>())
       return unchecked_factorial<T>(i);
    T result = boost::math::tgamma(i+1);
    if(result > tools::max_value(result))
       return result; // overflowed value, tgamma will have signalled the error already
-   return floor(result + 0.5)
+   return floor(result + 0.5);
 }
 
 template<>
diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp
index 53824e8c3..bc9fd3fc0 100644
--- a/include/boost/math/special_functions/gamma.hpp
+++ b/include/boost/math/special_functions/gamma.hpp
@@ -781,6 +781,8 @@ T gamma_P_imp(T a, T z, const L& l)
 template <class T, class L>
 T tgamma_delta_ratio_imp(T z, T delta, const L&)
 {
+   using namespace std;
+
    if((z <= 0) || (z + delta <= 0))
       tools::domain_error<T>(BOOST_CURRENT_FUNCTION, "Gamma function ratios only implemented for positive arguments.");
 
@@ -802,6 +804,8 @@ T tgamma_delta_ratio_imp(T z, T delta, const L&)
 template <class T>
 T tgamma_delta_ratio_imp(T z, T delta, const lanczos::undefined_lanczos&)
 {
+   using namespace std;
+
    if((z <= 0) || (z + delta <= 0))
       tools::domain_error<T>(BOOST_CURRENT_FUNCTION, "Gamma function ratios only implemented for positive arguments.");
 
diff --git a/include/boost/math/tools/test.hpp b/include/boost/math/tools/test.hpp
index 8cdc71dc5..3b97bd974 100644
--- a/include/boost/math/tools/test.hpp
+++ b/include/boost/math/tools/test.hpp
@@ -50,6 +50,9 @@ T relative_error(T a, T b)
    T min_val = (std::max)(
       tools::min_value(a), 
       static_cast<T>((std::numeric_limits<double>::min)()));
+   T max_val = (std::min)(
+      tools::max_value(a), 
+      static_cast<T>((std::numeric_limits<double>::max)()));
 #else
    T min_val = tools::min_value(a);
    T max_val = tools::max_value(a);