Imported Upstream version 1.60.0

Change-Id: Ie709530d6d5841088ceaba025cbe175a4ef43050
Signed-off-by: DongHun Kwak <dh0128.kwak@samsung.com>

26 files changed, 688 insertions, 318 deletions

diff --git a/boost/math/special_functions/bernoulli.hpp b/boost/math/special_functions/bernoulli.hpp
index 5f2b7e1422..e5323fbd42 100644
--- a/boost/math/special_functions/bernoulli.hpp
+++ b/boost/math/special_functions/bernoulli.hpp
@@ -63,7 +63,7 @@ inline T bernoulli_b2n(const int i, const Policy &pol)
    if(i < 0)
       return policies::raise_domain_error<T>("boost::math::bernoulli_b2n<%1%>", "Index should be >= 0 but got %1%", T(i), pol);
 
-   T result = 0; // The = 0 is just to silence compiler warings :-(
+   T result = static_cast<T>(0); // The = 0 is just to silence compiler warnings :-(
    boost::math::detail::bernoulli_number_imp<T>(&result, static_cast<std::size_t>(i), 1u, pol, tag_type());
    return result;
 }
diff --git a/boost/math/special_functions/bessel.hpp b/boost/math/special_functions/bessel.hpp
index eab723d8ba..1b57e0a873 100644
--- a/boost/math/special_functions/bessel.hpp
+++ b/boost/math/special_functions/bessel.hpp
@@ -302,20 +302,9 @@ inline T cyl_neumann_imp(T v, T x, const bessel_maybe_int_tag&, const Policy& po
 
    if(floor(v) == v)
    {
-      if(asymptotic_bessel_large_x_limit(v, x))
-      {
-         T r = asymptotic_bessel_y_large_x_2(static_cast<T>(abs(v)), x);
-         if((v < 0) && (itrunc(v, pol) & 1))
-            r = -r;
-         BOOST_MATH_INSTRUMENT_VARIABLE(r);
-         return r;
-      }
-      else
-      {
-         T r = bessel_yn(itrunc(v, pol), x, pol);
-         BOOST_MATH_INSTRUMENT_VARIABLE(r);
-         return r;
-      }
+      T r = bessel_yn(itrunc(v, pol), x, pol);
+      BOOST_MATH_INSTRUMENT_VARIABLE(r);
+      return r;
    }
    T r = cyl_neumann_imp<T>(v, x, bessel_no_int_tag(), pol);
    BOOST_MATH_INSTRUMENT_VARIABLE(r);
@@ -325,20 +314,7 @@ inline T cyl_neumann_imp(T v, T x, const bessel_maybe_int_tag&, const Policy& po
 template <class T, class Policy>
 inline T cyl_neumann_imp(int v, T x, const bessel_int_tag&, const Policy& pol)
 {
-   BOOST_MATH_STD_USING
-
-   BOOST_MATH_INSTRUMENT_VARIABLE(v);
-   BOOST_MATH_INSTRUMENT_VARIABLE(x);
-
-   if(asymptotic_bessel_large_x_limit(T(v), x))
-   {
-      T r = asymptotic_bessel_y_large_x_2(static_cast<T>(abs(v)), x);
-      if((v < 0) && (v & 1))
-         r = -r;
-      return r;
-   }
-   else
-      return bessel_yn(v, x, pol);
+   return bessel_yn(v, x, pol);
 }
 
 template <class T, class Policy>
@@ -582,7 +558,7 @@ inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_i(
       policies::promote_double<false>, 
       policies::discrete_quantile<>,
       policies::assert_undefined<> >::type forwarding_policy;
-   return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_i_imp<value_type>(v, static_cast<value_type>(x), forwarding_policy()), "boost::math::cyl_bessel_i<%1%>(%1%,%1%)");
+   return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_i_imp<value_type>(static_cast<value_type>(v), static_cast<value_type>(x), forwarding_policy()), "boost::math::cyl_bessel_i<%1%>(%1%,%1%)");
 }
 
 template <class T1, class T2>
diff --git a/boost/math/special_functions/bessel_prime.hpp b/boost/math/special_functions/bessel_prime.hpp
index c5f2d58c2b..2a09b2046b 100644
--- a/boost/math/special_functions/bessel_prime.hpp
+++ b/boost/math/special_functions/bessel_prime.hpp
@@ -237,7 +237,7 @@ inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_j_
       policies::promote_double<false>,
       policies::discrete_quantile<>,
       policies::assert_undefined<> >::type forwarding_policy;
-   return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_j_prime_imp<tag_type, value_type>(v, static_cast<value_type>(x), forwarding_policy()), "boost::math::cyl_bessel_j_prime<%1%,%1%>(%1%,%1%)");
+   return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_j_prime_imp<tag_type, value_type>(static_cast<value_type>(v), static_cast<value_type>(x), forwarding_policy()), "boost::math::cyl_bessel_j_prime<%1%,%1%>(%1%,%1%)");
 }
 
 template <class T1, class T2>
@@ -279,7 +279,7 @@ inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_i_
       policies::promote_double<false>,
       policies::discrete_quantile<>,
       policies::assert_undefined<> >::type forwarding_policy;
-   return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_i_prime_imp<value_type>(v, static_cast<value_type>(x), forwarding_policy()), "boost::math::cyl_bessel_i_prime<%1%>(%1%,%1%)");
+   return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_i_prime_imp<value_type>(static_cast<value_type>(v), static_cast<value_type>(x), forwarding_policy()), "boost::math::cyl_bessel_i_prime<%1%>(%1%,%1%)");
 }
 
 template <class T1, class T2>
@@ -301,7 +301,7 @@ inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_k_
       policies::promote_double<false>,
       policies::discrete_quantile<>,
       policies::assert_undefined<> >::type forwarding_policy;
-   return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_k_prime_imp<tag_type, value_type>(v, static_cast<value_type>(x), forwarding_policy()), "boost::math::cyl_bessel_k_prime<%1%,%1%>(%1%,%1%)");
+   return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_k_prime_imp<tag_type, value_type>(static_cast<value_type>(v), static_cast<value_type>(x), forwarding_policy()), "boost::math::cyl_bessel_k_prime<%1%,%1%>(%1%,%1%)");
 }
 
 template <class T1, class T2>
@@ -323,7 +323,7 @@ inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_neumann_p
       policies::promote_double<false>,
       policies::discrete_quantile<>,
       policies::assert_undefined<> >::type forwarding_policy;
-   return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_neumann_prime_imp<tag_type, value_type>(v, static_cast<value_type>(x), forwarding_policy()), "boost::math::cyl_neumann_prime<%1%,%1%>(%1%,%1%)");
+   return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_neumann_prime_imp<tag_type, value_type>(static_cast<value_type>(v), static_cast<value_type>(x), forwarding_policy()), "boost::math::cyl_neumann_prime<%1%,%1%>(%1%,%1%)");
 }
 
 template <class T1, class T2>
diff --git a/boost/math/special_functions/beta.hpp b/boost/math/special_functions/beta.hpp
index 98d8f7fa80..35b114ef15 100644
--- a/boost/math/special_functions/beta.hpp
+++ b/boost/math/special_functions/beta.hpp
@@ -47,13 +47,19 @@ T beta_imp(T a, T b, const Lanczos&, const Policy& pol)
 
    // Special cases:
    if((c == a) && (b < tools::epsilon<T>()))
-      return boost::math::tgamma(b, pol);
+      return 1 / b;
    else if((c == b) && (a < tools::epsilon<T>()))
-      return boost::math::tgamma(a, pol);
+      return 1 / a;
    if(b == 1)
       return 1/a;
    else if(a == 1)
       return 1/b;
