Imported Upstream version 1.61.0

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

9 files changed, 369 insertions, 59 deletions

diff --git a/boost/math/bindings/mpfr.hpp b/boost/math/bindings/mpfr.hpp
index 5edc0ae44c..4e4ebcfe53 100644
--- a/boost/math/bindings/mpfr.hpp
+++ b/boost/math/bindings/mpfr.hpp
@@ -297,7 +297,7 @@ struct promote_arg<__gmp_expr<T,U> >
 };
 
 template<>
-inline int digits<mpfr_class>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpfr_class))
+inline int digits<mpfr_class>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpfr_class)) BOOST_NOEXCEPT
 {
    return mpfr_class::get_dprec();
 }
diff --git a/boost/math/cstdfloat/cstdfloat_complex_std.hpp b/boost/math/cstdfloat/cstdfloat_complex_std.hpp
index f949a9aa46..affd92d2da 100644
--- a/boost/math/cstdfloat/cstdfloat_complex_std.hpp
+++ b/boost/math/cstdfloat/cstdfloat_complex_std.hpp
@@ -201,13 +201,13 @@
     };
 
     // Constructors from built-in complex representation of floating-point types.
-    complex<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>::complex(const complex<float>& f)        : re(BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE( f.real())), im(BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE( f.imag())) { }
-    complex<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>::complex(const complex<double>& d)       : re(BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE( d.real())), im(BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE( d.imag())) { }
-    complex<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>::complex(const complex<long double>& ld) : re(BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE(ld.real())), im(BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE(ld.imag())) { }
+    inline complex<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>::complex(const complex<float>& f)        : re(BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE( f.real())), im(BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE( f.imag())) { }
+    inline complex<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>::complex(const complex<double>& d)       : re(BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE( d.real())), im(BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE( d.imag())) { }
+    inline complex<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>::complex(const complex<long double>& ld) : re(BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE(ld.real())), im(BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE(ld.imag())) { }
   } // namespace std
 
   namespace boost { namespace math { namespace cstdfloat { namespace detail {
-  template<class float_type> std::complex<float_type> multiply_by_i(const std::complex<float_type>& x)
+  template<class float_type> inline std::complex<float_type> multiply_by_i(const std::complex<float_type>& x)
   {
     // Multiply x (in C) by I (the imaginary component), and return the result.
     return std::complex<float_type>(-x.imag(), x.real());
diff --git a/boost/math/special_functions/detail/polygamma.hpp b/boost/math/special_functions/detail/polygamma.hpp
index 0267bd3ac8..20a0292fb4 100644
--- a/boost/math/special_functions/detail/polygamma.hpp
+++ b/boost/math/special_functions/detail/polygamma.hpp
@@ -239,7 +239,7 @@
      if(boost::math::tools::max_value<T>() / scale < sum)
         return boost::math::policies::raise_overflow_error<T>(function, 0, pol);
      sum *= scale;
-     return n & 1 ? sum : -sum;
+     return n & 1 ? sum : T(-sum);
   }
 
   //
diff --git a/boost/math/special_functions/detail/t_distribution_inv.hpp b/boost/math/special_functions/detail/t_distribution_inv.hpp
index e8bd0db28f..ab5a8fbca6 100644
--- a/boost/math/special_functions/detail/t_distribution_inv.hpp
+++ b/boost/math/special_functions/detail/t_distribution_inv.hpp
@@ -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 = static_cast<T>(0.85498797333834849467655443627193L);
+            static const T c = static_cast<T>(0.85498797333834849467655443627193);
             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 0d160ecdf1..7c1eab45c2 100644
--- a/boost/math/special_functions/ellint_2.hpp
+++ b/boost/math/special_functions/ellint_2.hpp
@@ -74,7 +74,7 @@ T ellint_e_imp(T phi, T k, const Policy& pol)
     }
     else if(fabs(k) == 1)
     {
-       return invert ? T(-sin(phi)) : sin(phi);
+       return invert ? T(-sin(phi)) : T(sin(phi));
     }
     else
     {
diff --git a/boost/math/special_functions/ellint_3.hpp b/boost/math/special_functions/ellint_3.hpp
index 9ab0f8317a..b6670a196f 100644
--- a/boost/math/special_functions/ellint_3.hpp
+++ b/boost/math/special_functions/ellint_3.hpp
@@ -133,7 +133,7 @@ T ellint_pi_imp(T v, T phi, T k, T vc, const Policy& pol)
          if((m > 0) && (vc > 0))
             result += m * ellint_pi_imp(v, k, vc, pol);
       }
-      return phi < 0 ? -result : result;
+      return phi < 0 ? T(-result) : result;
    }
    if(k == 0)
    {
diff --git a/boost/math/special_functions/fpclassify.hpp b/boost/math/special_functions/fpclassify.hpp
index 0a4e1ac7aa..d83e111c48 100644
--- a/boost/math/special_functions/fpclassify.hpp
+++ b/boost/math/special_functions/fpclassify.hpp
@@ -81,7 +81,12 @@ is used.
 #include <float.h>
 #endif
 #ifdef BOOST_MATH_USE_FLOAT128
+#ifdef __has_include
+#if  __has_include("quadmath.h")
 #include "quadmath.h"
+#define BOOST_MATH_HAS_QUADMATH_H
+#endif
+#endif
 #endif
 
 #ifdef BOOST_NO_STDC_NAMESPACE
@@ -124,9 +129,18 @@ inline bool is_nan_helper(T, const boost::false_type&)
 {
    return false;
 }
-#ifdef BOOST_MATH_USE_FLOAT128
+#if defined(BOOST_MATH_USE_FLOAT128) 
+#if defined(BOOST_MATH_HAS_QUADMATH_H)
 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); }
+#elif defined(BOOST_GNU_STDLIB) && BOOST_GNU_STDLIB && \
+      _GLIBCXX_USE_C99_MATH && !_GLIBCXX_USE_C99_FP_MACROS_DYNAMIC
+inline bool is_nan_helper(__float128 f, const boost::true_type&) { return std::isnan(static_cast<double>(f)); }
+inline bool is_nan_helper(__float128 f, const boost::false_type&) { return std::isnan(static_cast<double>(f)); }
+#else
+inline bool is_nan_helper(__float128 f, const boost::true_type&) { return ::isnan(static_cast<double>(f)); }
+inline bool is_nan_helper(__float128 f, const boost::false_type&) { return ::isnan(static_cast<double>(f)); }
+#endif
 #endif
 }
 
