Polynomial: fixes for Stein gcd

include/boost/math/common_factor_rt.hpp

@@ -36,7 +36,10 @@ namespace boost {
       struct gcd_traits_abs_defaults
       {
          inline static const T& abs(const T& val) { return val; }
+
+         inline static void normalize(T&) {}
       };
+
       template <class T>
       struct gcd_traits_abs_defaults<T, false>
       {
@@ -45,10 +48,32 @@ namespace boost {
             using std::abs;
             return abs(val);
          }
+
+         inline static void normalize(T& val) { val = abs(val); }
+      };
+
+
+      template <typename T, bool = is_floating_point<T>::value || (std::numeric_limits<T>::is_specialized && !std::numeric_limits<T>::is_integer)>
+      struct gcd_traits_modulo_defaults
+      {
+          inline static void modulo(T& a, const T& b)
+          {
+              using std::fmod;
+              a = fmod(a, b);
+          }        
+      };
+
+      template <class T>
+      struct gcd_traits_modulo_defaults<T, false>
+      {
+          inline static void modulo(T& a, const T& b)
+          {
+              a %= b;
+          }
       };
 
       template <class T>
-      struct gcd_traits_defaults : public gcd_traits_abs_defaults<T>
+      struct gcd_traits_defaults : public gcd_traits_abs_defaults<T>, gcd_traits_modulo_defaults<T>
       {
          BOOST_FORCEINLINE static unsigned make_odd(T& val)
          {
@@ -60,11 +85,17 @@ namespace boost {
             }
             return r;
          }
+
          inline static bool less(const T& a, const T& b)
          {
             return a < b;
          }
 
+         inline static void subtract(T& a, const T& b)
+         {
+            a -= b;
+         }
+
          enum method_type
          {
             method_euclid = 0,
@@ -78,25 +109,13 @@ namespace boost {
             boost::has_right_shift_assign<T>::value && boost::has_left_shift_assign<T>::value && boost::has_less<T>::value
             ? method_binary : method_euclid;
       };
+
       //
       // Default gcd_traits just inherits from defaults:
       //
       template <class T>
       struct gcd_traits : public gcd_traits_defaults<T> {};
-      //
-      // Special handling for polynomials:
-      //
-      namespace tools {
-         template <class T>
-         class polynomial;
-      }
-
-      template <class T>
-      struct gcd_traits<boost::math::tools::polynomial<T> > : public gcd_traits_defaults<T>
-      {
-         static const boost::math::tools::polynomial<T>& abs(const boost::math::tools::polynomial<T>& val) { return val; }
-      };
-      //
+
       // Some platforms have fast bitscan operations, that allow us to implement
       // make_odd much more efficiently:
       //
@@ -253,6 +272,23 @@ namespace boost {
       };
 #endif
 
+template <typename T>
+inline 
+typename enable_if_c<std::numeric_limits<T>::is_integer, bool>::type
+odd(const T& x)
+{
+    return bool(x & 1u);
+}
+
+template <typename T>
+inline 
+typename enable_if_c<std::numeric_limits<T>::is_integer, bool>::type
+even(const T& x)
+{
+    return !odd(x);
+}
+
+
 namespace detail
 {
 
@@ -277,7 +313,7 @@ namespace detail
 
       shifts = (std::min)(gcd_traits<T>::make_odd(u), gcd_traits<T>::make_odd(v));
 
-      while(gcd_traits<T>::less(1, v))
+      while(u != v)
       {
          u %= v;
          v -= u;
@@ -290,36 +326,48 @@ namespace detail
          if(gcd_traits<T>::less(u, v))
             swap(u, v);
       }
-      return (v == 1 ? v : u) << shifts;
+      u <<= shifts;
+      return u;
    }
 
     /** Stein gcd (aka 'binary gcd')
      * 
      * From Mathematics to Generic Programming, Alexander Stepanov, Daniel Rose
+     * 
+     * What can we say about the Stein Domain? 
+     * - It must have zero.
+     * - It must have parity: being exclusively even or odd.
+     * - It must have a smallest prime (2 for integers, x for polynomials, etc).
+     * - It must have shift operations that add or remove factors of the smallest prime.
+     * - It does not work for negative integers.
      */
     template <typename SteinDomain>
     SteinDomain Stein_gcd(SteinDomain m, SteinDomain n)
     {
         using std::swap;
-        BOOST_ASSERT(m >= 0);
-        BOOST_ASSERT(n >= 0);
-        if (m == SteinDomain(0))
+        BOOST_ASSERT(m || n);
+        if (!m)
             return n;
-        if (n == SteinDomain(0))
+        if (!n)
             return m;
-        // m > 0 && n > 0
-        int d_m = gcd_traits<SteinDomain>::make_odd(m);
-        int d_n = gcd_traits<SteinDomain>::make_odd(n);
+        unsigned shifts = std::min(gcd_traits<SteinDomain>::make_odd(m), gcd_traits<SteinDomain>::make_odd(n));
         // odd(m) && odd(n)
         while (m != n)
         {
-            if (n > m)
+            if (gcd_traits<SteinDomain>::less(m, n))
                 swap(n, m);
-            m -= n;
+            gcd_traits<SteinDomain>::subtract(m, n);
+            BOOST_ASSERT(even(m));
+            // With polynomials for example, it is possible that m is now zero.
+            if (!m)
+            {
+                n <<= shifts;
+                return n;
+            }
             gcd_traits<SteinDomain>::make_odd(m);
         }
         // m == n
-        m <<= (std::min)(d_m, d_n);
+        m <<= shifts;
         return m;
     }
 
@@ -333,9 +381,9 @@ namespace detail
     inline EuclideanDomain Euclid_gcd(EuclideanDomain a, EuclideanDomain b)
     {
         using std::swap;
-        while (b != EuclideanDomain(0))
+        while (b)
         {
-            a %= b;
+            gcd_traits<EuclideanDomain>::modulo(a, b);
             swap(a, b);
         }
         return a;
@@ -380,7 +428,9 @@ namespace detail
 template <typename Integer>
 inline Integer gcd(Integer const &a, Integer const &b)
 {
-    return detail::optimal_gcd_select(static_cast<Integer>(gcd_traits<Integer>::abs(a)), static_cast<Integer>(gcd_traits<Integer>::abs(b)));
+    Integer result = detail::optimal_gcd_select(static_cast<Integer>(gcd_traits<Integer>::abs(a)), static_cast<Integer>(gcd_traits<Integer>::abs(b)));
+    gcd_traits<Integer>::normalize(result);
+    return result;
 }
 
 template <typename Integer>

include/boost/math/tools/polynomial.hpp

@@ -24,6 +24,7 @@
 #include <boost/math/special_functions/binomial.hpp>
 #include <boost/operators.hpp>
 
+#include <limits>
 #include <vector>
 #include <ostream>
 #include <algorithm>
@@ -277,6 +278,13 @@ class polynomial :
        normalize();
    }
 
+
+   explicit polynomial(std::vector<T> const &data)
+   : m_data(data)
+   {
+       normalize();
+   }
+
    template <class I>
    polynomial(I first, I last)
    : m_data(first, last)
@@ -444,8 +452,7 @@ class polynomial :
    template <typename U>
    polynomial& operator >>=(U const &n)
    {
-       BOOST_ASSERT(n <= m_data.size());
-       m_data.erase(m_data.begin(), m_data.begin() + n);
+       m_data.erase(m_data.begin(), m_data.begin() + std::min(size_type(n), m_data.size()));
        return *this;
    }
 
@@ -491,6 +498,12 @@ class polynomial :
        using namespace boost::lambda;
        m_data.erase(std::find_if(m_data.rbegin(), m_data.rend(), _1 != T(0)).base(), m_data.end());
    }
+
+   // Negate in-place.
+   void negate()
+   {
+       std::transform(data().begin(), data().end(), data().begin(), std::negate<T>());
+   }
 
 private:
     template <class U, class R1, class R2>
@@ -696,6 +709,7 @@ polynomial<T> pow(polynomial<T> base, int exp)
     return result;
 }
 
+
 template <class charT, class traits, class T>
 inline std::basic_ostream<charT, traits>& operator << (std::basic_ostream<charT, traits>& os, const polynomial<T>& poly)
 {
@@ -709,11 +723,23 @@ inline std::basic_ostream<charT, traits>& operator << (std::basic_ostream<charT,
    return os;
 }
 
+
+template <typename T>
+T constant_coefficient(polynomial<T> const &x)
+{
+    return x ? x.data().front() : T(0);
+}
+
+// Trivial but useful convenience function referred to simply as l() in Knuth.
+template <class T>
+T leading_coefficient(polynomial<T> const &x)
+{
+    BOOST_ASSERT(x);
+    return x.data().back();
+}
+
 } // namespace tools
 } // namespace math
 } // namespace boost
 
 #endif // BOOST_MATH_TOOLS_POLYNOMIAL_HPP
-
-
-