+   else if(c < tools::epsilon<T>())
+   {
+      result = c / a;
+      result /= b;
+      return result;
+   }
 
    /*
    //
@@ -82,11 +88,11 @@ T beta_imp(T a, T b, const Lanczos&, const Policy& pol)
       std::swap(a, b);
 
    // Lanczos calculation:
-   T agh = a + Lanczos::g() - T(0.5);
-   T bgh = b + Lanczos::g() - T(0.5);
-   T cgh = c + Lanczos::g() - T(0.5);
-   result = Lanczos::lanczos_sum_expG_scaled(a) * Lanczos::lanczos_sum_expG_scaled(b) / Lanczos::lanczos_sum_expG_scaled(c);
-   T ambh = a - T(0.5) - b;
+   T agh = static_cast<T>(a + Lanczos::g() - 0.5f);
+   T bgh = static_cast<T>(b + Lanczos::g() - 0.5f);
+   T cgh = static_cast<T>(c + Lanczos::g() - 0.5f);
+   result = Lanczos::lanczos_sum_expG_scaled(a) * (Lanczos::lanczos_sum_expG_scaled(b) / Lanczos::lanczos_sum_expG_scaled(c));
+   T ambh = a - 0.5f - b;
    if((fabs(b * ambh) < (cgh * 100)) && (a > 100))
    {
       // Special case where the base of the power term is close to 1
@@ -199,7 +205,9 @@ T ibeta_power_terms(T a,
                         T y,
                         const Lanczos&,
                         bool normalised,
-                        const Policy& pol)
+                        const Policy& pol,
+                        T prefix = 1,
+                        const char* function = "boost::math::ibeta<%1%>(%1%, %1%, %1%)")
 {
    BOOST_MATH_STD_USING
 
@@ -211,14 +219,17 @@ T ibeta_power_terms(T a,
 
    T result;
 
-   T prefix = 1;
    T c = a + b;
 
    // combine power terms with Lanczos approximation:
-   T agh = a + Lanczos::g() - T(0.5);
-   T bgh = b + Lanczos::g() - T(0.5);
-   T cgh = c + Lanczos::g() - T(0.5);
+   T agh = static_cast<T>(a + Lanczos::g() - 0.5f);
+   T bgh = static_cast<T>(b + Lanczos::g() - 0.5f);
+   T cgh = static_cast<T>(c + Lanczos::g() - 0.5f);
    result = Lanczos::lanczos_sum_expG_scaled(c) / (Lanczos::lanczos_sum_expG_scaled(a) * Lanczos::lanczos_sum_expG_scaled(b));
+   result *= prefix;
+   // combine with the leftover terms from the Lanczos approximation:
+   result *= sqrt(bgh / boost::math::constants::e<T>());
+   result *= sqrt(agh / cgh);
 
    // l1 and l2 are the base of the exponents minus one:
    T l1 = (x * b - y * agh) / agh;
@@ -308,7 +319,15 @@ T ibeta_power_terms(T a,
          // First base near 1 only:
          T l = a * boost::math::log1p(l1, pol)
             + b * log((y * cgh) / bgh);
-         result *= exp(l);
+         if((l <= tools::log_min_value<T>()) || (l >= tools::log_max_value<T>()))
+         {
+            l += log(result);
+            if(l >= tools::log_max_value<T>())
+               return policies::raise_overflow_error<T>(function, 0, pol);
+            result = exp(l);
+         }
+         else
+            result *= exp(l);
          BOOST_MATH_INSTRUMENT_VARIABLE(result);
       }
       else
@@ -316,7 +335,15 @@ T ibeta_power_terms(T a,
          // Second base near 1 only:
          T l = b * boost::math::log1p(l2, pol)
             + a * log((x * cgh) / agh);
-         result *= exp(l);
+         if((l <= tools::log_min_value<T>()) || (l >= tools::log_max_value<T>()))
+         {
+            l += log(result);
+            if(l >= tools::log_max_value<T>())
+               return policies::raise_overflow_error<T>(function, 0, pol);
+            result = exp(l);
+         }
+         else
+            result *= exp(l);
          BOOST_MATH_INSTRUMENT_VARIABLE(result);
       }
    }
@@ -337,11 +364,41 @@ T ibeta_power_terms(T a,
          || (l2 <= tools::log_min_value<T>())
          )
       {
-         // Oops, overflow, sidestep:
+         // Oops, under/overflow, sidestep if we can:
          if(a < b)
-            result *= pow(pow(b2, b/a) * b1, a);
+         {
+            T p1 = pow(b2, b / a);
+            T l3 = a * (log(b1) + log(p1));
+            if((l3 < tools::log_max_value<T>())
+               && (l3 > tools::log_min_value<T>()))
+            {
+               result *= pow(p1 * b1, a);
+            }
+            else
+            {
+               l2 += l1 + log(result);
+               if(l2 >= tools::log_max_value<T>())
+                  return policies::raise_overflow_error<T>(function, 0, pol);
+               result = exp(l2);
+            }
+         }
          else
-            result *= pow(pow(b1, a/b) * b2, b);
+         {
+            T p1 = pow(b1, a / b);
+            T l3 = (log(p1) + log(b2)) * b;
+            if((l3 < tools::log_max_value<T>())
+               && (l3 > tools::log_min_value<T>()))
+            {
+               result *= pow(p1 * b2, b);
+            }
+            else
+            {
+               l2 += l1 + log(result);
+               if(l2 >= tools::log_max_value<T>())
+                  return policies::raise_overflow_error<T>(function, 0, pol);
+               result = exp(l2);
+            }
+         }
          BOOST_MATH_INSTRUMENT_VARIABLE(result);
       }
       else
@@ -351,10 +408,6 @@ T ibeta_power_terms(T a,
          BOOST_MATH_INSTRUMENT_VARIABLE(result);
       }
    }
-   // combine with the leftover terms from the Lanczos approximation:
-   result *= sqrt(bgh / boost::math::constants::e<T>());
-   result *= sqrt(agh / cgh);
-   result *= prefix;
 
    BOOST_MATH_INSTRUMENT_VARIABLE(result);
 
@@ -380,7 +433,9 @@ T ibeta_power_terms(T a,
                         T y,
                         const boost::math::lanczos::undefined_lanczos&,
                         bool normalised,
-                        const Policy& pol)
+                        const Policy& pol, 
+                        T prefix = 1,
+                        const char* = "boost::math::ibeta<%1%>(%1%, %1%, %1%)")
 {
    BOOST_MATH_STD_USING
 
@@ -411,7 +466,7 @@ T ibeta_power_terms(T a,
 
    T b1 = (x * lc) / la;
    T b2 = (y * lc) / lb;
-   T e1 = lc - la - lb;
+   T e1 = -5; // lc - la - lb;
    T lb1 = a * log(b1);
    T lb2 = b * log(b2);
 
@@ -423,21 +478,23 @@ T ibeta_power_terms(T a,
       || (e1 <= tools::log_min_value<T>())
       )
    {
-      result = exp(lb1 + lb2 - e1);
+      result = exp(lb1 + lb2 - e1 + log(prefix));
    }
    else
    {
       T p1, p2;
-      if((fabs(b1 - 1) * a < 10) && (a > 1))
-         p1 = exp(a * boost::math::log1p((x * b - y * la) / la, pol));
+      p1 = (x * b - y * la) / la;
+      if(fabs(p1) < 0.5f)
+         p1 = exp(a * boost::math::log1p(p1, pol));
       else
          p1 = pow(b1, a);
-      if((fabs(b2 - 1) * b < 10) && (b > 1))
-         p2 = exp(b * boost::math::log1p((y * a - x * lb) / lb, pol));
+      p2 = (y * a - x * lb) / lb;
+      if(fabs(p2) < 0.5f)
+         p2 = exp(b * boost::math::log1p(p2, pol));
       else
          p2 = pow(b2, b);
       T p3 = exp(e1);
-      result = p1 * p2 / p3;
+      result = prefix * p1 * (p2 / p3);
    }
    // and combine with the remaining gamma function components:
    result /= sa * sb / sc;
@@ -480,21 +537,43 @@ T ibeta_series(T a, T b, T x, T s0, const Lanczos&, bool normalised, T* p_deriva
       T c = a + b;
 
       // incomplete beta power term, combined with the Lanczos approximation:
-      T agh = a + Lanczos::g() - T(0.5);
-      T bgh = b + Lanczos::g() - T(0.5);
-      T cgh = c + Lanczos::g() - T(0.5);
+      T agh = static_cast<T>(a + Lanczos::g() - 0.5f);
+      T bgh = static_cast<T>(b + Lanczos::g() - 0.5f);
+      T cgh = static_cast<T>(c + Lanczos::g() - 0.5f);
       result = Lanczos::lanczos_sum_expG_scaled(c) / (Lanczos::lanczos_sum_expG_scaled(a) * Lanczos::lanczos_sum_expG_scaled(b));
-      if(a * b < bgh * 10)
-         result *= exp((b - 0.5f) * boost::math::log1p(a / bgh, pol));
-      else
-         result *= pow(cgh / bgh, b - 0.5f);
-      result *= pow(x * cgh / agh, a);
-      result *= sqrt(agh / boost::math::constants::e<T>());
 
-      if(p_derivative)
+      T l1 = log(cgh / bgh) * (b - 0.5f);
+      T l2 = log(x * cgh / agh) * a;
+      //
+      // Check for over/underflow in the power terms:
+      //
+      if((l1 > tools::log_min_value<T>())
+         && (l1 < tools::log_max_value<T>())
+         && (l2 > tools::log_min_value<T>())
+         && (l2 < tools::log_max_value<T>()))
       {
-         *p_derivative = result * pow(y, b);
-         BOOST_ASSERT(*p_derivative >= 0);
+         if(a * b < bgh * 10)
+            result *= exp((b - 0.5f) * boost::math::log1p(a / bgh, pol));
+         else
+            result *= pow(cgh / bgh, b - 0.5f);
+         result *= pow(x * cgh / agh, a);
+         result *= sqrt(agh / boost::math::constants::e<T>());
+
+         if(p_derivative)
+         {
+            *p_derivative = result * pow(y, b);
+            BOOST_ASSERT(*p_derivative >= 0);
+         }
+      }
+      else
+      {
+         //
+         // Oh dear, we need logs, and this *will* cancel:
+         //
+         result = log(result) + l1 + l2 + (log(agh) - 1) / 2;
+         if(p_derivative)
+            *p_derivative = exp(result + b * log(y));
+         result = exp(result);
       }
    }
    else
@@ -606,7 +685,7 @@ struct ibeta_fraction2_t
       T denom = (a + 2 * m - 1);
       aN /= denom * denom;
 
-      T bN = m;
+      T bN = static_cast<T>(m);
       bN += (m * (b - m) * x) / (a + 2*m - 1);
       bN += ((a + m) * (a * y - b * x + 1 + m *(2 - x))) / (a + 2*m + 1);
 
@@ -943,12 +1022,12 @@ T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised, T* p_de
          if(b == 0)
             return policies::raise_domain_error<T>(function, "The arguments a and b to the incomplete beta function cannot both be zero, with x=%1%.", x, pol);
          if(b > 0)
-            return inv ? 0 : 1;
+            return static_cast<T>(inv ? 0 : 1);
       }
       else if(b == 0)
       {
          if(a > 0)
-            return inv ? 1 : 0;
+            return static_cast<T>(inv ? 1 : 0);
       }
    }
    else
@@ -1200,7 +1279,7 @@ T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised, T* p_de
       
       if(b < 40)
       {
-         if((floor(a) == a) && (floor(b) == b) && (a < (std::numeric_limits<int>::max)() - 100))
+         if((floor(a) == a) && (floor(b) == b) && (a < (std::numeric_limits<int>::max)() - 100) && (y != 1))
          {
             // relate to the binomial distribution and use a finite sum:
             T k = a - 1;
@@ -1342,20 +1421,8 @@ T ibeta_derivative_imp(T a, T b, T x, const Policy& pol)
    // Now the regular cases:
    //
    typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
-   T f1 = ibeta_power_terms<T>(a, b, x, 1 - x, lanczos_type(), true, pol);
    T y = (1 - x) * x;
-
-   if(f1 == 0)
-      return 0;
-   
-   if((tools::max_value<T>() * y < f1))
-   {
-      // overflow:
-      return policies::raise_overflow_error<T>(function, 0, pol);
-   }
-
-   f1 /= y;
-
+   T f1 = ibeta_power_terms<T>(a, b, x, 1 - x, lanczos_type(), true, pol, 1 / y, function);
    return f1;
 }
 //
diff --git a/boost/math/special_functions/detail/bernoulli_details.hpp b/boost/math/special_functions/detail/bernoulli_details.hpp
index 525c1fccf4..c12a172b60 100644
--- a/boost/math/special_functions/detail/bernoulli_details.hpp
+++ b/boost/math/special_functions/detail/bernoulli_details.hpp
@@ -202,9 +202,13 @@ struct bernoulli_initializer
          // initialize our dymanic table:
          //
          boost::math::bernoulli_b2n<T>(2, Policy());
+#ifndef BOOST_NO_EXCEPTIONS
          try{
+#endif
             boost::math::bernoulli_b2n<T>(max_bernoulli_b2n<T>::value + 1, Policy());
+#ifndef BOOST_NO_EXCEPTIONS
          } catch(const std::overflow_error&){}
+#endif
          boost::math::tangent_t2n<T>(2, Policy());
       }
       void force_instantiate()const{}
@@ -254,7 +258,9 @@ struct fixed_vector : private std::allocator<T>
    void resize(unsigned n, const T& val)
    {
       if(n > m_capacity)
-         throw std::runtime_error("Exhausted storage for Bernoulli numbers.");
+      {
+         BOOST_THROW_EXCEPTION(std::runtime_error("Exhausted storage for Bernoulli numbers."));
+      }
       for(unsigned i = m_used; i < n; ++i)
          new (m_data + i) T(val);
       m_used = n;
@@ -432,11 +438,11 @@ public:
       //
       if(start + n >= bn.size())
       {
-         std::size_t new_size = (std::min)((std::max)((std::max)(start + n, std::size_t(bn.size() + 20)), std::size_t(50)), std::size_t(bn.capacity()));
+         std::size_t new_size = (std::min)((std::max)((std::max)(std::size_t(start + n), std::size_t(bn.size() + 20)), std::size_t(50)), std::size_t(bn.capacity()));
          tangent_numbers_series(new_size);
       }
 
-      for(std::size_t i = (std::max)(max_bernoulli_b2n<T>::value + 1, start); i < start + n; ++i)
+      for(std::size_t i = (std::max)(std::size_t(max_bernoulli_b2n<T>::value + 1), start); i < start + n; ++i)
       {
          *out = (i >= m_overflow_limit) ? policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol) : bn[i];
          ++out;
@@ -448,11 +454,11 @@ public:
       boost::detail::lightweight_mutex::scoped_lock l(m_mutex);
       if(start + n >= bn.size())
       {
-         std::size_t new_size = (std::min)((std::max)((std::max)(start + n, std::size_t(bn.size() + 20)), std::size_t(50)), std::size_t(bn.capacity()));
+         std::size_t new_size = (std::min)((std::max)((std::max)(std::size_t(start + n), std::size_t(bn.size() + 20)), std::size_t(50)), std::size_t(bn.capacity()));
          tangent_numbers_series(new_size);
       }
 
-      for(std::size_t i = (std::max)(max_bernoulli_b2n<T>::value + 1, start); i < start + n; ++i)
+      for(std::size_t i = (std::max)(std::size_t(max_bernoulli_b2n<T>::value + 1), start); i < start + n; ++i)
       {
          *out = (i >= m_overflow_limit) ? policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol) : bn[i];
          ++out;
@@ -473,7 +479,7 @@ public:
          {
             if(start + n >= bn.size())
             {
-               std::size_t new_size = (std::min)((std::max)((std::max)(start + n, std::size_t(bn.size() + 20)), std::size_t(50)), std::size_t(bn.capacity()));
+               std::size_t new_size = (std::min)((std::max)((std::max)(std::size_t(start + n), std::size_t(bn.size() + 20)), std::size_t(50)), std::size_t(bn.capacity()));
                tangent_numbers_series(new_size);
             }
             m_counter.store(static_cast<atomic_integer_type>(bn.size()), BOOST_MATH_ATOMIC_NS::memory_order_release);
diff --git a/boost/math/special_functions/detail/bessel_derivatives_linear.hpp b/boost/math/special_functions/detail/bessel_derivatives_linear.hpp
index 2ee86a03ee..72fddc61c6 100644
--- a/boost/math/special_functions/detail/bessel_derivatives_linear.hpp
+++ b/boost/math/special_functions/detail/bessel_derivatives_linear.hpp
@@ -37,13 +37,25 @@ inline T sph_bessel_j_derivative_linear(unsigned v, T x, Policy pol)
 template <class T, class Policy>
 inline T bessel_i_derivative_linear(T v, T x, Policy pol)
 {
-   return (boost::math::detail::cyl_bessel_i_imp<T>(v-1, x, pol) + boost::math::detail::cyl_bessel_i_imp<T>(v+1, x, pol)) / 2;
+   T result = boost::math::detail::cyl_bessel_i_imp<T>(v - 1, x, pol);
+   if(result >= tools::max_value<T>())
+      return result;  // result is infinite
+   T result2 = boost::math::detail::cyl_bessel_i_imp<T>(v + 1, x, pol);
+   if(result2 >= tools::max_value<T>() - result)
+      return result2;  // result is infinite
+   return (result + result2) / 2;
 }
 
 template <class T, class Tag, class Policy>
 inline T bessel_k_derivative_linear(T v, T x, Tag tag, Policy pol)
 {
-   return (boost::math::detail::cyl_bessel_k_imp<T>(v-1, x, tag, pol) + boost::math::detail::cyl_bessel_k_imp<T>(v+1, x, tag, pol)) / -2;
+   T result = boost::math::detail::cyl_bessel_k_imp<T>(v - 1, x, tag, pol);
+   if(result >= tools::max_value<T>())
+      return -result;  // result is infinite
+   T result2 = boost::math::detail::cyl_bessel_k_imp<T>(v + 1, x, tag, pol);
+   if(result2 >= tools::max_value<T>() - result)
+      return -result2;  // result is infinite
+   return (result + result2) / -2;
 }
 
 template <class T, class Policy>
diff --git a/boost/math/special_functions/detail/bessel_ik.hpp b/boost/math/special_functions/detail/bessel_ik.hpp
index 10118d9715..aac1781e10 100644
--- a/boost/math/special_functions/detail/bessel_ik.hpp
+++ b/boost/math/special_functions/detail/bessel_ik.hpp
@@ -374,6 +374,7 @@ int bessel_ik(T v, T x, T* I, T* K, int kind, const Policy& pol)
     prev = Ku;
     current = Ku1;
     T scale = 1;
+    T scale_sign = 1;
     for (k = 1; k <= n; k++)                   // forward recurrence for K
     {
         T fact = 2 * (u + k) / x;
@@ -381,6 +382,7 @@ int bessel_ik(T v, T x, T* I, T* K, int kind, const Policy& pol)
         {
            prev /= current;
            scale /= current;
+           scale_sign *= boost::math::sign(current);
            current = 1;
         }
         next = fact * current + prev;
@@ -426,7 +428,7 @@ int bessel_ik(T v, T x, T* I, T* K, int kind, const Policy& pol)
         if(fact == 0)
            *I = Iv;
         else if(tools::max_value<T>() * scale < fact)
-           *I = (org_kind & need_i) ? T(sign(fact) * sign(scale) * policies::raise_overflow_error<T>(function, 0, pol)) : T(0);
+           *I = (org_kind & need_i) ? T(sign(fact) * scale_sign * policies::raise_overflow_error<T>(function, 0, pol)) : T(0);
         else
          *I = Iv + fact / scale;   // reflection formula
     }
@@ -435,7 +437,7 @@ int bessel_ik(T v, T x, T* I, T* K, int kind, const Policy& pol)
         *I = Iv;
     }
     if(tools::max_value<T>() * scale < Kv)
-      *K = (org_kind & need_k) ? T(sign(Kv) * sign(scale) * policies::raise_overflow_error<T>(function, 0, pol)) : T(0);
+       *K = (org_kind & need_k) ? T(sign(Kv) * scale_sign * policies::raise_overflow_error<T>(function, 0, pol)) : T(0);
     else
       *K = Kv / scale;
     BOOST_MATH_INSTRUMENT_VARIABLE(*I);
diff --git a/boost/math/special_functions/detail/bessel_jn.hpp b/boost/math/special_functions/detail/bessel_jn.hpp
index 3f15f9cd87..2413630637 100644
--- a/boost/math/special_functions/detail/bessel_jn.hpp
+++ b/boost/math/special_functions/detail/bessel_jn.hpp
@@ -35,7 +35,7 @@ T bessel_jn(int n, T x, const Policy& pol)
     //
     if (n < 0)
     {
-        factor = (n & 0x1) ? -1 : 1;  // J_{-n}(z) = (-1)^n J_n(z)
+        factor = static_cast<T>((n & 0x1) ? -1 : 1);  // J_{-n}(z) = (-1)^n J_n(z)
         n = -n;
     }
     else
@@ -50,6 +50,8 @@ T bessel_jn(int n, T x, const Policy& pol)
     //
     // Special cases:
     //
+    if(asymptotic_bessel_large_x_limit(T(n), x))
+       return factor * asymptotic_bessel_j_large_x_2<T>(T(n), x);
     if (n == 0)
     {
         return factor * bessel_j0(x);
@@ -64,9 +66,6 @@ T bessel_jn(int n, T x, const Policy& pol)
         return static_cast<T>(0);
     }
 
-    if(asymptotic_bessel_large_x_limit(T(n), x))
-      return factor * asymptotic_bessel_j_large_x_2<T>(n, x);
-
     BOOST_ASSERT(n > 1);
     T scale = 1;
     if (n < abs(x))                         // forward recurrence
diff --git a/boost/math/special_functions/detail/bessel_jy_asym.hpp b/boost/math/special_functions/detail/bessel_jy_asym.hpp
index 81f6238e58..4d7ac485ad 100644
--- a/boost/math/special_functions/detail/bessel_jy_asym.hpp
+++ b/boost/math/special_functions/detail/bessel_jy_asym.hpp
@@ -120,6 +120,24 @@ inline T asymptotic_bessel_j_large_x_2(T v, T x)
 }
 
 template <class T>
+inline bool asymptotic_bessel_large_x_limit(int v, const T& x)
+{
+   BOOST_MATH_STD_USING
+      //
+      // Determines if x is large enough compared to v to take the asymptotic
+      // forms above.  From A&S 9.2.28 we require: 
+      //    v < x * eps^1/8
+      // and from A&S 9.2.29 we require:
+      //    v^12/10 < 1.5 * x * eps^1/10
+      // using the former seems to work OK in practice with broadly similar 
+      // error rates either side of the divide for v < 10000.
+      // At double precision eps^1/8 ~= 0.01.
+      //
+      BOOST_ASSERT(v >= 0);
+      return (v ? v : 1) < x * 0.004f;
+}
+
+template <class T>
 inline bool asymptotic_bessel_large_x_limit(const T& v, const T& x)
 {
    BOOST_MATH_STD_USING
diff --git a/boost/math/special_functions/detail/bessel_jy_derivatives_series.hpp b/boost/math/special_functions/detail/bessel_jy_derivatives_series.hpp
index 0dc68fc73c..5425caf689 100644
--- a/boost/math/special_functions/detail/bessel_jy_derivatives_series.hpp
+++ b/boost/math/special_functions/detail/bessel_jy_derivatives_series.hpp
@@ -157,7 +157,8 @@ inline T bessel_y_derivative_small_z_series(T v, T x, const Policy& pol)
          gam = 1;
          if (boost::math::tools::max_value<T>() * p < gam)
          {
-            return -boost::math::policies::raise_overflow_error<T>(function, 0, pol);
+            // This term will overflow to -INF, when combined with the series below it becomes +INF:
+            return boost::math::policies::raise_overflow_error<T>(function, 0, pol);
          }
       }
       prefix = -gam / (boost::math::constants::pi<T>() * p);
@@ -173,7 +174,7 @@ inline T bessel_y_derivative_small_z_series(T v, T x, const Policy& pol)
          scale /= (boost::math::tools::max_value<T>() / 4);
          if (boost::math::tools::log_max_value<T>() < prefix)
          {
-            return -boost::math::policies::raise_overflow_error<T>(function, 0, pol);
+            return boost::math::policies::raise_overflow_error<T>(function, 0, pol);
          }
       }
       prefix = -exp(prefix);
diff --git a/boost/math/special_functions/detail/bessel_yn.hpp b/boost/math/special_functions/detail/bessel_yn.hpp
index 0509062bbd..62d7377e4f 100644
--- a/boost/math/special_functions/detail/bessel_yn.hpp
+++ b/boost/math/special_functions/detail/bessel_yn.hpp
@@ -45,14 +45,13 @@ T bessel_yn(int n, T x, const Policy& pol)
     //
     if (n < 0)
     {
-        factor = (n & 0x1) ? -1 : 1;  // Y_{-n}(z) = (-1)^n Y_n(z)
+        factor = static_cast<T>((n & 0x1) ? -1 : 1);  // Y_{-n}(z) = (-1)^n Y_n(z)
         n = -n;
     }
     else
     {
         factor = 1;
     }
-
     if(x < policies::get_epsilon<T, Policy>())
     {
        T scale = 1;
@@ -61,6 +60,10 @@ T bessel_yn(int n, T x, const Policy& pol)
           return boost::math::sign(scale) * boost::math::sign(value) * policies::raise_overflow_error<T>(function, 0, pol);
        value /= scale;
     }
+    else if(asymptotic_bessel_large_x_limit(n, x))
+    {
+       value = factor * asymptotic_bessel_y_large_x_2(static_cast<T>(abs(n)), x);
+    }
     else if (n == 0)
     {
         value = bessel_y0(x, pol);
@@ -76,23 +79,28 @@ T bessel_yn(int n, T x, const Policy& pol)
        int k = 1;
        BOOST_ASSERT(k < n);
        policies::check_series_iterations<T>("boost::math::bessel_y_n<%1%>(%1%,%1%)", n, pol);
-       do
+       T mult = 2 * k / x;
+       value = mult * current - prev;
+       prev = current;
+       current = value;
+       ++k;
+       if((mult > 1) && (fabs(current) > 1))
+       {
+          prev /= current;
+          factor /= current;
+          value /= current;
+          current = 1;
+       }
+       while(k < n)
        {
-           T fact = 2 * k / x;
-           if((fact > 1) && ((tools::max_value<T>() - fabs(prev)) / fact < fabs(current)))
-           {
-              prev /= current;
-              factor /= current;
-              current = 1;
-           }
-           value = fact * current - prev;
+           mult = 2 * k / x;
+           value = mult * current - prev;
            prev = current;
            current = value;
            ++k;
        }
-       while(k < n);
        if(fabs(tools::max_value<T>() * factor) < fabs(value))
-          return sign(value) * sign(value) * policies::raise_overflow_error<T>(function, 0, pol);
+          return sign(value) * sign(factor) * policies::raise_overflow_error<T>(function, 0, pol);
        value /= factor;
     }
     return value;
diff --git a/boost/math/special_functions/detail/lanczos_sse2.hpp b/boost/math/special_functions/detail/lanczos_sse2.hpp
index edef3a0412..b2dd7ea292 100644
--- a/boost/math/special_functions/detail/lanczos_sse2.hpp
+++ b/boost/math/special_functions/detail/lanczos_sse2.hpp
@@ -12,7 +12,7 @@
 
 #include <emmintrin.h>
 
-#if defined(__GNUC__) || defined(__PGI)
+#if defined(__GNUC__) || defined(__PGI) || defined(__SUNPRO_CC)
 #define ALIGN16 __attribute__((__aligned__(16)))
 #else
 #define ALIGN16 __declspec(align(16))
diff --git a/boost/math/special_functions/detail/polygamma.hpp b/boost/math/special_functions/detail/polygamma.hpp
index 4ef503bf7c..0267bd3ac8 100644
--- a/boost/math/special_functions/detail/polygamma.hpp
+++ b/boost/math/special_functions/detail/polygamma.hpp
@@ -44,7 +44,7 @@
         // x is crazy large, just concentrate on the first part of the expression and use logs:
         if(n == 1) return 1 / x;
         T nlx = n * log(x);
-        if((nlx < tools::log_max_value<T>()) && (n < max_factorial<T>::value))
+        if((nlx < tools::log_max_value<T>()) && (n < (int)max_factorial<T>::value))
            return ((n & 1) ? 1 : -1) * boost::math::factorial<T>(n - 1) * pow(x, -n);
         else
          return ((n & 1) ? 1 : -1) * exp(boost::math::lgamma(T(n), pol) - n * log(x));
@@ -71,13 +71,13 @@
      // or the power term underflows, this just gets set to 0 and then we
      // know that we have to use logs for the initial terms:
      //
-     part_term = ((n > boost::math::max_factorial<T>::value) && (T(n) * n > tools::log_max_value<T>())) 
+     part_term = ((n > (int)boost::math::max_factorial<T>::value) && (T(n) * n > tools::log_max_value<T>())) 
         ? T(0) : static_cast<T>(boost::math::factorial<T>(n - 1, pol) * pow(x, -n - 1));
      if(part_term == 0)
      {
         // Either n is very large, or the power term underflows,
         // set the initial values of part_term, term and sum via logs:
-        part_term = boost::math::lgamma(n, pol) - (n + 1) * log(x);
+        part_term = static_cast<T>(boost::math::lgamma(n, pol) - (n + 1) * log(x));
         sum = exp(part_term + log(n + 2 * x) - boost::math::constants::ln_two<T>());
         part_term += log(T(n) * (n + 1)) - boost::math::constants::ln_two<T>() - log(x);
         part_term = exp(part_term);
@@ -505,7 +505,7 @@
     // would mean setting the limit to ~ 1 / n,
     // but we can tolerate a small amount of divergence:
     //
-    T small_x_limit = std::min(T(T(5) / n), T(0.25f));
+    T small_x_limit = (std::min)(T(T(5) / n), T(0.25f));
     if(x < small_x_limit)
     {
       return polygamma_nearzero(n, x, pol, function);
diff --git a/boost/math/special_functions/detail/t_distribution_inv.hpp b/boost/math/special_functions/detail/t_distribution_inv.hpp
index 72f6f0c646..e8bd0db28f 100644
--- a/boost/math/special_functions/detail/t_distribution_inv.hpp
+++ b/boost/math/special_functions/detail/t_distribution_inv.hpp
@@ -56,9 +56,9 @@ T inverse_students_t_hill(T ndf, T u, const Policy& pol)
    }
    else
    {
-      y = ((1 / (((ndf + 6) / (ndf * y) - 0.089f * d - 0.822f)
+      y = static_cast<T>(((1 / (((ndf + 6) / (ndf * y) - 0.089f * d - 0.822f)
               * (ndf + 2) * 3) + 0.5 / (ndf + 4)) * y - 1)
-              * (ndf + 1) / (ndf + 2) + 1 / y;
+              * (ndf + 1) / (ndf + 2) + 1 / y);
    }
    q = sqrt(ndf * y);
 
@@ -144,41 +144,41 @@ T inverse_students_t_body_series(T df, T u, const Policy& pol)
    // only on the degrees of freedom (Eq 57 of Shaw):
    //
    T in = 1 / df;
-   c[2] = 0.16666666666666666667 + 0.16666666666666666667 * in;
-   c[3] = (0.0083333333333333333333 * in 
+   c[2] = static_cast<T>(0.16666666666666666667 + 0.16666666666666666667 * in);
+   c[3] = static_cast<T>((0.0083333333333333333333 * in
       + 0.066666666666666666667) * in 
-      + 0.058333333333333333333;
-   c[4] = ((0.00019841269841269841270 * in 
+      + 0.058333333333333333333);
+   c[4] = static_cast<T>(((0.00019841269841269841270 * in
       + 0.0017857142857142857143) * in 
       + 0.026785714285714285714) * in 
-      + 0.025198412698412698413;
-   c[5] = (((2.7557319223985890653e-6 * in 
+      + 0.025198412698412698413);
+   c[5] = static_cast<T>((((2.7557319223985890653e-6 * in
       + 0.00037477954144620811287) * in 
       - 0.0011078042328042328042) * in 
       + 0.010559964726631393298) * in 
-      + 0.012039792768959435626;
-   c[6] = ((((2.5052108385441718775e-8 * in 
+      + 0.012039792768959435626);
+   c[6] = static_cast<T>(((((2.5052108385441718775e-8 * in
       - 0.000062705427288760622094) * in 
       + 0.00059458674042007375341) * in 
       - 0.0016095979637646304313) * in 
       + 0.0061039211560044893378) * in 
-      + 0.0038370059724226390893;
-   c[7] = (((((1.6059043836821614599e-10 * in 
+      + 0.0038370059724226390893);
+   c[7] = static_cast<T>((((((1.6059043836821614599e-10 * in
       + 0.000015401265401265401265) * in 
       - 0.00016376804137220803887) * in
       + 0.00069084207973096861986) * in 
       - 0.0012579159844784844785) * in 
       + 0.0010898206731540064873) * in 
-      + 0.0032177478835464946576;
-   c[8] = ((((((7.6471637318198164759e-13 * in
+      + 0.0032177478835464946576);
+   c[8] = static_cast<T>(((((((7.6471637318198164759e-13 * in
       - 3.9851014346715404916e-6) * in
       + 0.000049255746366361445727) * in
       - 0.00024947258047043099953) * in 
       + 0.00064513046951456342991) * in
       - 0.00076245135440323932387) * in
       + 0.000033530976880017885309) * in 
-      + 0.0017438262298340009980;
-   c[9] = (((((((2.8114572543455207632e-15 * in
+      + 0.0017438262298340009980);
+   c[9] = static_cast<T>((((((((2.8114572543455207632e-15 * in
       + 1.0914179173496789432e-6) * in
       - 0.000015303004486655377567) * in
       + 0.000090867107935219902229) * in
@@ -186,8 +186,8 @@ T inverse_students_t_body_series(T df, T u, const Policy& pol)
       + 0.00051406605788341121363) * in
       - 0.00036307660358786885787) * in
       - 0.00031101086326318780412) * in 
-      + 0.00096472747321388644237;
-   c[10] = ((((((((8.2206352466243297170e-18 * in
+      + 0.00096472747321388644237);
+   c[10] = static_cast<T>(((((((((8.2206352466243297170e-18 * in
       - 3.1239569599829868045e-7) * in
       + 4.8903045291975346210e-6) * in
       - 0.000033202652391372058698) * in
@@ -196,7 +196,7 @@ T inverse_students_t_body_series(T df, T u, const Policy& pol)
       + 0.00035764655430568632777) * in
       - 0.00010230378073700412687) * in
       - 0.00036942667800009661203) * in
-      + 0.00054229262813129686486;
+      + 0.00054229262813129686486);
    //
    // The result is then a polynomial in v (see Eq 56 of Shaw):
    //
@@ -285,7 +285,7 @@ T inverse_students_t(T df, T u, T v, const Policy& pol, bool* pexact = 0)
             //
             T a = 4 * (u - u * u);//1 - 4 * (u - 0.5f) * (u - 0.5f);
             T b = boost::math::cbrt(a);
-            static const T c = 0.85498797333834849467655443627193;
+            static const T c = static_cast<T>(0.85498797333834849467655443627193L);
             T p = 6 * (1 + c * (1 / b - 1));
             T p0;
             do{
diff --git a/boost/math/special_functions/ellint_2.hpp b/boost/math/special_functions/ellint_2.hpp
index f4f65cc11d..0d160ecdf1 100644
--- a/boost/math/special_functions/ellint_2.hpp
+++ b/boost/math/special_functions/ellint_2.hpp
@@ -103,17 +103,10 @@ T ellint_e_imp(T phi, T k, const Policy& pol)
        {
           return policies::raise_domain_error<T>("boost::math::ellint_2<%1%>(%1%, %1%)", "The parameter k is out of range, got k = %1%", k, pol);
        }
-       else if(rphi == 0)
+       else if(rphi < tools::root_epsilon<T>())
        {
-          result = 0;
-       }
-       else if(sinp * sinp < tools::min_value<T>())
-       {
-          T x = cosp * cosp;
-          T t = k * k * sinp * sinp;
-          T y = 1 - t;
-          T z = 1;
-          result = s * sinp * (ellint_rf_imp(x, y, z, pol) - t * ellint_rd_imp(x, y, z, pol) / 3);
+          // See http://functions.wolfram.com/EllipticIntegrals/EllipticE2/06/01/03/0001/
+          result = s * rphi;
        }
        else
        {
diff --git a/boost/math/special_functions/ellint_d.hpp b/boost/math/special_functions/ellint_d.hpp
index 5bd065d6e3..bc5a4b2a56 100644
--- a/boost/math/special_functions/ellint_d.hpp
+++ b/boost/math/special_functions/ellint_d.hpp
@@ -60,7 +60,7 @@ T ellint_d_imp(T phi, T k, const Policy& pol)
     if(phi >= tools::max_value<T>())
     {
        // Need to handle infinity as a special case:
-       result = policies::raise_overflow_error<T>("boost::math::ellint_e<%1%>(%1%,%1%)", 0, pol);
+       result = policies::raise_overflow_error<T>("boost::math::ellint_d<%1%>(%1%,%1%)", 0, pol);
     }
     else if(phi > 1 / tools::epsilon<T>())
     {
@@ -113,15 +113,11 @@ T ellint_d_imp(T k, const Policy& pol)
     BOOST_MATH_STD_USING
     using namespace boost::math::tools;
 
-    if (abs(k) > 1)
+    if (abs(k) >= 1)
     {
-       return policies::raise_domain_error<T>("boost::math::ellint_e<%1%>(%1%)",
+       return policies::raise_domain_error<T>("boost::math::ellint_d<%1%>(%1%)",
             "Got k = %1%, function requires |k| <= 1", k, pol);
     }
-    if (abs(k) == 1)
-    {
-        return static_cast<T>(1);
-    }
     if(fabs(k) <= tools::root_epsilon<T>())
        return constants::pi<T>() / 4;
 
diff --git a/boost/math/special_functions/fpclassify.hpp b/boost/math/special_functions/fpclassify.hpp
index 8e75fae0f2..0a4e1ac7aa 100644
--- a/boost/math/special_functions/fpclassify.hpp
+++ b/boost/math/special_functions/fpclassify.hpp
@@ -80,6 +80,9 @@ is used.
 #if defined(_MSC_VER) || defined(__BORLANDC__)
 #include <float.h>
 #endif
+#ifdef BOOST_MATH_USE_FLOAT128
+#include "quadmath.h"
+#endif
 
 #ifdef BOOST_NO_STDC_NAMESPACE
   namespace std{ using ::abs; using ::fabs; }
@@ -121,7 +124,10 @@ inline bool is_nan_helper(T, const boost::false_type&)
 {
    return false;
 }
-
+#ifdef BOOST_MATH_USE_FLOAT128
+inline bool is_nan_helper(__float128 f, const boost::true_type&) { return ::isnanq(f); }
+inline bool is_nan_helper(__float128 f, const boost::false_type&) { return ::isnanq(f); }
+#endif
 }
 
 namespace math{
@@ -513,6 +519,13 @@ inline bool (isinf)(long double x)
    return detail::isinf_impl(static_cast<value_type>(x), method());
 }
 #endif
+#ifdef BOOST_MATH_USE_FLOAT128
+template<>
+inline bool (isinf)(__float128 x)
+{
+   return ::isinfq(x);
+}
+#endif
 
 //------------------------------------------------------------------------------
 
@@ -598,6 +611,13 @@ inline bool (isnan)(long double x)
    return detail::isnan_impl(x, method());
 }
 #endif
+#ifdef BOOST_MATH_USE_FLOAT128
+template<>
+inline bool (isnan)(__float128 x)
+{
+   return ::isnanq(x);
+}
+#endif
 
 } // namespace math
 } // namespace boost
diff --git a/boost/math/special_functions/gamma.hpp b/boost/math/special_functions/gamma.hpp
index 1db1e4c2d3..3a3191a807 100644
--- a/boost/math/special_functions/gamma.hpp
+++ b/boost/math/special_functions/gamma.hpp
@@ -162,7 +162,7 @@ T gamma_imp(T z, const Policy& pol, const Lanczos& l)
    {
       if (z < 1 / tools::max_value<T>())
          result = policies::raise_overflow_error<T>(function, 0, pol);
-	   result *= 1 / z - constants::euler<T>();
+      result *= 1 / z - constants::euler<T>();
    }
    else
    {
@@ -238,13 +238,13 @@ T lgamma_imp(T z, const Policy& pol, const Lanczos& l, int* sign = 0)
    }
    else if (z < tools::root_epsilon<T>())
    {
-	   if (0 == z)
-		   return policies::raise_pole_error<T>(function, "Evaluation of lgamma at %1%.", z, pol);
+     if (0 == z)
+	      return policies::raise_pole_error<T>(function, "Evaluation of lgamma at %1%.", z, pol);
       if (fabs(z) < 1 / tools::max_value<T>())
          result = -log(fabs(z));
       else
-	      result = log(fabs(1 / z - constants::euler<T>()));
-	   if (z < 0)
+         result = log(fabs(1 / z - constants::euler<T>()));
+      if (z < 0)
 		sresult = -1;
    }
    else if(z < 15)
@@ -526,9 +526,9 @@ T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sig
 
    if (zz < min_arg_for_recursion)
    {
-	   // Here we simply take the logarithm of tgamma(). This is somewhat
-	   // inefficient, but simple. The rationale is that the argument here
-	   // is relatively small and overflow is not expected to be likely.
+      // Here we simply take the logarithm of tgamma(). This is somewhat
+      // inefficient, but simple. The rationale is that the argument here
+      // is relatively small and overflow is not expected to be likely.
       if (z > -tools::root_epsilon<T>())
       {
          // Reflection formula may fail if z is very close to zero, let the series
@@ -540,7 +540,7 @@ T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sig
          }
          return log_gamma_value;
       }
-	   else
+      else
       {
          // No issue with spurious overflow in reflection formula, 
          // just fall through to regular code:
@@ -1394,7 +1394,7 @@ T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const Lanczos&
          return 1 / (z * boost::math::tgamma(z + delta, pol));
       }
    }
-   T zgh = z + Lanczos::g() - constants::half<T>();
+   T zgh = static_cast<T>(z + Lanczos::g() - constants::half<T>());
    T result;
    if(fabs(delta) < 10)
    {
diff --git a/boost/math/special_functions/math_fwd.hpp b/boost/math/special_functions/math_fwd.hpp
index 96f60726a6..ca8e58144e 100644
--- a/boost/math/special_functions/math_fwd.hpp
+++ b/boost/math/special_functions/math_fwd.hpp
@@ -181,19 +181,19 @@ namespace boost
    template <class T>
    typename tools::promote_args<T>::type
          legendre_p(int l, T x);
-
+#if !BOOST_WORKAROUND(BOOST_MSVC, <= 1310)
    template <class T, class Policy>
    typename boost::enable_if_c<policies::is_policy<Policy>::value, typename tools::promote_args<T>::type>::type
          legendre_p(int l, T x, const Policy& pol);
-
+#endif
    template <class T>
    typename tools::promote_args<T>::type
          legendre_q(unsigned l, T x);
-
+#if !BOOST_WORKAROUND(BOOST_MSVC, <= 1310)
    template <class T, class Policy>
    typename boost::enable_if_c<policies::is_policy<Policy>::value, typename tools::promote_args<T>::type>::type
          legendre_q(unsigned l, T x, const Policy& pol);
-
+#endif
    template <class T1, class T2, class T3>
    typename tools::promote_args<T1, T2, T3>::type
          legendre_next(unsigned l, unsigned m, T1 x, T2 Pl, T3 Plm1);
@@ -607,8 +607,10 @@ namespace boost
       template <class T1, class T2, class Policy>
       struct bessel_traits
       {
-         typedef typename tools::promote_args<
-            T1, T2
+         typedef typename mpl::if_<
+            is_integral<T1>,
+            typename tools::promote_args<T2>::type,
+            typename tools::promote_args<T1, T2>::type
          >::type result_type;
 
          typedef typename policies::precision<result_type, Policy>::type precision_type;
@@ -994,6 +996,16 @@ namespace boost
    template <class T>
    typename tools::promote_args<T>::type float_advance(const T& val, int distance);
 
+   template <class T, class Policy>
+   typename tools::promote_args<T>::type ulp(const T& val, const Policy& pol);
+   template <class T>
+   typename tools::promote_args<T>::type ulp(const T& val);
+
+   template <class T, class U>
+   typename tools::promote_args<T, U>::type relative_difference(const T&, const U&);
+   template <class T, class U>
+   typename tools::promote_args<T, U>::type epsilon_difference(const T&, const U&);
+
    template<class T>
    T unchecked_bernoulli_b2n(const std::size_t n);
    template <class T, class Policy>
@@ -1447,6 +1459,7 @@ template <class OutputIterator, class T>\
    template <class T> T float_next(const T& a){ return boost::math::float_next(a, Policy()); }\
    template <class T> T float_prior(const T& a){ return boost::math::float_prior(a, Policy()); }\
    template <class T> T float_distance(const T& a, const T& b){ return boost::math::float_distance(a, b, Policy()); }\
+   template <class T> T ulp(const T& a){ return boost::math::ulp(a, Policy()); }\
    \
    template <class RT1, class RT2>\
    inline typename boost::math::tools::promote_args<RT1, RT2>::type owens_t(RT1 a, RT2 z){ return boost::math::owens_t(a, z, Policy()); }\
diff --git a/boost/math/special_functions/nonfinite_num_facets.hpp b/boost/math/special_functions/nonfinite_num_facets.hpp
index 84d3f1070a..cd1eba9320 100644
--- a/boost/math/special_functions/nonfinite_num_facets.hpp
+++ b/boost/math/special_functions/nonfinite_num_facets.hpp
@@ -24,6 +24,7 @@
 #include <locale>
 
 #include <boost/version.hpp>
+#include <boost/throw_exception.hpp>
 
 #include <boost/math/special_functions/fpclassify.hpp>
 #include <boost/math/special_functions/sign.hpp>
@@ -103,7 +104,7 @@ namespace boost {
         case FP_INFINITE:
           if(flags_ & trap_infinity)
           {
-            throw std::ios_base::failure("Infinity");
+            BOOST_THROW_EXCEPTION(std::ios_base::failure("Infinity"));
           }
           else if((boost::math::signbit)(val))
           { // negative infinity.
@@ -122,7 +123,7 @@ namespace boost {
         case FP_NAN:
           if(flags_ & trap_nan)
           {
-            throw std::ios_base::failure("NaN");
+            BOOST_THROW_EXCEPTION(std::ios_base::failure("NaN"));
           }
           else if((boost::math::signbit)(val))
           { // negative so "-nan".
diff --git a/boost/math/special_functions/owens_t.hpp b/boost/math/special_functions/owens_t.hpp
index 6de93a4887..7fbd8918c0 100644
--- a/boost/math/special_functions/owens_t.hpp
+++ b/boost/math/special_functions/owens_t.hpp
@@ -61,9 +61,9 @@ namespace boost
          inline unsigned short owens_t_compute_code(const RealType h, const RealType a)
          {
             static const RealType hrange[] =
-            {0.02, 0.06, 0.09, 0.125, 0.26, 0.4,  0.6,  1.6,  1.7,  2.33,  2.4,  3.36, 3.4,  4.8};
+            { 0.02f, 0.06f, 0.09f, 0.125f, 0.26f, 0.4f,  0.6f,  1.6f,  1.7f,  2.33f,  2.4f,  3.36f, 3.4f,  4.8f };
 
-            static const RealType arange[] = {0.025, 0.09, 0.15, 0.36, 0.5, 0.9, 0.99999};
+            static const RealType arange[] = { 0.025f, 0.09f, 0.15f, 0.36f, 0.5f, 0.9f, 0.99999f };
             /*
             original select array from paper:
             1, 1, 2,13,13,13,13,13,13,13,13,16,16,16, 9
@@ -229,17 +229,17 @@ namespace boost
 
             static const RealType c2[] =
             {
-               0.99999999999999987510,
-               -0.99999999999988796462,      0.99999999998290743652,
-               -0.99999999896282500134,      0.99999996660459362918,
-               -0.99999933986272476760,      0.99999125611136965852,
-               -0.99991777624463387686,      0.99942835555870132569,
-               -0.99697311720723000295,      0.98751448037275303682,
-               -0.95915857980572882813,      0.89246305511006708555,
-               -0.76893425990463999675,      0.58893528468484693250,
-               -0.38380345160440256652,      0.20317601701045299653,
-               -0.82813631607004984866E-01,  0.24167984735759576523E-01,
-               -0.44676566663971825242E-02,  0.39141169402373836468E-03
+               static_cast<RealType>(0.99999999999999987510),
+               static_cast<RealType>(-0.99999999999988796462),      static_cast<RealType>(0.99999999998290743652),
+               static_cast<RealType>(-0.99999999896282500134),      static_cast<RealType>(0.99999996660459362918),
+               static_cast<RealType>(-0.99999933986272476760),      static_cast<RealType>(0.99999125611136965852),
+               static_cast<RealType>(-0.99991777624463387686),      static_cast<RealType>(0.99942835555870132569),
+               static_cast<RealType>(-0.99697311720723000295),      static_cast<RealType>(0.98751448037275303682),
+               static_cast<RealType>(-0.95915857980572882813),      static_cast<RealType>(0.89246305511006708555),
+               static_cast<RealType>(-0.76893425990463999675),      static_cast<RealType>(0.58893528468484693250),
+               static_cast<RealType>(-0.38380345160440256652),      static_cast<RealType>(0.20317601701045299653),
+               static_cast<RealType>(-0.82813631607004984866E-01),  static_cast<RealType>(0.24167984735759576523E-01),
+               static_cast<RealType>(-0.44676566663971825242E-02),  static_cast<RealType>(0.39141169402373836468E-03)
             };
 
             const RealType as = a*a;
@@ -402,20 +402,22 @@ namespace boost
              */
 
             const unsigned short m = 13;
