diff --git a/include/boost/math/common_factor_rt.hpp b/include/boost/math/common_factor_rt.hpp
index acde21d2c..e3a5c1c2a 100644
--- a/include/boost/math/common_factor_rt.hpp
+++ b/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>
diff --git a/include/boost/math/tools/polynomial.hpp b/include/boost/math/tools/polynomial.hpp
index 464c334d1..73210fb48 100644
--- a/include/boost/math/tools/polynomial.hpp
+++ b/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>
@@ -444,8 +445,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 +491,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 +702,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 +716,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
-
-
-
diff --git a/include/boost/math/tools/polynomial_gcd.hpp b/include/boost/math/tools/polynomial_gcd.hpp
new file mode 100644
index 000000000..7ae602c47
--- /dev/null
+++ b/include/boost/math/tools/polynomial_gcd.hpp
@@ -0,0 +1,153 @@
+//  (C) Copyright John Maddock and Jeremy William Murphy 2016.
+
+//  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_TOOLS_POLYNOMIAL_GCD_HPP
+#define BOOST_MATH_TOOLS_POLYNOMIAL_GCD_HPP
+
+#ifdef _MSC_VER
+#pragma once
+#endif
+
+#include <boost/math/tools/polynomial.hpp>
+#include <boost/math/common_factor_rt.hpp>
+
+namespace boost { namespace math {
+
+// Common gcd traits for polynomials.
+template <typename T>
+struct gcd_traits_polynomial_defaults : public gcd_traits_defaults< boost::math::tools::polynomial<T> >
+{
+    typedef boost::math::tools::polynomial<T> polynomial_type;
+
+    static polynomial_type
+    abs(const polynomial_type &val)
+    {
+        return leading_coefficient(val) < T(0) ? -val : val;
+    }
+    
+    inline static int make_odd(polynomial_type &x)
+    {
+        unsigned r = 0;
+        while (even(x))
+        {
+            x >>= 1;
+            r++;
+        }
+        return r;
+    }
+    
+    inline static bool
+    less(polynomial_type const &a, polynomial_type const &b)
+    {
+        return a.size() < b.size();
+    }
+
+    inline static 
+    void normalize(polynomial_type &x)
+    {
+        using boost::lambda::_1;
+        
+        // This does assume that the coefficients are totally ordered.
+        if (x)
+        {
+            if (leading_coefficient(x) < T(0))
+                x.negate();
+            // Find the first non-zero, because we can't do gcd(0, 0).
+            T const d = gcd_range(find_if(x.data().begin(), x.data().end(), _1 != T(0)), x.data().end()).first;
+            x /= d;
+        }
+    }
+
+    /**
+     * Alexander Stepanov's subtraction normalizes b with the ratio of the 
+     * constant coefficients before subtracting it from a.
+     * 
+     * It does not appear to be reliable due to the floating point arithmetic 
+     * inaccuracy in the division operation.
+     */
+    inline static void
+    Stepanov_subtraction(polynomial_type &a, polynomial_type const &b)
+    {
+        using std::modf;
+        
+        T r = constant_coefficient(a) / constant_coefficient(b);
+        a -= r * b;
+        if (a && modf(a.data().back(), &r))
+            a /= a.data().back();
+    }
+    
+    /**
+     * Joux subtraction multiplies each polynomial by the other's constant
+     * coefficient then subtracts one from the other.
+     * 
+     * Joux, Antoine. Algorithmic cryptanalysis. CRC Press, 2009.
+     */
+    inline static void
+    Joux_subtraction(polynomial_type &a, polynomial_type const &b)
+    {
+        T const a0 = constant_coefficient(a);
+        a *= constant_coefficient(b);
+        a -= a0 * b;        
+    }
+    
+    
+    /**
+     * Factored Joux subtraction first calculates the gcd of the constant 
+     * coefficients and removes that factor before multiplying the polynomials.
+     * The effect is to reduce intermediate value swell.
+     */
+    inline static void
+    factored_Joux_subtraction(polynomial_type &a, polynomial_type const &b)
+    {
+        T const a0 = constant_coefficient(a);
+        T const b0 = constant_coefficient(b);
+        T const gcd_a0b0 = gcd(a0, b0);
+        a *= b0 / gcd_a0b0;
+        a -= a0 / gcd_a0b0 * b;
+    }
+    
+    
+    inline static void
+    subtract(polynomial_type &a, polynomial_type const &b)
+    {
+        gcd_traits< boost::math::tools::polynomial<T> >::Joux_subtraction(a, b);
+    }
+};
+
+
+//
+// Special handling for polynomials:
+// Note that gcd_traits_polynomial is templated on the coefficient type T,
+// not on polynomial<T>.
+//
+template <typename T, typename Enabled = void>
+struct gcd_traits_polynomial 
+{};
+
+
+// Note: Use these trait classes for Z[x] and R[x], etc to customize for operations
+// on different kinds of polynomials.
+
+// gcd_traits for Z[x].
+template <typename T>
+struct gcd_traits_polynomial<T, typename enable_if_c<std::numeric_limits<T>::is_integer>::type> : public gcd_traits_polynomial_defaults<T>
+{};
+
+
+// gcd_traits for R[x].
+template <typename T>
+struct gcd_traits_polynomial<T, typename enable_if_c<is_floating_point<T>::value || (std::numeric_limits<T>::is_specialized && !std::numeric_limits<T>::is_integer) >::type> : public gcd_traits_polynomial_defaults<T>
+{};
+
+
+template <typename T>
+struct gcd_traits< boost::math::tools::polynomial<T> > : public gcd_traits_polynomial<T>
+{};
+
+} // namespace math
+} // namespace boost
+
+#endif // BOOST_MATH_TOOLS_POLYNOMIAL_GCD_HPP
diff --git a/reporting/performance/test_gcd.cpp b/reporting/performance/test_gcd.cpp
index 8e365db99..baf0924f0 100644
--- a/reporting/performance/test_gcd.cpp
+++ b/reporting/performance/test_gcd.cpp
@@ -169,47 +169,7 @@ T binary_textbook(T u, T v)
    return u + v;
 }
 
