Imported Upstream version 1.60.0

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

10 files changed, 100 insertions, 54 deletions

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{