-            static const RealType pts[] = {0.35082039676451715489E-02,
-               0.31279042338030753740E-01,  0.85266826283219451090E-01,
-               0.16245071730812277011,      0.25851196049125434828,
-               0.36807553840697533536,      0.48501092905604697475,
-               0.60277514152618576821,      0.71477884217753226516,
-               0.81475510988760098605,      0.89711029755948965867,
-               0.95723808085944261843,      0.99178832974629703586};
-            static const RealType wts[] = { 0.18831438115323502887E-01,
-               0.18567086243977649478E-01,  0.18042093461223385584E-01,
-               0.17263829606398753364E-01,  0.16243219975989856730E-01,
-               0.14994592034116704829E-01,  0.13535474469662088392E-01,
-               0.11886351605820165233E-01,  0.10070377242777431897E-01,
-               0.81130545742299586629E-02,  0.60419009528470238773E-02,
-               0.38862217010742057883E-02,  0.16793031084546090448E-02};
+            static const RealType pts[] = {
+               static_cast<RealType>(0.35082039676451715489E-02),
+               static_cast<RealType>(0.31279042338030753740E-01),  static_cast<RealType>(0.85266826283219451090E-01),
+               static_cast<RealType>(0.16245071730812277011),      static_cast<RealType>(0.25851196049125434828),
+               static_cast<RealType>(0.36807553840697533536),      static_cast<RealType>(0.48501092905604697475),
+               static_cast<RealType>(0.60277514152618576821),      static_cast<RealType>(0.71477884217753226516),
+               static_cast<RealType>(0.81475510988760098605),      static_cast<RealType>(0.89711029755948965867),
+               static_cast<RealType>(0.95723808085944261843),      static_cast<RealType>(0.99178832974629703586) };
+            static const RealType wts[] = { 
+               static_cast<RealType>(0.18831438115323502887E-01),
+               static_cast<RealType>(0.18567086243977649478E-01),  static_cast<RealType>(0.18042093461223385584E-01),
+               static_cast<RealType>(0.17263829606398753364E-01),  static_cast<RealType>(0.16243219975989856730E-01),
+               static_cast<RealType>(0.14994592034116704829E-01),  static_cast<RealType>(0.13535474469662088392E-01),
+               static_cast<RealType>(0.11886351605820165233E-01),  static_cast<RealType>(0.10070377242777431897E-01),
+               static_cast<RealType>(0.81130545742299586629E-02),  static_cast<RealType>(0.60419009528470238773E-02),
+               static_cast<RealType>(0.38862217010742057883E-02),  static_cast<RealType>(0.16793031084546090448E-02) };
 
             const RealType as = a*a;
             const RealType hs = -h*h*boost::math::constants::half<RealType>();