-//
-// The Mixed Binary Euclid Algorithm
-// Sidi Mohamed Sedjelmaci
-// Electronic Notes in Discrete Mathematics 35 (2009) 169–176
-//
-template <class T>
-T mixed_binary_gcd(T u, T v)
-{
-   using std::swap;
-   if(u < v)
-      swap(u, v);
-
-   unsigned shifts = 0;
-
-   if(!u)
-      return v;
-   if(!v)
-      return u;
-
-   while(even(u) && even(v))
-   {
-      u >>= 1u;
-      v >>= 1u;
-      ++shifts;
-   }
 
-   while(v > 1)
-   {
-      u %= v;
-      v -= u;
-      if(!u)
-         return v << shifts;
-      if(!v)
-         return u << shifts;
-      while(even(u)) u >>= 1u;
-      while(even(v)) v >>= 1u;
-      if(u < v)
-         swap(u, v);
-   }
-   return (v == 1 ? v : u) << shifts;
-}
 
 template <class T>
 void test_type(const char* name)
@@ -220,7 +180,12 @@ void test_type(const char* name)
 
    for(unsigned i = 0; i < 1000; ++i)
    {
-      data.push_back(std::make_pair(get_prime_products<T>(), get_prime_products<T>()));
+       data.push_back(pair<int_type, int_type>());
+       do
+      {     
+        data.back() = std::make_pair(get_prime_products<T>(), get_prime_products<T>());
+      }
+      while (data.back().first == 0 && data.back().second == 0);
    }
    std::string row_name("gcd<");
    row_name += name;
@@ -342,7 +307,8 @@ inline typename enable_if_c<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBit
    using default_ops::eval_lsb;
    using default_ops::eval_is_zero;
    using default_ops::eval_get_sign;
