From d484fa657c65cc8a233e1aeabb72c2e7199c581e Mon Sep 17 00:00:00 2001
From: "Jeremy W. Murphy" <jeremy.william.murphy@gmail.com>
Date: Wed, 20 Apr 2016 22:53:56 +1000
Subject: [PATCH] Remove inferior generalized gcd; include integral size
 promotion.

---
 include/boost/math/tools/polynomial.hpp | 67 +++++++++----------------
 test/test_polynomial.cpp                | 11 ++--
 2 files changed, 30 insertions(+), 48 deletions(-)

diff --git a/include/boost/math/tools/polynomial.hpp b/include/boost/math/tools/polynomial.hpp
index 7e5fea6a9..ba3d929bf 100644
--- a/include/boost/math/tools/polynomial.hpp
+++ b/include/boost/math/tools/polynomial.hpp
@@ -16,12 +16,14 @@
 #include <boost/assert.hpp>
 #include <boost/config.hpp>
 #include <boost/function.hpp>
+#include <boost/integer.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 <boost/math/common_factor_rt.hpp>
+#include <boost/mpl/if.hpp>
+#include <boost/operators.hpp>
 
 #include <vector>
 #include <ostream>
@@ -663,19 +665,33 @@ namespace detail
 
 /* 4.6.1C: Greatest common divisor over a unique factorization domain.
  * 
+ * Although step C3 keeps the coefficients to a "reasonable" size, they are
+ * still potentially several binary orders of magnitude larger than the inputs.
  */
 template <class T>
 typename enable_if_c< std::numeric_limits<T>::is_integer, polynomial<T> >::type
-gcd_ufd(polynomial<T> u, polynomial<T> v)
+gcd_ufd(polynomial<T> const &a, polynomial<T> const &b)
 {
-    BOOST_ASSERT(u);
-    BOOST_ASSERT(v);
+    BOOST_ASSERT(a);
+    BOOST_ASSERT(b);
     typedef typename polynomial<T>::size_type N;
+    // If T is a POD, double the output size, otherwise assume it is multi-prec.
+    // This precludes the use of long as a coefficient type in combination with
+    // gcd.
+    /*
+    typedef typename mpl::if_<is_pod<T>,
+                            typename mpl::if_<is_unsigned<T>, 
+                                    typename uint_t<sizeof(T) * 16>::fast, 
+                                    typename int_t<sizeof(T) * 16>::fast >::type,
+                            T>::type U;
+                            */
+    typedef T U;
     using boost::math::gcd;
     
-    T const d = detail::reduce_to_primitive(u, v);
-    T g = 1, h = 1;
-    polynomial<T> r;
+    polynomial<U> u(a), v(b);
+    U const d = detail::reduce_to_primitive(u, v);
+    U g = 1, h = 1;
+    polynomial<U> r;
     while (true)
     {
         // Pseudo-division.
@@ -685,7 +701,7 @@ gcd_ufd(polynomial<T> u, polynomial<T> v)
             goto gcd_end;
         if (r.degree() == 0)
         {
-            v = polynomial<T>(T(1));
+            v = polynomial<U>(U(1));
             goto gcd_end;
         }
         // Adjust remainder.
@@ -703,41 +719,6 @@ gcd_end:
     return d * primitive_part(v);
 }
 
-
-// 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)
diff --git a/test/test_polynomial.cpp b/test/test_polynomial.cpp
index fe898b5e2..d6cd74cf5 100644
--- a/test/test_polynomial.cpp
+++ b/test/test_polynomial.cpp
@@ -224,14 +224,17 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(test_cont_and_pp, T, integral_test_types)
     BOOST_CHECK_EQUAL(a, content(a) * primitive_part(a));
 }
 
-BOOST_AUTO_TEST_CASE_TEMPLATE(test_gcd_ufd, T, integral_test_types)
+
+typedef boost::mpl::list<short, int> sp_integral_gcd_test_types;
+typedef boost::mpl::joint_view<sp_integral_gcd_test_types, mp_integral_test_types> integral_gcd_test_types;
+
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(test_gcd_ufd, T, sp_integral_gcd_test_types)
 {
     polynomial<T> const a(d8.begin(), d8.end());
     polynomial<T> const b(d6.begin(), d6.end());
     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}};
@@ -241,6 +244,4 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(test_gcd_ufd, T, integral_test_types)
     polynomial<T> const w(i2.begin(), i2.end());
     d = gcd_ufd(u, v);
     BOOST_CHECK_EQUAL(d, w);
-    polynomial<T> d_g = generalized_gcd(u, v);
-    BOOST_CHECK_EQUAL(d, w);
 }