@@ -519,7 +533,7 @@ inline bool (isinf)(long double x)
    return detail::isinf_impl(static_cast<value_type>(x), method());
 }
 #endif
-#ifdef BOOST_MATH_USE_FLOAT128
+#if defined(BOOST_MATH_USE_FLOAT128) && defined(BOOST_MATH_HAS_QUADMATH_H)
 template<>
 inline bool (isinf)(__float128 x)
 {
@@ -611,7 +625,7 @@ inline bool (isnan)(long double x)
    return detail::isnan_impl(x, method());
 }
 #endif
-#ifdef BOOST_MATH_USE_FLOAT128
+#if defined(BOOST_MATH_USE_FLOAT128) && defined(BOOST_MATH_HAS_QUADMATH_H)
 template<>
 inline bool (isnan)(__float128 x)
 {
diff --git a/boost/math/tools/config.hpp b/boost/math/tools/config.hpp
index ffd0ab43a6..75d29b646e 100644
--- a/boost/math/tools/config.hpp
+++ b/boost/math/tools/config.hpp
@@ -265,18 +265,6 @@
 #  define BOOST_MATH_INT_VALUE_SUFFIX(RV, SUF) RV##SUF
 #endif
 //
-// Test whether to support __float128, if we don't have quadmath.h then this can't currently work:
-//
-#ifndef BOOST_MATH_USE_FLOAT128
-#ifdef __has_include
-#if ! __has_include("quadmath.h")
-#define BOOST_MATH_DISABLE_FLOAT128
-#endif
-#elif !defined(BOOST_ARCH_X86)
-#define BOOST_MATH_DISABLE_FLOAT128
-#endif
-#endif
-//
 // And then the actual configuration:
 //
 #if defined(_GLIBCXX_USE_FLOAT128) && defined(BOOST_GCC) && !defined(__STRICT_ANSI__) \
diff --git a/boost/math/tools/polynomial.hpp b/boost/math/tools/polynomial.hpp
index 9858880a4d..73980d6aa8 100644
--- a/boost/math/tools/polynomial.hpp
+++ b/boost/math/tools/polynomial.hpp
@@ -1,4 +1,7 @@
 //  (C) Copyright John Maddock 2006.
+//  (C) Copyright Jeremy William Murphy 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)
@@ -11,13 +14,20 @@
 #endif
 
 #include <boost/assert.hpp>
+#include <boost/config.hpp>
+#include <boost/function.hpp>
+#include <boost/lambda/lambda.hpp>
 #include <boost/math/tools/rational.hpp>
 #include <boost/math/tools/real_cast.hpp>
 #include <boost/math/special_functions/binomial.hpp>
+#include <boost/operators.hpp>
 
 #include <vector>
 #include <ostream>
 #include <algorithm>
+#ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
+#include <initializer_list>
+#endif
 
 namespace boost{ namespace math{ namespace tools{
 
@@ -98,8 +108,157 @@ T evaluate_chebyshev(const Seq& a, const T& x)
    return a[0] / 2 + yk * x - yk1;
 }
 
+
+template <typename T>
+class polynomial;
+
+namespace detail {
+
+/**
+* Knuth, The Art of Computer Programming: Volume 2, Third edition, 1998
+* Chapter 4.6.1, Algorithm D: Division of polynomials over a field.
+*
+* @tparam  T   Coefficient type, must be not be an integer.
+*
+* Template-parameter T actually must be a field but we don't currently have that
+* subtlety of distinction.
+*/
+template <typename T, typename N>
+BOOST_DEDUCED_TYPENAME disable_if_c<std::numeric_limits<T>::is_integer, void >::type
+division_impl(polynomial<T> &q, polynomial<T> &u, const polynomial<T>& v, N n, N k)
+{
+    q[k] = u[n + k] / v[n];
+    for (N j = n + k; j > k;)
+    {
+        j--;
+        u[j] -= q[k] * v[j - k];
+    }
+}
+
+template <class T, class N>
+T integer_power(T t, N n)
+{
+   switch(n)
+   {
+   case 0:
+      return static_cast<T>(1u);
+   case 1:
+      return t;
+   case 2:
+      return t * t;
+   case 3:
+      return t * t * t;
+   }
+   T result = integer_power(t, n / 2);
+   result *= result;
+   if(n & 1)
+      result *= t;
+   return result;
+}
+
+
+/**
+* Knuth, The Art of Computer Programming: Volume 2, Third edition, 1998
+* Chapter 4.6.1, Algorithm R: Pseudo-division of polynomials.
+*
+* @tparam  T   Coefficient type, must be an integer.
+*
+* Template-parameter T actually must be a unique factorization domain but we
+* don't currently have that subtlety of distinction.
+*/
+template <typename T, typename N>
+BOOST_DEDUCED_TYPENAME enable_if_c<std::numeric_limits<T>::is_integer, void >::type
+division_impl(polynomial<T> &q, polynomial<T> &u, const polynomial<T>& v, N n, N k)
+{
+    q[k] = u[n + k] * integer_power(v[n], k);
+    for (N j = n + k; j > 0;)
+    {
+        j--;
+        u[j] = v[n] * u[j] - (j < k ? T(0) : u[n + k] * v[j - k]);
+    }
+}
+
+
+/**
+ * Knuth, The Art of Computer Programming: Volume 2, Third edition, 1998
+ * Chapter 4.6.1, Algorithm D and R: Main loop.
+ *
+ * @param   u   Dividend.
+ * @param   v   Divisor.
+ */
+template <typename T>
+std::pair< polynomial<T>, polynomial<T> >
+division(polynomial<T> u, const polynomial<T>& v)
+{
+    BOOST_ASSERT(v.size() <= u.size());
+    BOOST_ASSERT(v != zero_element(std::multiplies< polynomial<T> >()));
+    BOOST_ASSERT(u != zero_element(std::multiplies< polynomial<T> >()));
+
+    typedef typename polynomial<T>::size_type N;
+    
+    N const m = u.size() - 1, n = v.size() - 1;
+    N k = m - n;
+    polynomial<T> q;
+    q.data().resize(m - n + 1);
+
+    do
+    {
+        division_impl(q, u, v, n, k);
+    }
+    while (k-- != 0);
+    u.data().resize(n);
+    u.normalize(); // Occasionally, the remainder is zeroes.
+    return std::make_pair(q, u);
+}
+
 template <class T>
-class polynomial
+struct identity
+{
+    T operator()(T const &x) const
+    {
+        return x;
+    }
+};
+
+} // namespace detail
+
+/**
+ * Returns the zero element for multiplication of polynomials.
+ */
+template <class T>
+polynomial<T> zero_element(std::multiplies< polynomial<T> >)
+{
+    return polynomial<T>();
+}
+
+template <class T>
+polynomial<T> identity_element(std::multiplies< polynomial<T> >)
+{
+    return polynomial<T>(T(1));
+}
+
+/* Calculates a / b and a % b, returning the pair (quotient, remainder) together
+ * because the same amount of computation yields both.
+ * This function is not defined for division by zero: user beware.
+ */
+template <typename T>
+std::pair< polynomial<T>, polynomial<T> >
+quotient_remainder(const polynomial<T>& dividend, const polynomial<T>& divisor)
+{
+    BOOST_ASSERT(divisor != zero_element(std::multiplies< polynomial<T> >()));
+    if (dividend.size() < divisor.size())
+        return std::make_pair(zero_element(std::multiplies< polynomial<T> >()), dividend);
+    return detail::division(dividend, divisor);
+}
+
+
+template <class T>
+class polynomial :
+    equality_comparable< polynomial<T>,
+    dividable< polynomial<T>,
+    dividable2< polynomial<T>, T,
+    modable< polynomial<T>,
+    modable2< polynomial<T>, T > > > > >
 {
 public:
    // typedefs:
@@ -108,15 +267,26 @@ public:
 
    // construct:
    polynomial(){}
+
    template <class U>
    polynomial(const U* data, unsigned order)
       : m_data(data, data + order + 1)
    {
+       normalize();
    }
+
+   template <class I>
+   polynomial(I first, I last)
+   : m_data(first, last)
+   {
+       normalize();
+   }
+
    template <class U>
-   polynomial(const U& point)
+   explicit polynomial(const U& point)
    {
-      m_data.push_back(point);
+       if (point != U(0))
+          m_data.push_back(point);
    }
 
    // copy:
@@ -131,10 +301,29 @@ public:
          m_data.push_back(boost::math::tools::real_cast<T>(p[i]));
       }
    }
+   
+#ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
+    polynomial(std::initializer_list<T> l) : polynomial(std::begin(l), std::end(l))
+    {
+    }
+    
+    polynomial&
+    operator=(std::initializer_list<T> l)
+    {
+        m_data.assign(std::begin(l), std::end(l));
+        return *this;
+    }
+#endif
+
 
    // access:
    size_type size()const { return m_data.size(); }
-   size_type degree()const { return m_data.size() - 1; }
+   size_type degree()const
+   {
+       if (size() == 0)
+           throw std::logic_error("degree() is undefined for the zero polynomial.");
+       return m_data.size() - 1;
+    }
    value_type& operator[](size_type i)
    {
       return m_data[i];
@@ -152,67 +341,88 @@ public:
       return polynomial_to_chebyshev(m_data);
    }
 
+   std::vector<T> const& data() const
+   {
+       return m_data;
+   }
+
+   std::vector<T> & data()
+   {
+       return m_data;
+   }
+
    // operators:
    template <class U>
    polynomial& operator +=(const U& value)
    {
-      if(m_data.size() == 0)
-         m_data.push_back(value);
-      else
-      {
-         m_data[0] += value;
-      }
-      return *this;
+       addition(value);
+       normalize();
+       return *this;
    }
+
    template <class U>
    polynomial& operator -=(const U& value)
    {
-      if(m_data.size() == 0)
-         m_data.push_back(-value);
-      else
-      {
-         m_data[0] -= value;
-      }
-      return *this;
+       subtraction(value);
+       normalize();
+       return *this;
    }
+
    template <class U>
    polynomial& operator *=(const U& value)
    {
-      for(size_type i = 0; i < m_data.size(); ++i)
-         m_data[i] *= value;
+      multiplication(value);
+      normalize();
       return *this;
    }
+
+   template <class U>
+   polynomial& operator /=(const U& value)
+   {
+       division(value);
+       normalize();
+       return *this;
+   }
+
+   template <class U>
+   polynomial& operator %=(const U& /*value*/)
+   {
+       // We can always divide by a scalar, so there is no remainder:
+       *this = zero_element(std::multiplies<polynomial>());
+       return *this;
+   }
+
    template <class U>
    polynomial& operator +=(const polynomial<U>& value)
    {
-      size_type s1 = (std::min)(m_data.size(), value.size());
-      for(size_type i = 0; i < s1; ++i)
-         m_data[i] += value[i];
-      for(size_type i = s1; i < value.size(); ++i)
-         m_data.push_back(value[i]);
+      addition(value);
+      normalize();
       return *this;
    }
+
    template <class U>
    polynomial& operator -=(const polynomial<U>& value)
    {
-      size_type s1 = (std::min)(m_data.size(), value.size());
-      for(size_type i = 0; i < s1; ++i)
-         m_data[i] -= value[i];
-      for(size_type i = s1; i < value.size(); ++i)
-         m_data.push_back(-value[i]);
-      return *this;
+       subtraction(value);
+       normalize();
+       return *this;
    }
    template <class U>
    polynomial& operator *=(const polynomial<U>& value)
    {
       // TODO: FIXME: use O(N log(N)) algorithm!!!
-      BOOST_ASSERT(value.size());
+      polynomial const zero = zero_element(std::multiplies<polynomial>());
+      if (value == zero)
+      {
+          *this = zero;
+          return *this;
+      }
       polynomial base(*this);
-      *this *= value[0];
+      this->multiplication(value[0]);
       for(size_type i = 1; i < value.size(); ++i)
       {
          polynomial t(base);
-         t *= value[i];
+         t.multiplication(value[i]);
          size_type s = size() - i;
          for(size_type j = 0; j < s; ++j)
          {
@@ -224,10 +434,94 @@ public:
       return *this;
    }
 
+   template <typename U>
+   polynomial& operator /=(const polynomial<U>& value)
+   {
+       *this = quotient_remainder(*this, value).first;
+       return *this;
+   }
+
+   template <typename U>
+   polynomial& operator %=(const polynomial<U>& value)
+   {
+       *this = quotient_remainder(*this, value).second;
+       return *this;
+   }
+
+   /** Remove zero coefficients 'from the top', that is for which there are no
+    *        non-zero coefficients of higher degree. */
+   void normalize()
+   {
+       using namespace boost::lambda;
+       m_data.erase(std::find_if(m_data.rbegin(), m_data.rend(), _1 != T(0)).base(), m_data.end());
+   }
+
 private:
-   std::vector<T> m_data;
+    template <class U, class R1, class R2>
+    polynomial& addition(const U& value, R1 sign, R2 op)
+    {
+        if(m_data.size() == 0)
+            m_data.push_back(sign(value));
+        else
+            m_data[0] = op(m_data[0], value);
+        return *this;
+    }
+
+    template <class U>
+    polynomial& addition(const U& value)
+    {
+        return addition(value, detail::identity<U>(), std::plus<U>());
+    }
+
+    template <class U>
+    polynomial& subtraction(const U& value)
+    {
+        return addition(value, std::negate<U>(), std::minus<U>());
+    }
+
+    template <class U, class R1, class R2>
+    polynomial& addition(const polynomial<U>& value, R1 sign, R2 op)
+    {
+        size_type s1 = (std::min)(m_data.size(), value.size());
+        for(size_type i = 0; i < s1; ++i)
+            m_data[i] = op(m_data[i], value[i]);
+        for(size_type i = s1; i < value.size(); ++i)
+            m_data.push_back(sign(value[i]));
+        return *this;
+    }
+
+    template <class U>
+    polynomial& addition(const polynomial<U>& value)
+    {
+        return addition(value, detail::identity<U>(), std::plus<U>());
+    }
+
+    template <class U>
+    polynomial& subtraction(const polynomial<U>& value)
+    {
+        return addition(value, std::negate<U>(), std::minus<U>());
+    }
+
+    template <class U>
+    polynomial& multiplication(const U& value)
+    {
+        using namespace boost::lambda;
+        std::transform(m_data.begin(), m_data.end(), m_data.begin(), ret<T>(_1 * value));
+        return *this;
+    }
+
+    template <class U>
+    polynomial& division(const U& value)
+    {
+        using namespace boost::lambda;
+        std::transform(m_data.begin(), m_data.end(), m_data.begin(), ret<T>(_1 / value));
+        return *this;
+    }
+
+    std::vector<T> m_data;
 };
 
+
 template <class T>
 inline polynomial<T> operator + (const polynomial<T>& a, const polynomial<T>& b)
 {
@@ -300,6 +594,20 @@ inline polynomial<T> operator * (const U& a, const polynomial<T>& b)
    return result;
 }
 
+template <class T>
+bool operator == (const polynomial<T> &a, const polynomial<T> &b)
+{
+    return a.data() == b.data();
+}
+
+// Unary minus (negate).
+template <class T>
+polynomial<T> operator - (polynomial<T> a)
+{
+    std::transform(a.data().begin(), a.data().end(), a.data().begin(), std::negate<T>());
+    return a;
+}
+
 template <class charT, class traits, class T>
 inline std::basic_ostream<charT, traits>& operator << (std::basic_ostream<charT, traits>& os, const polynomial<T>& poly)
 {