@@ -575,14 +577,18 @@ namespace boost
             // when the last accelerated term was small enough...
             //
             int n;
+#ifndef BOOST_NO_EXCEPTIONS
             try
             {
+#endif
                n = itrunc(T(tools::log_max_value<T>() / 6));
+#ifndef BOOST_NO_EXCEPTIONS
             }
             catch(...)
             {
                n = (std::numeric_limits<int>::max)();
             }
+#endif
             n = (std::min)(n, 1500);
             T d = pow(3 + sqrt(T(8)), n);
             d = (d + 1 / d) / 2;
@@ -689,14 +695,18 @@ namespace boost
             // when the last accelerated term was small enough...
             //
             int n;
+#ifndef BOOST_NO_EXCEPTIONS
             try
             {
+#endif
                n = itrunc(RealType(tools::log_max_value<RealType>() / 6));
+#ifndef BOOST_NO_EXCEPTIONS
             }
             catch(...)
             {
                n = (std::numeric_limits<int>::max)();
             }
+#endif
             n = (std::min)(n, 1500);
             RealType d = pow(3 + sqrt(RealType(8)), n);
             d = (d + 1 / d) / 2;
@@ -858,25 +868,33 @@ namespace boost
             bool have_t1(false), have_t2(false);
             if(ah < 3)
             {
+#ifndef BOOST_NO_EXCEPTIONS
                try
                {
+#endif
                   have_t1 = true;
                   p1 = owens_t_T1_accelerated(h, a, forwarding_policy());
                   if(p1.second < target_precision)
                      return p1.first;
+#ifndef BOOST_NO_EXCEPTIONS
                }
                catch(const boost::math::evaluation_error&){}  // T1 may fail and throw, that's OK
+#endif
             }
             if(ah > 1)
             {
+#ifndef BOOST_NO_EXCEPTIONS
                try
                {
+#endif
                   have_t2 = true;
                   p2 = owens_t_T2_accelerated(h, a, ah, forwarding_policy());
                   if(p2.second < target_precision)
                      return p2.first;
+#ifndef BOOST_NO_EXCEPTIONS
                }
                catch(const boost::math::evaluation_error&){}  // T2 may fail and throw, that's OK
+#endif
             }
             //
             // If we haven't tried T1 yet, do it now - sometimes it succeeds and the number of iterations
