From dc81bc8e6e97fc5d1eff6578e69eacb85b47e3ce Mon Sep 17 00:00:00 2001
From: "Jeremy W. Murphy" <jeremy.william.murphy@gmail.com>
Date: Wed, 20 Apr 2016 01:38:05 +1000
Subject: [PATCH] Add 4.6.1E, generalized gcd.

Should be useful for cross-validating 4.6.1C.
---
 include/boost/math/tools/polynomial.hpp | 62 ++++++++++++++++++++++---
 test/test_polynomial.cpp                | 12 +++--
 2 files changed, 64 insertions(+), 10 deletions(-)

diff --git a/include/boost/math/tools/polynomial.hpp b/include/boost/math/tools/polynomial.hpp
index c198814e3..7e5fea6a9 100644
--- a/include/boost/math/tools/polynomial.hpp
+++ b/include/boost/math/tools/polynomial.hpp
@@ -649,23 +649,31 @@ T leading_coefficient(polynomial<T> const &x)
     return x.data().back();
 }
 
+namespace detail
+{
+    template <class T>
+    T reduce_to_primitive(polynomial<T> &u, polynomial<T> &v)
+    {
+        T const u_cont = content(u), v_cont = content(v);
+        u /= u_cont;
+        v /= v_cont;
+        return gcd(u_cont, v_cont);
+    }
+}
+
 /* 4.6.1C: Greatest common divisor over a unique factorization domain.
  * 
  */
 template <class T>
 typename enable_if_c< std::numeric_limits<T>::is_integer, polynomial<T> >::type
-gcd(polynomial<T> u, polynomial<T> v)
+gcd_ufd(polynomial<T> u, polynomial<T> v)
 {
     BOOST_ASSERT(u);
     BOOST_ASSERT(v);
     typedef typename polynomial<T>::size_type N;
     using boost::math::gcd;
     
-    // Reduce to primitive.
-    T const u_cont = content(u), v_cont = content(v);
-    T const d = gcd(u_cont, v_cont);
-    u /= u_cont;
-    v /= v_cont;
+    T const d = detail::reduce_to_primitive(u, v);
     T g = 1, h = 1;
     polynomial<T> r;
     while (true)
@@ -696,6 +704,48 @@ gcd_end:
 }
 
 
+// Knuth, 4.6.1E
+template <class T>
+typename enable_if_c< std::numeric_limits<T>::is_integer, polynomial<T> >::type
+generalized_gcd(polynomial<T> u, polynomial<T> v)
+{
+    BOOST_ASSERT(u);
+    BOOST_ASSERT(v);
+    using boost::math::gcd;
+
+    // Reduce to primitive.
+    T const d = detail::reduce_to_primitive(u, v);
+
+    polynomial<T> r;
+    while (true)
+    {
+        // Pseudo-division.
+        r = u % v;
+        if (!r)
+            goto gcd_end;
+        if (r.degree() == 0)
+        {
+            v = polynomial<T>(T(1));
+            goto gcd_end;
+        }
+        // Make remainder primitive.
+        u = v;
+        v = primitive_part(r);
+    }
+
+gcd_end:
+    return d * v;
+}
+
+
+template <class T>
+typename enable_if_c< std::numeric_limits<T>::is_integer, polynomial<T> >::type
+gcd(polynomial<T> const &u, polynomial<T> const &v)
+{
+    return gcd_ufd(u, v);
+}
+
+
 template <class charT, class traits, class T>
 inline std::basic_ostream<charT, traits>& operator << (std::basic_ostream<charT, traits>& os, const polynomial<T>& poly)
 {
diff --git a/test/test_polynomial.cpp b/test/test_polynomial.cpp
index 29d293bf2..fe898b5e2 100644
--- a/test/test_polynomial.cpp
+++ b/test/test_polynomial.cpp
@@ -228,15 +228,19 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(test_gcd_ufd, T, integral_test_types)
 {
     polynomial<T> const a(d8.begin(), d8.end());
     polynomial<T> const b(d6.begin(), d6.end());
-    polynomial<T> d = gcd(a, b);
+    polynomial<T> d = gcd_ufd(a, b);
     BOOST_CHECK_EQUAL(d, polynomial<T>(1));
-
+    d = generalized_gcd(a, b);
+    BOOST_CHECK_EQUAL(d, polynomial<T>(1));
+    
     boost::array<T, 5> const i4a = {{105, 278, -88, -56, 16}};
     boost::array<T, 5> const i4b = {{70, 232, -44, -64, 16}};
     boost::array<T, 3> const i2 = {{-35, 24, -4}};
     polynomial<T> const u(i4a.begin(), i4a.end());
     polynomial<T> const v(i4b.begin(), i4b.end());
     polynomial<T> const w(i2.begin(), i2.end());
-    d = gcd(u, v);
-    BOOST_CHECK_EQUAL(w, d);
+    d = gcd_ufd(u, v);
+    BOOST_CHECK_EQUAL(d, w);
+    polynomial<T> d_g = generalized_gcd(u, v);
+    BOOST_CHECK_EQUAL(d, w);
 }