-
+   using boost::math::detail::mixed_binary_gcd;
+   
    if(a.size() == 1)
    {
       eval_gcd(result, b, *a.limbs());
diff --git a/reporting/performance/test_polynomial_gcd.cpp b/reporting/performance/test_polynomial_gcd.cpp
new file mode 100644
index 000000000..7e50dd9f6
--- /dev/null
+++ b/reporting/performance/test_polynomial_gcd.cpp
@@ -0,0 +1,373 @@
+//  Copyright Jeremy Murphy 2017.
+//  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)
+
+#ifdef _MSC_VER
+#  pragma warning (disable : 4224)
+#endif
+
+#include <boost/math/common_factor_rt.hpp>
+#include <boost/math/special_functions/prime.hpp>
+#include <boost/math/tools/polynomial.hpp>
+#include <boost/math/tools/polynomial_gcd.hpp>
+#include <boost/multiprecision/cpp_bin_float.hpp>
+#include <boost/multiprecision/cpp_int.hpp>
+#include <boost/multiprecision/integer.hpp>
+#include <boost/random.hpp>
+#include <boost/array.hpp>
+#include <boost/type_traits.hpp>
+
+#include <iostream>
+#include <algorithm>
+#include <numeric>
+#include <string>
+#include <tuple>
+#include <type_traits>
+#include <vector>
+#include <functional>
+#include <typeinfo>
+
+#include "fibonacci.hpp"
+#include "../../test/table_type.hpp"
+#include "table_helper.hpp"
+#include "performance.hpp"
+
+
+using namespace std;
+using namespace boost::math::tools;
+using boost::math::gcd_traits;
+
+polynomial<double> total_sum(0);
+polynomial<long int> total_sum_int(0);
+typedef boost::multiprecision::cpp_bin_float_quad mp_float_type;
+polynomial<mp_float_type> total_sum_mp(0);
+
+
+template <typename T>
+typename boost::enable_if_c<boost::is_floating_point<T>::value, void>::type
+add_total(polynomial<T> const &x)
+{
+    total_sum += x;
+}
+
+
+template <typename T>
+typename boost::enable_if_c<boost::is_integral<T>::value, void>::type
+add_total(polynomial<T> const &x)
+{
+    total_sum += x;
+}
+
+void add_total(polynomial<mp_float_type> const &x)
+{
+    total_sum_mp += x;
+}
+
+template <typename SteinDomain>
+SteinDomain Stein_gcd_Joux(SteinDomain m, SteinDomain n)
+{
+    using std::swap;
+    BOOST_ASSERT(m || n);
+    if (!m)
+        return n;
+    if (!n)
+        return m;
+    unsigned shifts = std::min(gcd_traits<SteinDomain>::make_odd(m), gcd_traits<SteinDomain>::make_odd(n));
+    // odd(m) && odd(n)
+    while (m != n)
+    {
+        if (gcd_traits<SteinDomain>::less(m, n))
+            swap(n, m);
+        gcd_traits<SteinDomain>::Joux_subtraction(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 <<= shifts;
+    return m;
+}
+
+
+template <typename SteinDomain>
+SteinDomain Stein_gcd_factored_Joux(SteinDomain m, SteinDomain n)
+{
+    using std::swap;
+    BOOST_ASSERT(m || n);
+    if (!m)
+        return n;
+    if (!n)
+        return m;
+    unsigned shifts = std::min(gcd_traits<SteinDomain>::make_odd(m), gcd_traits<SteinDomain>::make_odd(n));
+    // odd(m) && odd(n)
+    while (m != n)
+    {
+        if (gcd_traits<SteinDomain>::less(m, n))
+            swap(n, m);
+        gcd_traits<SteinDomain>::factored_Joux_subtraction(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 <<= shifts;
+    return m;
+}
+
+
+template <typename Func, class Table>
+double exec_timed_test_foo(Func f, const Table& data, double min_elapsed = 0.5)
+{
+    double t = 0;
+    unsigned repeats = 1;
+    typename Table::value_type::first_type sum{0};
+    stopwatch<boost::chrono::high_resolution_clock> w;
+    do
+    {
+       for(unsigned count = 0; count < repeats; ++count)
+       {
+          for(typename Table::size_type n = 0; n < data.size(); ++n)
+            sum += f(data[n].first, data[n].second);
+       }
+
+        t = boost::chrono::duration_cast<boost::chrono::duration<double>>(w.elapsed()).count();
+        if(t < min_elapsed)
+            repeats *= 2;
+    }
+    while(t < min_elapsed);
+    add_total(sum);
+    return t / repeats;
+}
+
+
+template <typename T>
+struct test_function_template
+{
+    vector<pair<T, T> > const & data;
+    const char* data_name;
+    
+    test_function_template(vector<pair<T, T> > const &data, const char* name) : data(data), data_name(name)
+    {
+        cerr << "Testing: " << name << endl;
+    }
+    
+    template <typename Function>
+    void operator()(pair<Function, string> const &f) const
+    {
+        cerr << "Algorithm: " << f.second << endl;
+        auto result = exec_timed_test_foo(f.first, data);
+        auto table_name = string("polynomial gcd method comparison with ") + compiler_name() + string(" on ") + platform_name();
+
+        report_execution_time(result, 
+                            table_name,
+                            string(data_name), 
+                            string(f.second) + "\n" + boost_name());
+    }
+};
+
+boost::random::mt19937 rng;
+boost::random::uniform_int_distribution<> d_0_6(0, 6);
+boost::random::uniform_int_distribution<> d_1_20(1, 20);
+
+template <class T>
+T get_prime_products()
+{
+   int n_primes = d_0_6(rng);
+   switch(n_primes)
+   {
+   case 0:
+      // Generate a power of 2:
+      // return static_cast<T>(1u) << d_1_20(rng);
+      return std::pow(2, d_1_20(rng));
+   case 1:
+      // prime number:
+      return boost::math::prime(d_1_20(rng) + 3);
+   }
+   T result = 1;
+   for(int i = 0; i < n_primes; ++i)
+      result *= boost::math::prime(d_1_20(rng) + 3) * boost::math::prime(d_1_20(rng) + 3) * boost::math::prime(d_1_20(rng) + 3) * boost::math::prime(d_1_20(rng) + 3) * boost::math::prime(d_1_20(rng) + 3);
+   return result;
+}
+
+
+// T is integral.
+template <class T>
+typename boost::enable_if_c<boost::is_integral<T>::value, T>::type
+get_uniform_random()
+{
+    static boost::random::uniform_int_distribution<T> minmax(0, 255);
+    return minmax(rng);
+}
+
+
+// T is floating point, limit to small values so that IVS does not crash test.
+template <class T>
+typename boost::enable_if_c<boost::is_floating_point<T>::value, T>::type
+get_uniform_random()
+{
+    using std::round;
+    static boost::random::uniform_real_distribution<T> minmax(0, 3);
+    return round(minmax(rng));
+}
+
+
+// T is not a POD.
+template <class T>
+typename boost::enable_if_c<!boost::is_pod<T>::value, T>::type
+get_uniform_random()
+{
+    using std::round;
+    static boost::random::uniform_real_distribution<T> minmax(0, 3);
+    return round(minmax(rng));
+}
+
+
+template <class T>
+inline bool even(T const& val)
+{
+   return !(val & 1u);
+}
+
+
+template <class T, class OutputIterator>
+typename boost::enable_if_c<boost::is_integral<T>::value, void>::type
+gcd_algorithms(OutputIterator output)
+{
+    *output++ = std::make_pair(boost::math::detail::Stein_gcd< polynomial<T> >, "Stein_gcd");
+}
+
+template <class T, class OutputIterator>
+typename boost::enable_if_c<boost::is_same<T, boost::multiprecision::cpp_int>::value, void>::type
+gcd_algorithms(OutputIterator output)
+{
+    *output++ = std::make_pair(boost::math::detail::Stein_gcd< polynomial<T> >, "Stein_gcd");
+    // *output++ = std::make_pair(boost::math::detail::subresultant_gcd< polynomial<T> >, "subresultant gcd");    
+}
+
+
+template <class T, class OutputIterator>
+typename boost::enable_if_c<boost::is_floating_point<T>::value || !std::numeric_limits<T>::is_integer, void>::type
+gcd_algorithms(OutputIterator output)
+{
+    // *output++ = std::make_pair(boost::math::detail::Euclid_gcd< polynomial<T> >, "Euclid_gcd");
+    *output++ = std::make_pair(boost::math::detail::Stein_gcd< polynomial<T> >, "Stein_gcd (Stepanov-factored Joux)");
+    *output++ = std::make_pair(Stein_gcd_factored_Joux< polynomial<T> >, "Stein_gcd (factored Joux)");
+    *output++ = std::make_pair(Stein_gcd_Joux< polynomial<T> >, "Stein_gcd (Joux)");
+}
+
+
+template <class T>
+void test_type(const std::string name)
+{
+   using namespace boost::math::detail;
+   std::vector<pair< polynomial<T>, polynomial<T> > > data;
+
+   for(unsigned i = 0; i < 10; ++i)
+   {
+       data.push_back(pair< polynomial<T>, polynomial<T> >());
+       for (unsigned j = 0; j != 5; j++)
+       {
+           data.back().first.data().push_back(get_uniform_random<T>());
+           data.back().second.data().push_back(get_uniform_random<T>());
+       }
+       data.back().first.normalize();
+       data.back().second.normalize();
+   }
+   std::string row_name("gcd<");
+   row_name += name;
+   row_name += "> (random prime number products)";
+
+   typedef pair< function<polynomial<T>(polynomial<T>, polynomial<T>)>, string> f_test;
+   vector<f_test> test_functions;
+   gcd_algorithms<T>(back_inserter(test_functions));
+   for_each(begin(test_functions), end(test_functions), test_function_template< polynomial<T> >(data, row_name.c_str()));
+
+   data.clear();
+#if 0
+   for(unsigned i = 0; i < 1000; ++i)
+   {
+      data.push_back(std::make_pair(get_uniform_random<T>(), get_uniform_random<T>()));
+   }
+   row_name.erase();
+   row_name += "gcd<";
+   row_name += name;
+   row_name += "> (uniform random numbers)";
+   for_each(begin(test_functions), end(test_functions), test_function_template<int_type>(data, row_name.c_str()));
+
+   // Fibonacci number tests:
+   row_name.erase();
+   row_name += "gcd<";
+   row_name += name;
+   row_name += "> (adjacent Fibonacci numbers)";
+   for_each(begin(test_functions), end(test_functions), test_function_template<int_type>(fibonacci_numbers_permution_1<T>(), row_name.c_str()));
+
+   row_name.erase();
+   row_name += "gcd<";
+   row_name += name;
+   row_name += "> (permutations of Fibonacci numbers)";
+   for_each(begin(test_functions), end(test_functions), test_function_template<int_type>(fibonacci_numbers_permution_2<T>(), row_name.c_str()));
+
+   row_name.erase();
+   row_name += "gcd<";
+   row_name += name;
+   row_name += "> (Trivial cases)";
+   for_each(begin(test_functions), end(test_functions), test_function_template<int_type>(trivial_gcd_test_cases<T>(), row_name.c_str()));
+#endif
+}
+
+/*******************************************************************************************************************/
+
+template <class T>
+T generate_random(unsigned bits_wanted)
+{
+   static boost::random::mt19937 gen;
+   typedef boost::random::mt19937::result_type random_type;
+
+   T max_val;
+   unsigned digits;
+   if(std::numeric_limits<T>::is_bounded && (bits_wanted == (unsigned)std::numeric_limits<T>::digits))
+   {
+      max_val = (std::numeric_limits<T>::max)();
+      digits = std::numeric_limits<T>::digits;
+   }
+   else
+   {
+      max_val = T(1) << bits_wanted;
+      digits = bits_wanted;
+   }
+
+   unsigned bits_per_r_val = std::numeric_limits<random_type>::digits - 1;
+   while((random_type(1) << bits_per_r_val) > (gen.max)()) --bits_per_r_val;
+
+   unsigned terms_needed = digits / bits_per_r_val + 1;
+
+   T val = 0;
+   for(unsigned i = 0; i < terms_needed; ++i)
+   {
+      val *= (gen.max)();
+      val += gen();
+   }
+   val %= max_val;
+   return val;
+}
+
+
+
+int main()
+{
+    // test_type<int>("polynomial<int>");
+    test_type<float>("polynomial<float>");
+    test_type<double>("polynomial<double>");
+    test_type<mp_float_type>("polynomial<cpp_bin_float_quad>");
+}
diff --git a/test/test_polynomial.cpp b/test/test_polynomial.cpp
index 892991400..7029e2353 100644
--- a/test/test_polynomial.cpp
+++ b/test/test_polynomial.cpp
@@ -7,6 +7,7 @@
 #define BOOST_TEST_MAIN
 #include <boost/array.hpp>
 #include <boost/math/tools/polynomial.hpp>
+#include <boost/math/tools/polynomial_gcd.hpp>
 #include <boost/math/common_factor_rt.hpp>
 #include <boost/mpl/list.hpp>
 #include <boost/mpl/joint_view.hpp>
@@ -15,10 +16,16 @@
 #include <boost/multiprecision/cpp_int.hpp>
 #include <boost/multiprecision/cpp_bin_float.hpp>
 #include <boost/multiprecision/cpp_dec_float.hpp>
+#include <boost/tuple/tuple.hpp>
+
 #include <utility>
 
+using namespace boost::math;
 using namespace boost::math::tools;
 using namespace std;
+using boost::math::detail::Stein_gcd;
+using boost::math::detail::Euclid_gcd;
+
 
 template <typename T>
 struct answer
@@ -146,38 +153,291 @@ BOOST_AUTO_TEST_CASE( test_division_over_ufd )
 }
 
 
-BOOST_AUTO_TEST_CASE( test_gcd )
-{
-    /* NOTE: Euclidean gcd is not yet customized to return THE greatest 
-     * common polynomial divisor. If d is THE greatest common divisior of u and
-     * v, then gcd(u, v) will return d or -d according to the algorithm.
-     * By convention, it should return d, as for example Maxima and Wolfram 
-     * Alpha do.
-     * This test is an example of the fact that it returns -d.
-     */
-    boost::array<double, 9> const d8 = {{105, 278, -88, -56, 16}};
-    boost::array<double, 7> const d6 = {{70, 232, -44, -64, 16}};
-    boost::array<double, 7> const d2 = {{-35, 24, -4}};
-    polynomial<double> const u(d8.begin(), d8.end());
-    polynomial<double> const v(d6.begin(), d6.end());
-    polynomial<double> const w(d2.begin(), d2.end());
-    polynomial<double> const d = boost::math::gcd(u, v);
-    BOOST_CHECK_EQUAL(w, d);
-}
+typedef boost::mpl::list<int, long> pod_integral_test_types;
 
-// Sanity checks to make sure I didn't break it.
 typedef boost::mpl::list<int, long
 #if !BOOST_WORKAROUND(BOOST_MSVC, <= 1500)
-   , boost::multiprecision::cpp_int
+    , boost::multiprecision::cpp_int
 #endif
 > integral_test_types;
-typedef boost::mpl::list<double
+
+typedef boost::mpl::list<double> pod_floating_point_types;
+typedef boost::mpl::joint_view<pod_integral_test_types, pod_floating_point_types> pod_test_types;
+
+typedef boost::mpl::list<
+#if !BOOST_WORKAROUND(BOOST_MSVC, <= 1500)
+boost::multiprecision::cpp_bin_float_single, boost::multiprecision::cpp_dec_float_50
+#endif
+> mp_floating_point_types;
+
+typedef boost::mpl::list<
 #if !BOOST_WORKAROUND(BOOST_MSVC, <= 1500)
-   , boost::multiprecision::cpp_rational, boost::multiprecision::cpp_bin_float_single, boost::multiprecision::cpp_dec_float_50
+boost::multiprecision::cpp_rational, boost::multiprecision::cpp_bin_float_single, boost::multiprecision::cpp_dec_float_50
 #endif
-> non_integral_test_types;
+> mp_non_integral_test_types;
+
+typedef boost::mpl::joint_view<pod_floating_point_types, mp_floating_point_types> floating_point_types;
+typedef boost::mpl::joint_view<pod_floating_point_types, mp_non_integral_test_types> non_integral_test_types;
 typedef boost::mpl::joint_view<integral_test_types, non_integral_test_types> all_test_types;
 
+// Test data.
+/*
+typedef boost::tuple< boost::array<int, 5>, boost::array<int, 5>, boost::array<int, 3> > test_datum;
+
+std::vector<test_datum> make_test_data()
+{
+    std::vector< test_datum > test_data;
+    boost::array<int, 5> x_data = {{105, 278, -88, -56, 16}};
+    boost::array<int, 5> y_data = {{70, 232, -44, -64, 16}};
+    boost::array<int, 3> z_data = {{35, -24, 4}};
+    test_data.push_back(boost::make_tuple(x_data, y_data, z_data));
+
+    return test_data;
+}
+
+std::vector<test_datum> foo = make_test_data();
+*/
+
+template <typename T>
+struct FM2GP_Ex_8_3__1
+{
+    polynomial<T> x;
+    polynomial<T> y;
+    polynomial<T> z;
+    
+    FM2GP_Ex_8_3__1()
+    {
+        boost::array<T, 5> const x_data = {{105, 278, -88, -56, 16}};
+        boost::array<T, 5> const y_data = {{70, 232, -44, -64, 16}};
+        boost::array<T, 3> const z_data = {{35, -24, 4}};
+        x = polynomial<T>(x_data.begin(), x_data.end());
+        y = polynomial<T>(y_data.begin(), y_data.end());
+        z = polynomial<T>(z_data.begin(), z_data.end());
+    }
+};
+
+template <typename T>
+struct FM2GP_Ex_8_3__2
+{
+    polynomial<T> x;
+    polynomial<T> y;
+    polynomial<T> z;
+    
+    FM2GP_Ex_8_3__2()
+    {
+        boost::array<T, 5> const x_data = {{1, -6, -8, 6, 7}};
+        boost::array<T, 5> const y_data = {{1, -5, -2, 15, 11}};
+        boost::array<T, 3> const z_data = {{1, 2, 1}};
+        x = polynomial<T>(x_data.begin(), x_data.end());
+        y = polynomial<T>(y_data.begin(), y_data.end());
+        z = polynomial<T>(z_data.begin(), z_data.end());
+    }
+};
+
+
+template <typename T>
+struct FM2GP_mixed
+{
+    polynomial<T> x;
+    polynomial<T> y;
+    polynomial<T> z;
+    
+    FM2GP_mixed()
+    {
+        boost::array<T, 4> const x_data = {{-2.2, -3.3, 0, 1}};
+        boost::array<T, 3> const y_data = {{-4.4, 0, 1}};
+        boost::array<T, 2> const z_data= {{-2, 1}};
+        x = polynomial<T>(x_data.begin(), x_data.end());
+        y = polynomial<T>(y_data.begin(), y_data.end());
+        z = polynomial<T>(z_data.begin(), z_data.end());
+    }
+};
+
+
+template <typename T>
+struct FM2GP_trivial
+{
+    polynomial<T> x;
+    polynomial<T> y;
+    polynomial<T> z;
+    
+    FM2GP_trivial()
+    {
+        boost::array<T, 4> const x_data = {{-2, -3, 0, 1}};
+        boost::array<T, 3> const y_data = {{-4, 0, 1}};
+        boost::array<T, 2> const z_data= {{-2, 1}};
+        x = polynomial<T>(x_data.begin(), x_data.end());
+        y = polynomial<T>(y_data.begin(), y_data.end());
+        z = polynomial<T>(z_data.begin(), z_data.end());
+    }
+};
+
+BOOST_AUTO_TEST_SUITE(Stein_gcd_ufd)
+
+BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_Stein_gcd_1_int, T, integral_test_types, FM2GP_Ex_8_3__1<T> )
+{
+    typedef FM2GP_Ex_8_3__1<T> fixture_type;
+    polynomial<T> w;
+    w = Stein_gcd(fixture_type::x, fixture_type::y);
+    gcd_traits< polynomial<T> >::normalize(w);
+    BOOST_CHECK_EQUAL(w, fixture_type::z);
+    w = Stein_gcd(fixture_type::y, fixture_type::x);
+    gcd_traits< polynomial<T> >::normalize(w);
+    BOOST_CHECK_EQUAL(w, fixture_type::z);
+}
+
+BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_Stein_gcd_2_int, T, integral_test_types, FM2GP_Ex_8_3__2<T> )
+{
+    typedef FM2GP_Ex_8_3__2<T> fixture_type;
+    polynomial<T> w;
+    w = Stein_gcd(fixture_type::x, fixture_type::y);
+    gcd_traits< polynomial<T> >::normalize(w);
+    BOOST_CHECK_EQUAL(w, fixture_type::z);
+    w = Stein_gcd(fixture_type::y, fixture_type::x);
+    gcd_traits< polynomial<T> >::normalize(w);
+    BOOST_CHECK_EQUAL(w, fixture_type::z);
+}
+
+BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_Stein_gcd_trivial, T, integral_test_types, FM2GP_trivial<T> )
+{
+    typedef FM2GP_trivial<T> fixture_type;
+    polynomial<T> w;
+    w = Stein_gcd(fixture_type::x, fixture_type::y);
+    gcd_traits< polynomial<T> >::normalize(w);
+    BOOST_CHECK_EQUAL(w, fixture_type::z);
+    w = Stein_gcd(fixture_type::y, fixture_type::x);
+    gcd_traits< polynomial<T> >::normalize(w);
+    BOOST_CHECK_EQUAL(w, fixture_type::z);
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+
+
+BOOST_AUTO_TEST_SUITE(Stein_gcd_field)
+
+BOOST_FIXTURE_TEST_CASE_TEMPLATE( Ex_8_3__1, T, pod_floating_point_types, FM2GP_Ex_8_3__1<T> )
+{
+    typedef FM2GP_Ex_8_3__1<T> fixture_type;
+    polynomial<T> w;
+    w = Stein_gcd(fixture_type::x, fixture_type::y);
+    gcd_traits< polynomial<T> >::normalize(w);
+    BOOST_CHECK_EQUAL(w, fixture_type::z);
+    w = Stein_gcd(fixture_type::y, fixture_type::x);
+    gcd_traits< polynomial<T> >::normalize(w);
+    BOOST_CHECK_EQUAL(w, fixture_type::z);
+}
+
+// TODO
+// BOOST_AUTO_TEST_CASE_EXPECTED_FAILURES(Ex_8_3__2, 6)
+
+BOOST_FIXTURE_TEST_CASE_TEMPLATE( Ex_8_3__2, T, pod_floating_point_types, FM2GP_Ex_8_3__2<T> )
+{
+    typedef FM2GP_Ex_8_3__2<T> fixture_type;
+    polynomial<T> w;
+    w = Stein_gcd(fixture_type::x, fixture_type::y);
+    gcd_traits< polynomial<T> >::normalize(w);
+    BOOST_CHECK_EQUAL(w, fixture_type::z);
+    w = Stein_gcd(fixture_type::y, fixture_type::x);
+    gcd_traits< polynomial<T> >::normalize(w);
+    BOOST_CHECK_EQUAL(w, fixture_type::z);
+}
+
+
+BOOST_FIXTURE_TEST_CASE_TEMPLATE( trivial, T, pod_floating_point_types, FM2GP_trivial<T> )
+{
+    typedef FM2GP_trivial<T> fixture_type;
+    polynomial<T> w;
+    w = Stein_gcd(fixture_type::x, fixture_type::y);
+    gcd_traits< polynomial<T> >::normalize(w);
+    BOOST_CHECK_EQUAL(w, fixture_type::z);
+    w = Stein_gcd(fixture_type::y, fixture_type::x);
+    gcd_traits< polynomial<T> >::normalize(w);
+    BOOST_CHECK_EQUAL(w, fixture_type::z);
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+
+
+BOOST_AUTO_TEST_SUITE(Euclid_gcd_field)
+
+BOOST_FIXTURE_TEST_CASE_TEMPLATE( Ex_8_3__1, T, floating_point_types, FM2GP_Ex_8_3__1<T> )
+{
+    typedef FM2GP_Ex_8_3__1<T> fixture_type;
+    polynomial<T> w;
+    w = gcd_traits< polynomial<T> >::abs(Euclid_gcd(fixture_type::x, fixture_type::y));
+    BOOST_CHECK_EQUAL(w, fixture_type::z);
+    w = gcd_traits< polynomial<T> >::abs(Euclid_gcd(fixture_type::y, fixture_type::x));
+    BOOST_CHECK_EQUAL(w, fixture_type::z);
+}
+
+
+BOOST_FIXTURE_TEST_CASE_TEMPLATE( Ex_8_3__2, T, floating_point_types, FM2GP_Ex_8_3__2<T> )
+{
+    typedef FM2GP_Ex_8_3__2<T> fixture_type;
+    polynomial<T> w;
+    w = gcd_traits< polynomial<T> >::abs(Euclid_gcd(fixture_type::x, fixture_type::y));
+    BOOST_CHECK_EQUAL(w, fixture_type::z);
+    w = gcd_traits< polynomial<T> >::abs(Euclid_gcd(fixture_type::y, fixture_type::x));
+    BOOST_CHECK_EQUAL(w, fixture_type::z);
+}
+
+
+BOOST_FIXTURE_TEST_CASE_TEMPLATE( trivial, T, floating_point_types, FM2GP_trivial<T> )
+{
+    typedef FM2GP_trivial<T> fixture_type;
+    polynomial<T> w;
+    w = gcd_traits< polynomial<T> >::abs(Euclid_gcd(fixture_type::x, fixture_type::y));
+    BOOST_CHECK_EQUAL(w, fixture_type::z);
+    w = gcd_traits< polynomial<T> >::abs(Euclid_gcd(fixture_type::y, fixture_type::x));
+    BOOST_CHECK_EQUAL(w, fixture_type::z);
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+
+
+BOOST_AUTO_TEST_SUITE(Euclid_gcd_ufd)
+
+BOOST_FIXTURE_TEST_CASE_TEMPLATE( Ex_8_3__1, T, integral_test_types, FM2GP_Ex_8_3__1<T> )
+{
+    typedef FM2GP_Ex_8_3__1<T> fixture_type;
+    polynomial<T> w;
+    w = Euclid_gcd(fixture_type::x, fixture_type::y);
+    gcd_traits< polynomial<T> >::normalize(w);
+    BOOST_CHECK_EQUAL(w, fixture_type::z);
+    w = Euclid_gcd(fixture_type::y, fixture_type::x);
+    gcd_traits< polynomial<T> >::normalize(w);
+    BOOST_CHECK_EQUAL(w, fixture_type::z);
+}
+
+
+BOOST_FIXTURE_TEST_CASE_TEMPLATE( Ex_8_3__2, T, integral_test_types, FM2GP_Ex_8_3__2<T> )
+{
+    typedef FM2GP_Ex_8_3__2<T> fixture_type;
+    polynomial<T> w;
+    w = Euclid_gcd(fixture_type::x, fixture_type::y);
+    gcd_traits< polynomial<T> >::normalize(w);
+    BOOST_CHECK_EQUAL(w, fixture_type::z);
+    w = Euclid_gcd(fixture_type::y, fixture_type::x);
+    gcd_traits< polynomial<T> >::normalize(w);
+    BOOST_CHECK_EQUAL(w, fixture_type::z);
+}
+
+
+BOOST_FIXTURE_TEST_CASE_TEMPLATE( trivial, T, integral_test_types, FM2GP_trivial<T> )
+{
+    typedef FM2GP_trivial<T> fixture_type;
+    polynomial<T> w;
+    w = Euclid_gcd(fixture_type::x, fixture_type::y);
+    gcd_traits< polynomial<T> >::normalize(w);
+    BOOST_CHECK_EQUAL(w, fixture_type::z);
+    w = Euclid_gcd(fixture_type::y, fixture_type::x);
+    gcd_traits< polynomial<T> >::normalize(w);
+    BOOST_CHECK_EQUAL(w, fixture_type::z);
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+
+
 BOOST_AUTO_TEST_CASE_TEMPLATE( test_addition, T, all_test_types )
 {
     polynomial<T> const a(d3a.begin(), d3a.end());