@@ -884,14 +902,18 @@ namespace boost
             //
             if(!have_t1)
             {
+#ifndef BOOST_NO_EXCEPTIONS
                try
                {
+#endif
                   have_t1 = true;
                   p1 = owens_t_T1_accelerated(h, a, forwarding_policy());
                   if(p1.second < target_precision)
                      return p1.first;
+#ifndef BOOST_NO_EXCEPTIONS
                }
                catch(const boost::math::evaluation_error&){}  // T1 may fail and throw, that's OK
+#endif
             }
             //
             // If we haven't tried T2 yet, do it now - sometimes it succeeds and the number of iterations
@@ -899,24 +921,32 @@ namespace boost
             //
             if(!have_t2)
             {
+#ifndef BOOST_NO_EXCEPTIONS
                try
                {
+#endif
                   have_t2 = true;
                   p2 = owens_t_T2_accelerated(h, a, ah, forwarding_policy());
                   if(p2.second < target_precision)
                      return p2.first;
+#ifndef BOOST_NO_EXCEPTIONS
                }
                catch(const boost::math::evaluation_error&){}  // T2 may fail and throw, that's OK
+#endif
             }
             //
             // OK, nothing left to do but try the most expensive option which is T4,
             // this is often slow to converge, but when it does converge it tends to
             // be accurate:
+#ifndef BOOST_NO_EXCEPTIONS
             try
             {
+#endif
                return T4_mp(h, a, pol);
+#ifndef BOOST_NO_EXCEPTIONS
             }
             catch(const boost::math::evaluation_error&){}  // T4 may fail and throw, that's OK
+#endif
             //
             // Now look back at the results from T1 and T2 and see if either gave better
             // results than we could get from the 64-bit precision versions.
diff --git a/boost/math/special_functions/powm1.hpp b/boost/math/special_functions/powm1.hpp
index f3af3d6e59..fe2fce35d8 100644
--- a/boost/math/special_functions/powm1.hpp
+++ b/boost/math/special_functions/powm1.hpp
@@ -13,23 +13,40 @@
 #include <boost/math/special_functions/math_fwd.hpp>
 #include <boost/math/special_functions/log1p.hpp>
 #include <boost/math/special_functions/expm1.hpp>
+#include <boost/math/special_functions/trunc.hpp>
 #include <boost/assert.hpp>
 
 namespace boost{ namespace math{ namespace detail{
 
 template <class T, class Policy>
-inline T powm1_imp(const T a, const T z, const Policy& pol)
+inline T powm1_imp(const T x, const T y, const Policy& pol)
 {
    BOOST_MATH_STD_USING
+   static const char* function = "boost::math::powm1<%1%>(%1%, %1%)";
 
-   if((fabs(a) < 1) || (fabs(z) < 1))
+   if (x > 0)
    {
-      T p = log(a) * z;
-      if(fabs(p) < 2)
-         return boost::math::expm1(p, pol);
-      // otherwise fall though:
+      if ((fabs(y * (x - 1)) < 0.5) || (fabs(y) < 0.2))
+      {
+         // We don't have any good/quick approximation for log(x) * y
+         // so just try it and see:
+         T l = y * log(x);
+         if (l < 0.5)
+            return boost::math::expm1(l);
+         if (l > boost::math::tools::log_max_value<T>())
+            return boost::math::policies::raise_overflow_error<T>(function, 0, pol);
+         // fall through....
+      }
    }
-   return pow(a, z) - 1;
+   else
+   {
+      // y had better be an integer:
+      if (boost::math::trunc(y) != y)
+         return boost::math::policies::raise_domain_error<T>(function, "For non-integral exponent, expected base > 0 but got %1%", x, pol);
+      if (boost::math::trunc(y / 2) == y / 2)
+         return powm1_imp(T(-x), y, pol);
+   }
+   return pow(x, y) - 1;
 }
 
 } // detail
diff --git a/boost/math/special_functions/relative_difference.hpp b/boost/math/special_functions/relative_difference.hpp
new file mode 100644
index 0000000000..544feda61b
--- /dev/null
+++ b/boost/math/special_functions/relative_difference.hpp
@@ -0,0 +1,134 @@
+//  (C) Copyright John Maddock 2006, 2015
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_MATH_RELATIVE_ERROR
+#define BOOST_MATH_RELATIVE_ERROR
+
+#include <boost/math/special_functions/fpclassify.hpp>
+#include <boost/math/tools/promotion.hpp>
+#include <boost/math/tools/precision.hpp>
+
+namespace boost{
+   namespace math{
+
+      template <class T, class U>
+      typename boost::math::tools::promote_args<T,U>::type relative_difference(const T& arg_a, const U& arg_b)
+      {
+         typedef typename boost::math::tools::promote_args<T, U>::type result_type;
+         result_type a = arg_a;
+         result_type b = arg_b;
+         BOOST_MATH_STD_USING
+#ifdef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
+         //
+         // If math.h has no long double support we can't rely
+         // on the math functions generating exponents outside
+         // the range of a double:
+         //
+         result_type min_val = (std::max)(
+         tools::min_value<result_type>(),
+         static_cast<result_type>((std::numeric_limits<double>::min)()));
+         result_type max_val = (std::min)(
+            tools::max_value<result_type>(),
+            static_cast<result_type>((std::numeric_limits<double>::max)()));
+#else
+         result_type min_val = tools::min_value<result_type>();
+         result_type max_val = tools::max_value<result_type>();
+#endif
+         // Screen out NaN's first, if either value is a NaN then the distance is "infinite":
+         if((boost::math::isnan)(a) || (boost::math::isnan)(b))
+            return max_val;
+         // Screen out infinites:
+         if(fabs(b) > max_val)
+         {
+            if(fabs(a) > max_val)
+               return (a < 0) == (b < 0) ? 0 : max_val;  // one infinity is as good as another!
+            else
+               return max_val;  // one infinity and one finite value implies infinite difference
+         }
+         else if(fabs(a) > max_val)
+            return max_val;    // one infinity and one finite value implies infinite difference
+
+         //
+         // If the values have different signs, treat as infinite difference:
+         //
+         if(((a < 0) != (b < 0)) && (a != 0) && (b != 0))
+            return max_val;
+         a = fabs(a);
+         b = fabs(b);
+         //
+         // Now deal with zero's, if one value is zero (or denorm) then treat it the same as
+         // min_val for the purposes of the calculation that follows:
+         //
+         if(a < min_val)
+            a = min_val;
+         if(b < min_val)
+            b = min_val;
+
+         return (std::max)(fabs((a - b) / a), fabs((a - b) / b));
+      }
+
+#if defined(macintosh) || defined(__APPLE__) || defined(__APPLE_CC__)
+      template <>
+      inline boost::math::tools::promote_args<double, double>::type relative_difference(const double& arg_a, const double& arg_b)
+      {
+         BOOST_MATH_STD_USING
+         double a = arg_a;
+         double b = arg_b;
+         //
+         // On Mac OS X we evaluate "double" functions at "long double" precision,
+         // but "long double" actually has a very slightly narrower range than "double"!  
+         // Therefore use the range of "long double" as our limits since results outside
+         // that range may have been truncated to 0 or INF:
+         //
+         double min_val = (std::max)((double)tools::min_value<long double>(), tools::min_value<double>());
+         double max_val = (std::min)((double)tools::max_value<long double>(), tools::max_value<double>());
+
+         // Screen out NaN's first, if either value is a NaN then the distance is "infinite":
+         if((boost::math::isnan)(a) || (boost::math::isnan)(b))
+            return max_val;
+         // Screen out infinites:
+         if(fabs(b) > max_val)
+         {
+            if(fabs(a) > max_val)
+               return 0;  // one infinity is as good as another!
+            else
+               return max_val;  // one infinity and one finite value implies infinite difference
+         }
+         else if(fabs(a) > max_val)
+            return max_val;    // one infinity and one finite value implies infinite difference
+
+         //
+         // If the values have different signs, treat as infinite difference:
+         //
+         if(((a < 0) != (b < 0)) && (a != 0) && (b != 0))
+            return max_val;
+         a = fabs(a);
+         b = fabs(b);
+         //
+         // Now deal with zero's, if one value is zero (or denorm) then treat it the same as
+         // min_val for the purposes of the calculation that follows:
+         //
+         if(a < min_val)
+            a = min_val;
+         if(b < min_val)
+            b = min_val;
+
+         return (std::max)(fabs((a - b) / a), fabs((a - b) / b));
+      }
+#endif
+
+      template <class T, class U>
+      inline typename boost::math::tools::promote_args<T, U>::type epsilon_difference(const T& arg_a, const U& arg_b)
+      {
+         typedef typename boost::math::tools::promote_args<T, U>::type result_type;
+         result_type r = relative_difference(arg_a, arg_b);
+         if(tools::max_value<result_type>() * boost::math::tools::epsilon<result_type>() < r)
+            return tools::max_value<result_type>();
+         return r / boost::math::tools::epsilon<result_type>();
+      }
+} // namespace math
+} // namespace boost
+
+#endif
diff --git a/boost/math/special_functions/trigamma.hpp b/boost/math/special_functions/trigamma.hpp
index 6fccb36a3a..17c9174c78 100644
--- a/boost/math/special_functions/trigamma.hpp
+++ b/boost/math/special_functions/trigamma.hpp
@@ -69,22 +69,22 @@ T trigamma_prec(T x, const mpl::int_<53>*, const Policy&)
    // Expected Error Term :                       -6.895e-018
    // Maximum Relative Change in Control Points :  8.497e-004
    static const T P_4_inf[] = {
-      0.68947581948701249e-17L,
-      0.49999999999998975L,
-      1.0177274392923795L,
-      2.498208511343429L,
-      2.1921221359427595L,
-      1.5897035272532764L,
-      0.40154388356961734L,
+      static_cast<T>(0.68947581948701249e-17L),
+      static_cast<T>(0.49999999999998975L),
+      static_cast<T>(1.0177274392923795L),
+      static_cast<T>(2.498208511343429L),
+      static_cast<T>(2.1921221359427595L),
+      static_cast<T>(1.5897035272532764L),
+      static_cast<T>(0.40154388356961734L),
    };
    static const T Q_4_inf[] = {
-      1.0L,
-      1.7021215452463932L,
-      4.4290431747556469L,
-      2.9745631894384922L,
-      2.3013614809773616L,
-      0.28360399799075752L,
-      0.022892987908906897L,
+      static_cast<T>(1.0L),
+      static_cast<T>(1.7021215452463932L),
+      static_cast<T>(4.4290431747556469L),
+      static_cast<T>(2.9745631894384922L),
+      static_cast<T>(2.3013614809773616L),
+      static_cast<T>(0.28360399799075752L),
+      static_cast<T>(0.022892987908906897L),
    };
 
    if(x <= 2)
diff --git a/boost/math/special_functions/ulp.hpp b/boost/math/special_functions/ulp.hpp
new file mode 100644
index 0000000000..3d78a1c7c2
--- /dev/null
+++ b/boost/math/special_functions/ulp.hpp
@@ -0,0 +1,71 @@
+//  (C) Copyright John Maddock 2015.
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_MATH_SPECIAL_ULP_HPP
+#define BOOST_MATH_SPECIAL_ULP_HPP
+
+#ifdef _MSC_VER
+#pragma once
+#endif
+
+#include <boost/math/special_functions/math_fwd.hpp>
+#include <boost/math/policies/error_handling.hpp>
+#include <boost/math/special_functions/fpclassify.hpp>
+#include <boost/math/special_functions/next.hpp>
+
+namespace boost{ namespace math{ namespace detail{
+
+template <class T, class Policy>
+T ulp_imp(const T& val, const Policy& pol)
+{
+   BOOST_MATH_STD_USING
+   int expon;
+   static const char* function = "ulp<%1%>(%1%)";
+
+   int fpclass = (boost::math::fpclassify)(val);
+
+   if(fpclass == (int)FP_NAN)
+   {
+      return policies::raise_domain_error<T>(
+         function,
+         "Argument must be finite, but got %1%", val, pol);
+   }
+   else if((fpclass == (int)FP_INFINITE) || (fabs(val) >= tools::max_value<T>()))
+   {
+      return (val < 0 ? -1 : 1) * policies::raise_overflow_error<T>(function, 0, pol);
+   }
+   else if(fpclass == FP_ZERO)
+      return detail::get_smallest_value<T>();
+   //
+   // This code is almost the same as that for float_next, except for negative integers,
+   // where we preserve the relation ulp(x) == ulp(-x) as does Java:
+   //
+   frexp(fabs(val), &expon);
+   T diff = ldexp(T(1), expon - tools::digits<T>());
+   if(diff == 0)
+      diff = detail::get_smallest_value<T>();
+   return diff;
+}
+
+}
+
+template <class T, class Policy>
+inline typename tools::promote_args<T>::type ulp(const T& val, const Policy& pol)
+{
+   typedef typename tools::promote_args<T>::type result_type;
+   return detail::ulp_imp(static_cast<result_type>(val), pol);
+}
+
+template <class T>
+inline typename tools::promote_args<T>::type ulp(const T& val)
+{
+   return ulp(val, policies::policy<>());
+}
+
+
+}} // namespaces
+
+#endif // BOOST_MATH_SPECIAL_ULP_HPP
+
diff --git a/boost/math/special_functions/zeta.hpp b/boost/math/special_functions/zeta.hpp
index 1ba282f4d1..616bb0cccc 100644
--- a/boost/math/special_functions/zeta.hpp
+++ b/boost/math/special_functions/zeta.hpp
@@ -200,20 +200,20 @@ inline T zeta_imp_prec(T s, T sc, const Policy&, const mpl::int_<53>&)
       // Expected Error Term:                         -2.020e-18
       // Max error found at double precision:         3.994987e-17
       static const T P[6] = {    
-         0.24339294433593750202L,
-         -0.49092470516353571651L,
-         0.0557616214776046784287L,
-         -0.00320912498879085894856L,
-         0.000451534528645796438704L,
-         -0.933241270357061460782e-5L,
+         static_cast<T>(0.24339294433593750202L),
+         static_cast<T>(-0.49092470516353571651L),
+         static_cast<T>(0.0557616214776046784287L),
+         static_cast<T>(-0.00320912498879085894856L),
+         static_cast<T>(0.000451534528645796438704L),
+         static_cast<T>(-0.933241270357061460782e-5L),
         };
       static const T Q[6] = {    
-         1L,
-         -0.279960334310344432495L,
-         0.0419676223309986037706L,
-         -0.00413421406552171059003L,
-         0.00024978985622317935355L,
-         -0.101855788418564031874e-4L,
+         static_cast<T>(1L),
+         static_cast<T>(-0.279960334310344432495L),
+         static_cast<T>(0.0419676223309986037706L),
+         static_cast<T>(-0.00413421406552171059003L),
+         static_cast<T>(0.00024978985622317935355L),
+         static_cast<T>(-0.101855788418564031874e-4L),
       };
       result = tools::evaluate_polynomial(P, sc) / tools::evaluate_polynomial(Q, sc);
       result -= 1.2433929443359375F;
@@ -225,20 +225,20 @@ inline T zeta_imp_prec(T s, T sc, const Policy&, const mpl::int_<53>&)
       // Maximum Deviation Found:        9.007e-20
       // Expected Error Term:            9.007e-20
       static const T P[6] = {    
-         0.577215664901532860516,
-         0.243210646940107164097,
-         0.0417364673988216497593,
-         0.00390252087072843288378,
-         0.000249606367151877175456,
-         0.110108440976732897969e-4,
+         static_cast<T>(0.577215664901532860516L),
+         static_cast<T>(0.243210646940107164097L),
+         static_cast<T>(0.0417364673988216497593L),
+         static_cast<T>(0.00390252087072843288378L),
+         static_cast<T>(0.000249606367151877175456L),
+         static_cast<T>(0.110108440976732897969e-4L),
       };
       static const T Q[6] = {    
-         1,
-         0.295201277126631761737,
-         0.043460910607305495864,
-         0.00434930582085826330659,
-         0.000255784226140488490982,
-         0.10991819782396112081e-4,
+         static_cast<T>(1.0),
+         static_cast<T>(0.295201277126631761737L),
+         static_cast<T>(0.043460910607305495864L),
+         static_cast<T>(0.00434930582085826330659L),
+         static_cast<T>(0.000255784226140488490982L),
+         static_cast<T>(0.10991819782396112081e-4L),
       };
       result = tools::evaluate_polynomial(P, T(-sc)) / tools::evaluate_polynomial(Q, T(-sc));
       result += 1 / (-sc);
@@ -249,21 +249,21 @@ inline T zeta_imp_prec(T s, T sc, const Policy&, const mpl::int_<53>&)
       // Expected Error Term:              -5.946e-22
       static const float Y = 0.6986598968505859375;
       static const T P[6] = {    
-         -0.0537258300023595030676,
-         0.0445163473292365591906,
-         0.0128677673534519952905,
-         0.00097541770457391752726,
-         0.769875101573654070925e-4,
-         0.328032510000383084155e-5,
+         static_cast<T>(-0.0537258300023595030676L),
+         static_cast<T>(0.0445163473292365591906L),
+         static_cast<T>(0.0128677673534519952905L),
+         static_cast<T>(0.00097541770457391752726L),
+         static_cast<T>(0.769875101573654070925e-4L),
+         static_cast<T>(0.328032510000383084155e-5L),
       };
       static const T Q[7] = {    
-         1,
-         0.33383194553034051422,
-         0.0487798431291407621462,
-         0.00479039708573558490716,
-         0.000270776703956336357707,
-         0.106951867532057341359e-4,
-         0.236276623974978646399e-7,
+         1.0f,
+         static_cast<T>(0.33383194553034051422L),
+         static_cast<T>(0.0487798431291407621462L),
+         static_cast<T>(0.00479039708573558490716L),
+         static_cast<T>(0.000270776703956336357707L),
+         static_cast<T>(0.106951867532057341359e-4L),
+         static_cast<T>(0.236276623974978646399e-7L),
       };
       result = tools::evaluate_polynomial(P, T(s - 2)) / tools::evaluate_polynomial(Q, T(s - 2));
       result += Y + 1 / (-sc);
@@ -275,23 +275,23 @@ inline T zeta_imp_prec(T s, T sc, const Policy&, const mpl::int_<53>&)
       // Max error found at double precision:         2.009135e-16
 
       static const T P[6] = {    
-         -2.49710190602259410021,
-         -2.60013301809475665334,
-         -0.939260435377109939261,
-         -0.138448617995741530935,
-         -0.00701721240549802377623,
-         -0.229257310594893932383e-4,
+         static_cast<T>(-2.49710190602259410021L),
+         static_cast<T>(-2.60013301809475665334L),
+         static_cast<T>(-0.939260435377109939261L),
+         static_cast<T>(-0.138448617995741530935L),
+         static_cast<T>(-0.00701721240549802377623L),
+         static_cast<T>(-0.229257310594893932383e-4L),
       };
       static const T Q[9] = {    
-         1,
-         0.706039025937745133628,
-         0.15739599649558626358,
-         0.0106117950976845084417,
-         -0.36910273311764618902e-4,
-         0.493409563927590008943e-5,
-         -0.234055487025287216506e-6,
-         0.718833729365459760664e-8,
-         -0.1129200113474947419e-9,
+         1.0f,
+         static_cast<T>(0.706039025937745133628L),
+         static_cast<T>(0.15739599649558626358L),
+         static_cast<T>(0.0106117950976845084417L),
+         static_cast<T>(-0.36910273311764618902e-4L),
+         static_cast<T>(0.493409563927590008943e-5L),
+         static_cast<T>(-0.234055487025287216506e-6L),
+         static_cast<T>(0.718833729365459760664e-8L),
+         static_cast<T>(-0.1129200113474947419e-9L),
       };
       result = tools::evaluate_polynomial(P, T(s - 4)) / tools::evaluate_polynomial(Q, T(s - 4));
       result = 1 + exp(result);
@@ -302,24 +302,24 @@ inline T zeta_imp_prec(T s, T sc, const Policy&, const mpl::int_<53>&)
       // Expected Error Term:                         7.117e-16
       // Max error found at double precision:         9.387771e-16
       static const T P[7] = {    
-         -4.78558028495135619286,
-         -1.89197364881972536382,
-         -0.211407134874412820099,
-         -0.000189204758260076688518,
-         0.00115140923889178742086,
-         0.639949204213164496988e-4,
-         0.139348932445324888343e-5,
+         static_cast<T>(-4.78558028495135619286L),
+         static_cast<T>(-1.89197364881972536382L),
+         static_cast<T>(-0.211407134874412820099L),
+         static_cast<T>(-0.000189204758260076688518L),
+         static_cast<T>(0.00115140923889178742086L),
+         static_cast<T>(0.639949204213164496988e-4L),
+         static_cast<T>(0.139348932445324888343e-5L),
         };
       static const T Q[9] = {    
-         1,
-         0.244345337378188557777,
-         0.00873370754492288653669,
-         -0.00117592765334434471562,
-         -0.743743682899933180415e-4,
-         -0.21750464515767984778e-5,
-         0.471001264003076486547e-8,
-         -0.833378440625385520576e-10,
-         0.699841545204845636531e-12,
+         1.0f,
+         static_cast<T>(0.244345337378188557777L),
+         static_cast<T>(0.00873370754492288653669L),
+         static_cast<T>(-0.00117592765334434471562L),
+         static_cast<T>(-0.743743682899933180415e-4L),
+         static_cast<T>(-0.21750464515767984778e-5L),
+         static_cast<T>(0.471001264003076486547e-8L),
+         static_cast<T>(-0.833378440625385520576e-10L),
+         static_cast<T>(0.699841545204845636531e-12L),
         };
       result = tools::evaluate_polynomial(P, T(s - 7)) / tools::evaluate_polynomial(Q, T(s - 7));
       result = 1 + exp(result);
@@ -329,24 +329,24 @@ inline T zeta_imp_prec(T s, T sc, const Policy&, const mpl::int_<53>&)
       // Max error in interpolated form:             1.668e-17
       // Max error found at long double precision:   1.669714e-17
       static const T P[8] = {    
-         -10.3948950573308896825,
-         -2.85827219671106697179,
-         -0.347728266539245787271,
-         -0.0251156064655346341766,
-         -0.00119459173416968685689,
-         -0.382529323507967522614e-4,
-         -0.785523633796723466968e-6,
-         -0.821465709095465524192e-8,
+         static_cast<T>(-10.3948950573308896825L),
+         static_cast<T>(-2.85827219671106697179L),
+         static_cast<T>(-0.347728266539245787271L),
+         static_cast<T>(-0.0251156064655346341766L),
+         static_cast<T>(-0.00119459173416968685689L),
+         static_cast<T>(-0.382529323507967522614e-4L),
+         static_cast<T>(-0.785523633796723466968e-6L),
+         static_cast<T>(-0.821465709095465524192e-8L),
       };
       static const T Q[10] = {    
-         1,
-         0.208196333572671890965,
-         0.0195687657317205033485,
-         0.00111079638102485921877,
-         0.408507746266039256231e-4,
-         0.955561123065693483991e-6,
-         0.118507153474022900583e-7,
-         0.222609483627352615142e-14,
+         1.0f,
+         static_cast<T>(0.208196333572671890965L),
+         static_cast<T>(0.0195687657317205033485L),
+         static_cast<T>(0.00111079638102485921877L),
+         static_cast<T>(0.408507746266039256231e-4L),
+         static_cast<T>(0.955561123065693483991e-6L),
+         static_cast<T>(0.118507153474022900583e-7L),
+         static_cast<T>(0.222609483627352615142e-14L),
       };
       result = tools::evaluate_polynomial(P, T(s - 15)) / tools::evaluate_polynomial(Q, T(s - 15));
       result = 1 + exp(result);
@@ -883,14 +883,14 @@ T zeta_imp_odd_integer(int s, const T& sc, const Policy& pol, const mpl::false_&
    if(!is_init)
    {
       is_init = true;
-      for(int k = 0; k < sizeof(results) / sizeof(results[0]); ++k)
+      for(unsigned k = 0; k < sizeof(results) / sizeof(results[0]); ++k)
       {
          T arg = k * 2 + 3;
          T c_arg = 1 - arg;
          results[k] = zeta_polynomial_series(arg, c_arg, pol);
       }
    }
-   int index = (s - 3) / 2;
+   unsigned index = (s - 3) / 2;
    return index >= sizeof(results) / sizeof(results[0]) ? zeta_polynomial_series(T(s), sc, pol): results[index];
 }
 
@@ -915,8 +915,12 @@ T zeta_imp(T s, T sc, const Policy& pol, const Tag& tag)
    //
    if(floor(s) == s)
    {
+#ifndef BOOST_NO_EXCEPTIONS
+      // Without exceptions we expect itrunc to return INT_MAX on overflow
+      // and we fall through anyway.
       try
       {
+#endif
          int v = itrunc(s);
          if(v == s)
          {
@@ -925,12 +929,12 @@ T zeta_imp(T s, T sc, const Policy& pol, const Tag& tag)
                if(((-v) & 1) == 0)
                   return 0;
                int n = (-v + 1) / 2;
-               if(n <= boost::math::max_bernoulli_b2n<T>::value)
+               if(n <= (int)boost::math::max_bernoulli_b2n<T>::value)
                   return T((-v & 1) ? -1 : 1) * boost::math::unchecked_bernoulli_b2n<T>(n) / (1 - v);
             }
             else if((v & 1) == 0)
             {
-               if(((v / 2) <= boost::math::max_bernoulli_b2n<T>::value) && (v <= boost::math::max_factorial<T>::value))
+               if(((v / 2) <= (int)boost::math::max_bernoulli_b2n<T>::value) && (v <= (int)boost::math::max_factorial<T>::value))
                   return T(((v / 2 - 1) & 1) ? -1 : 1) * ldexp(T(1), v - 1) * pow(constants::pi<T, Policy>(), v) *
                      boost::math::unchecked_bernoulli_b2n<T>(v / 2) / boost::math::unchecked_factorial<T>(v);
                return T(((v / 2 - 1) & 1) ? -1 : 1) * ldexp(T(1), v - 1) * pow(constants::pi<T, Policy>(), v) *
@@ -939,9 +943,11 @@ T zeta_imp(T s, T sc, const Policy& pol, const Tag& tag)
             else
                return zeta_imp_odd_integer(v, sc, pol, mpl::bool_<(Tag::value <= 113) && Tag::value>());
          }
+#ifndef BOOST_NO_EXCEPTIONS
       }
       catch(const boost::math::rounding_error&){} // Just fall through, s is too large to round
       catch(const std::overflow_error&){}
+#endif
    }
 
    if(fabs(s) < tools::root_epsilon<T>())