From a1b34b8cdfd489f3e1f7c32d0edb3b31f4c5e792 Mon Sep 17 00:00:00 2001
From: Kolya Matteo <kundor@kundor.org>
Date: Fri, 13 May 2016 07:33:23 -0400
Subject: [PATCH 01/15] Some polynomial doc fixes

---
 doc/internals/polynomial.qbk | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/doc/internals/polynomial.qbk b/doc/internals/polynomial.qbk
index 9f1179694..e7fbb0476 100644
--- a/doc/internals/polynomial.qbk
+++ b/doc/internals/polynomial.qbk
@@ -47,6 +47,10 @@
       template <class U>
       polynomial& operator %=(const U& value);
       template <class U>
+      polynomial& operator >>=(const U& a);
+      template <class U>
+      polynomial& operator <<=(const U& a);
+      template <class U>
       polynomial& operator +=(const polynomial<U>& value);
       template <class U>
       polynomial& operator -=(const polynomial<U>& value);
@@ -75,10 +79,10 @@
    polynomial<T> operator - (const polynomial<T>& a, const U& b);
    template <class T, class U>
    polynomial<T> operator * (const polynomial<T>& a, const U& b);
-   template <class T, class U>
-   polynomial<T> operator / (const polynomial<T>& a, const U& b);
-   template <class T, class U>
-   polynomial<T> operator % (const polynomial<T>& a, const U& b);
+   template <class T>
+   polynomial<T> operator / (const polynomial<T>& a, const T& b);
+   template <class T>
+   polynomial<T> operator % (const polynomial<T>& a, const T& b);
 
    template <class U, class T>
    polynomial<T> operator + (const U& a, const polynomial<T>& b);
@@ -88,15 +92,11 @@
    polynomial<T> operator * (const U& a, const polynomial<T>& b);
 
    template <class T>
-   polynomial<T> operator - (const polynomial<T>& a);
+   polynomial<T> operator - (const polynomial<T>& a); // Unary minus
 
    template <class T>
-   polynomial<T> operator >>= (const U& a);
-   template <class T>
    polynomial<T> operator >> (polynomial<T> const &a, const U& b);
    template <class T>
-   polynomial<T> operator <<= (const U& a);
-   template <class T>
    polynomial<T> operator << (polynomial<T> const &a, const U& b);   
    
    template <class T>

From bf46cd37d8fe26631e6ddaab2b21ae7ecd18adcc Mon Sep 17 00:00:00 2001
From: Kolya Matteo <kundor@kundor.org>
Date: Fri, 13 May 2016 07:34:23 -0400
Subject: [PATCH 02/15] Fix polynomial<T> /= polynomial<U> (call to
 quotient_remainder was ambiguous)

---
 include/boost/math/tools/polynomial.hpp | 17 +++++++++++++++--
 1 file changed, 15 insertions(+), 2 deletions(-)

diff --git a/include/boost/math/tools/polynomial.hpp b/include/boost/math/tools/polynomial.hpp
index ae62e8aea..0ce2a1fcd 100644
--- a/include/boost/math/tools/polynomial.hpp
+++ b/include/boost/math/tools/polynomial.hpp
@@ -427,17 +427,30 @@ class polynomial :
       return *this;
    }
 
+   polynomial& operator /=(const polynomial& value)
+   {
+       *this = quotient_remainder(*this, value).first;
+       return *this;
+   }
+
    template <typename U>
    polynomial& operator /=(const polynomial<U>& value)
    {
-       *this = quotient_remainder(*this, value).first;
+       *this = quotient_remainder(*this, polynomial(value)).first;
        return *this;
    }
 
+   polynomial& operator %=(const polynomial& value)
+   {
+       *this = quotient_remainder(*this, value).second;
+       return *this;
+   }
+
+
    template <typename U>
    polynomial& operator %=(const polynomial<U>& value)
    {
-       *this = quotient_remainder(*this, value).second;
+       *this = quotient_remainder(*this, polynomial(value)).second;
        return *this;
    }
 

From 127181801dd979bda3353523298bc8929e500799 Mon Sep 17 00:00:00 2001
From: Kolya Matteo <kundor@kundor.org>
Date: Sat, 14 May 2016 13:31:12 -0400
Subject: [PATCH 03/15] Use Boost.operators for all poly-poly operations, and
 no poly-scalar operations.

---
 include/boost/math/tools/polynomial.hpp | 48 ++++++++++++++++-----------------
 1 file changed, 24 insertions(+), 24 deletions(-)

diff --git a/include/boost/math/tools/polynomial.hpp b/include/boost/math/tools/polynomial.hpp
index 0ce2a1fcd..fa7f72795 100644
--- a/include/boost/math/tools/polynomial.hpp
+++ b/include/boost/math/tools/polynomial.hpp
@@ -254,12 +254,11 @@ quotient_remainder(const polynomial<T>& dividend, const polynomial<T>& divisor)
 
 
 template <class T>
-class polynomial :
-    equality_comparable< polynomial<T>,
-    dividable< polynomial<T>,
-    dividable2< polynomial<T>, T,
-    modable< polynomial<T>,
-    modable2< polynomial<T>, T > > > > >
+class polynomial : ordered_euclidean_ring_operators< polynomial<T> >
+    /* Provides operators +, -, *, /, % for two polynomials (using member += etc.);
+     * operators >, <=, >= using operator < ; and operator != using operator ==.
+     * Note that polynomials are not actually a euclidean ring when T is not
+     * a field type, in particular, if T is integral. */
 {
 public:
    // typedefs:
@@ -544,51 +543,43 @@ class polynomial :
 };
 
 
-template <class T>
-inline polynomial<T> operator + (const polynomial<T>& a, const polynomial<T>& b)
+template <class T, class U>
+inline polynomial<T> operator + (const polynomial<T>& a, const U& b)
 {
    polynomial<T> result(a);
    result += b;
    return result;
 }
 
-template <class T>
-inline polynomial<T> operator - (const polynomial<T>& a, const polynomial<T>& b)
+template <class T, class U>
+inline polynomial<T> operator - (const polynomial<T>& a, const U& b)
 {
    polynomial<T> result(a);
    result -= b;
    return result;
 }
 
-template <class T>
-inline polynomial<T> operator * (const polynomial<T>& a, const polynomial<T>& b)
-{
-   polynomial<T> result(a);
-   result *= b;
-   return result;
-}
-
 template <class T, class U>
-inline polynomial<T> operator + (const polynomial<T>& a, const U& b)
+inline polynomial<T> operator * (const polynomial<T>& a, const U& b)
 {
    polynomial<T> result(a);
-   result += b;
+   result *= b;
    return result;
 }
 
 template <class T, class U>
-inline polynomial<T> operator - (const polynomial<T>& a, const U& b)
+inline polynomial<T> operator / (const polynomial<T>& a, const U& b)
 {
    polynomial<T> result(a);
-   result -= b;
+   result /= b;
    return result;
 }
 
 template <class T, class U>
-inline polynomial<T> operator * (const polynomial<T>& a, const U& b)
+inline polynomial<T> operator % (const polynomial<T>& a, const U& b)
 {
    polynomial<T> result(a);
-   result *= b;
+   result %= b;
    return result;
 }
 
@@ -622,6 +613,15 @@ bool operator == (const polynomial<T> &a, const polynomial<T> &b)
     return a.data() == b.data();
 }
 
+template <class T>
+bool operator < (const polynomial<T> &a, const polynomial<T> &b)
+{
+    if (a.size() != b.size())
+        return a.size() < b.size();
+    return std::lexicographical_compare(a.data().rbegin(), a.data().rend(),
+                                        b.data().rbegin(), b.data().rend());
+}
+
 template <typename T, typename U>
 polynomial<T> operator >> (const polynomial<T>& a, const U& b)
 {

From 9c7d39e2c928739a44ec3d2b7dc19aacfea37126 Mon Sep 17 00:00:00 2001
From: Kolya Matteo <kundor@kundor.org>
Date: Sat, 14 May 2016 14:04:39 -0400
Subject: [PATCH 04/15] Consistently use `inline`

It shouldn't make any difference for template functions, but it's silly
to only have it sometimes
---
 include/boost/math/tools/polynomial.hpp | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/include/boost/math/tools/polynomial.hpp b/include/boost/math/tools/polynomial.hpp
index fa7f72795..260c731fc 100644
--- a/include/boost/math/tools/polynomial.hpp
+++ b/include/boost/math/tools/polynomial.hpp
@@ -608,13 +608,13 @@ inline polynomial<T> operator * (const U& a, const polynomial<T>& b)
 }
 
 template <class T>
-bool operator == (const polynomial<T> &a, const polynomial<T> &b)
+inline bool operator == (const polynomial<T> &a, const polynomial<T> &b)
 {
     return a.data() == b.data();
 }
 
 template <class T>
-bool operator < (const polynomial<T> &a, const polynomial<T> &b)
+inline bool operator < (const polynomial<T> &a, const polynomial<T> &b)
 {
     if (a.size() != b.size())
         return a.size() < b.size();
@@ -623,7 +623,7 @@ bool operator < (const polynomial<T> &a, const polynomial<T> &b)
 }
 
 template <typename T, typename U>
-polynomial<T> operator >> (const polynomial<T>& a, const U& b)
+inline polynomial<T> operator >> (const polynomial<T>& a, const U& b)
 {
     polynomial<T> result(a);
     result >>= b;
@@ -631,7 +631,7 @@ polynomial<T> operator >> (const polynomial<T>& a, const U& b)
 }
 
 template <typename T, typename U>
-polynomial<T> operator << (const polynomial<T>& a, const U& b)
+inline polynomial<T> operator << (const polynomial<T>& a, const U& b)
 {
     polynomial<T> result(a);
     result <<= b;
@@ -640,20 +640,20 @@ polynomial<T> operator << (const polynomial<T>& a, const U& b)
 
 // Unary minus (negate).
 template <class T>
-polynomial<T> operator - (polynomial<T> a)
+inline polynomial<T> operator - (polynomial<T> a)
 {
     std::transform(a.data().begin(), a.data().end(), a.data().begin(), std::negate<T>());
     return a;
 }
 
 template <class T>
-bool odd(polynomial<T> const &a)
+inline bool odd(polynomial<T> const &a)
 {
     return a.size() > 0 && a[0] != static_cast<T>(0);
 }
 
 template <class T>
-bool even(polynomial<T> const &a)
+inline bool even(polynomial<T> const &a)
 {
     return !odd(a);
 }

From 78d7dc08f82bc15f764c9f1c825f2c6fec52b89f Mon Sep 17 00:00:00 2001
From: Kolya Matteo <kundor@kundor.org>
Date: Sat, 14 May 2016 15:08:27 -0400
Subject: [PATCH 05/15] Prevent binding polynomial<U>s to the scalar operators.

---
 include/boost/math/tools/polynomial.hpp | 22 +++++-----------------
 1 file changed, 5 insertions(+), 17 deletions(-)

diff --git a/include/boost/math/tools/polynomial.hpp b/include/boost/math/tools/polynomial.hpp
index 260c731fc..a567987f8 100644
--- a/include/boost/math/tools/polynomial.hpp
+++ b/include/boost/math/tools/polynomial.hpp
@@ -22,6 +22,7 @@
 #include <boost/math/policies/error_handling.hpp>
 #include <boost/math/special_functions/binomial.hpp>
 #include <boost/operators.hpp>
+#include <boost/type_traits/is_convertible.hpp>
 
 #include <vector>
 #include <ostream>
@@ -378,7 +379,8 @@ class polynomial : ordered_euclidean_ring_operators< polynomial<T> >
    }
 
    template <class U>
-   polynomial& operator /=(const U& value)
+   BOOST_DEDUCED_TYPENAME enable_if<boost::is_convertible<U,T>, polynomial >::type&
+   operator /=(const U& value)
    {
        division(value);
        normalize();
@@ -386,7 +388,8 @@ class polynomial : ordered_euclidean_ring_operators< polynomial<T> >
    }
 
    template <class U>
-   polynomial& operator %=(const U& /*value*/)
+   BOOST_DEDUCED_TYPENAME enable_if<boost::is_convertible<U,T>, polynomial >::type&
+   operator %=(const U& /*value*/)
    {
        // We can always divide by a scalar, so there is no remainder:
        *this = zero_element(std::multiplies<polynomial>());
@@ -432,27 +435,12 @@ class polynomial : ordered_euclidean_ring_operators< polynomial<T> >
        return *this;
    }
 
-   template <typename U>
-   polynomial& operator /=(const polynomial<U>& value)
-   {
-       *this = quotient_remainder(*this, polynomial(value)).first;
-       return *this;
-   }
-
    polynomial& operator %=(const polynomial& value)
    {
        *this = quotient_remainder(*this, value).second;
        return *this;
    }
 
-
-   template <typename U>
-   polynomial& operator %=(const polynomial<U>& value)
-   {
-       *this = quotient_remainder(*this, polynomial(value)).second;
-       return *this;
-   }
-
    template <typename U>
    polynomial& operator >>=(U const &n)
    {

From fb6b8997db77160d8256ccf3def50b175958a387 Mon Sep 17 00:00:00 2001
From: Kolya Matteo <kundor@kundor.org>
Date: Sat, 14 May 2016 21:31:04 -0400
Subject: [PATCH 06/15] For integer types, give a remainder with scalar values

---
 include/boost/math/tools/polynomial.hpp | 27 ++++++++++++++++++++++++---
 1 file changed, 24 insertions(+), 3 deletions(-)

diff --git a/include/boost/math/tools/polynomial.hpp b/include/boost/math/tools/polynomial.hpp
index a567987f8..222dc2567 100644
--- a/include/boost/math/tools/polynomial.hpp
+++ b/include/boost/math/tools/polynomial.hpp
@@ -388,11 +388,23 @@ class polynomial : ordered_euclidean_ring_operators< polynomial<T> >
    }
 
    template <class U>
-   BOOST_DEDUCED_TYPENAME enable_if<boost::is_convertible<U,T>, polynomial >::type&
+   BOOST_DEDUCED_TYPENAME enable_if_c<boost::is_convertible<U,T>::value &&
+                           std::numeric_limits<T>::is_integer, polynomial >::type&
+   operator %=(const U& value)
+   {
+       // In the case that T is integral, this preserves the semantics
+       // p == r*(p/r) + (p % r), for polynomial<T> p and U r.
+       modulus(value);
+       normalize();
+       return *this;
+   }
+
+   template <class U>
+   BOOST_DEDUCED_TYPENAME enable_if_c<boost::is_convertible<U,T>::value &&
+                           !std::numeric_limits<T>::is_integer, polynomial >::type&
    operator %=(const U& /*value*/)
    {
-       // We can always divide by a scalar, so there is no remainder:
-       *this = zero_element(std::multiplies<polynomial>());
+       m_data.clear();
        return *this;
    }
 
@@ -527,6 +539,15 @@ class polynomial : ordered_euclidean_ring_operators< polynomial<T> >
         return *this;
     }
 
+    template <class U>
+    polynomial& modulus(const U& value)
+    {
+       typedef typename std::vector<T>::iterator iter;
+       for (iter it = m_data.begin(); it != m_data.end(); ++it)
+           *it -= T(value * T(*it / value));
+       return *this;
+    }
+    
     std::vector<T> m_data;
 };
 

From 89331d01158df9a16b5d63e23ec3a7dde12a571d Mon Sep 17 00:00:00 2001
From: Kolya Matteo <kundor@kundor.org>
Date: Sat, 14 May 2016 21:31:47 -0400
Subject: [PATCH 07/15] Update polynomial docs

---
 doc/internals/polynomial.qbk | 25 ++++++++++++++++---------
 1 file changed, 16 insertions(+), 9 deletions(-)

diff --git a/doc/internals/polynomial.qbk b/doc/internals/polynomial.qbk
index e7fbb0476..1491c8033 100644
--- a/doc/internals/polynomial.qbk
+++ b/doc/internals/polynomial.qbk
@@ -56,10 +56,8 @@
       polynomial& operator -=(const polynomial<U>& value);
       template <class U>
       polynomial& operator *=(const polynomial<U>& value);
-      template <class U>
-      polynomial& operator /=(const polynomial<U>& value);
-      template <class U>
-      polynomial& operator %=(const polynomial<U>& value);
+      polynomial& operator /=(const polynomial<T>& value);
+      polynomial& operator %=(const polynomial<T>& value);
    };
 
    template <class T>
@@ -79,10 +77,10 @@
    polynomial<T> operator - (const polynomial<T>& a, const U& b);
    template <class T, class U>
    polynomial<T> operator * (const polynomial<T>& a, const U& b);
-   template <class T>
-   polynomial<T> operator / (const polynomial<T>& a, const T& b);
-   template <class T>
-   polynomial<T> operator % (const polynomial<T>& a, const T& b);
+   template <class T, class U>
+   polynomial<T> operator / (const polynomial<T>& a, const U& b);
+   template <class T, class U>
+   polynomial<T> operator % (const polynomial<T>& a, const U& b);
 
    template <class U, class T>
    polynomial<T> operator + (const U& a, const polynomial<T>& b);
@@ -105,12 +103,21 @@
    bool operator != (const polynomial<T> &a, const polynomial<T> &b);
 
    template <class T>
-   polynomial<T> pow(polynomial<T> base, int exp);
+   bool operator < (const polynomial<T> &a, const polynomial<T> &b);
+   template <class T>
+   bool operator <= (const polynomial<T> &a, const polynomial<T> &b);
+   template <class T>
+   bool operator > (const polynomial<T> &a, const polynomial<T> &b);
+   template <class T>
+   bool operator >= (const polynomial<T> &a, const polynomial<T> &b);
 
    template <class charT, class traits, class T>
    std::basic_ostream<charT, traits>& operator <<
       (std::basic_ostream<charT, traits>& os, const polynomial<T>& poly);
 
+   template <class T>
+   polynomial<T> pow(polynomial<T> base, int exp);
+
    template <typename T>
    std::pair< polynomial<T>, polynomial<T> >
    quotient_remainder(const polynomial<T>& a, const polynomial<T>& b);

From a1f74c1d3a98a3a5b4be07e4ba843cfeb49d1a4d Mon Sep 17 00:00:00 2001
From: Kolya Matteo <kundor@kundor.org>
Date: Sat, 14 May 2016 21:41:51 -0400
Subject: [PATCH 08/15] Taking a const reference just to copy it is silly

---
 include/boost/math/tools/polynomial.hpp | 63 ++++++++++++++-------------------
 1 file changed, 27 insertions(+), 36 deletions(-)

diff --git a/include/boost/math/tools/polynomial.hpp b/include/boost/math/tools/polynomial.hpp
index 222dc2567..68dd078a3 100644
--- a/include/boost/math/tools/polynomial.hpp
+++ b/include/boost/math/tools/polynomial.hpp
@@ -553,51 +553,45 @@ class polynomial : ordered_euclidean_ring_operators< polynomial<T> >
 
 
 template <class T, class U>
-inline polynomial<T> operator + (const polynomial<T>& a, const U& b)
+inline polynomial<T> operator + (polynomial<T> a, const U& b)
 {
-   polynomial<T> result(a);
-   result += b;
-   return result;
+   a += b;
+   return a;
 }
 
 template <class T, class U>
-inline polynomial<T> operator - (const polynomial<T>& a, const U& b)
+inline polynomial<T> operator - (polynomial<T> a, const U& b)
 {
-   polynomial<T> result(a);
-   result -= b;
-   return result;
+   a -= b;
+   return a;
 }
 
 template <class T, class U>
-inline polynomial<T> operator * (const polynomial<T>& a, const U& b)
+inline polynomial<T> operator * (polynomial<T> a, const U& b)
 {
-   polynomial<T> result(a);
-   result *= b;
-   return result;
+   a *= b;
+   return a;
 }
 
 template <class T, class U>
-inline polynomial<T> operator / (const polynomial<T>& a, const U& b)
+inline polynomial<T> operator / (polynomial<T> a, const U& b)
 {
-   polynomial<T> result(a);
-   result /= b;
-   return result;
+   a /= b;
+   return a;
 }
 
 template <class T, class U>
-inline polynomial<T> operator % (const polynomial<T>& a, const U& b)
+inline polynomial<T> operator % (polynomial<T> a, const U& b)
 {
-   polynomial<T> result(a);
-   result %= b;
-   return result;
+   a %= b;
+   return a;
 }
 
 template <class U, class T>
-inline polynomial<T> operator + (const U& a, const polynomial<T>& b)
+inline polynomial<T> operator + (const U& a, polynomial<T> b)
 {
-   polynomial<T> result(b);
-   result += a;
-   return result;
+   b += a;
+   return b;
 }
 
 template <class U, class T>
@@ -609,11 +603,10 @@ inline polynomial<T> operator - (const U& a, const polynomial<T>& b)
 }
 
 template <class U, class T>
-inline polynomial<T> operator * (const U& a, const polynomial<T>& b)
+inline polynomial<T> operator * (const U& a, polynomial<T> b)
 {
-   polynomial<T> result(b);
-   result *= a;
-   return result;
+   b *= a;
+   return b;
 }
 
 template <class T>
@@ -632,19 +625,17 @@ inline bool operator < (const polynomial<T> &a, const polynomial<T> &b)
 }
 
 template <typename T, typename U>
-inline polynomial<T> operator >> (const polynomial<T>& a, const U& b)
+inline polynomial<T> operator >> (polynomial<T> a, const U& b)
 {
-    polynomial<T> result(a);
-    result >>= b;
-    return result;
+    a >>= b;
+    return a;
 }
 
 template <typename T, typename U>
-inline polynomial<T> operator << (const polynomial<T>& a, const U& b)
+inline polynomial<T> operator << (polynomial<T> a, const U& b)
 {
-    polynomial<T> result(a);
-    result <<= b;
-    return result;
+    a <<= b;
+    return a;
 }
 
 // Unary minus (negate).

From f4d73feea02edaaffb66516a0f6f69e36f6ec059 Mon Sep 17 00:00:00 2001
From: Kolya Matteo <kundor@kundor.org>
Date: Sun, 15 May 2016 13:46:31 -0400
Subject: [PATCH 09/15] Tests for mixed-type arithmetic, poly<int> % scalar,
 ordering

---
 test/test_polynomial.cpp | 59 ++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 59 insertions(+)

diff --git a/test/test_polynomial.cpp b/test/test_polynomial.cpp
index 65498d21c..28b50a7d6 100644
--- a/test/test_polynomial.cpp
+++ b/test/test_polynomial.cpp
@@ -321,3 +321,62 @@ BOOST_AUTO_TEST_CASE_TEMPLATE( test_pow, T, all_test_types )
     BOOST_CHECK_THROW(pow(a, -1), std::domain_error);
     BOOST_CHECK_EQUAL(pow(one, 137), one);
 }
+
+BOOST_AUTO_TEST_CASE( test_mixed_consistency )
+{
+    polynomial<int> const a(d3a.begin(), d3a.end()); //3x^3 - 4x^2 - 6x + 10
+
+    BOOST_CHECK_EQUAL(a * 0.66667, a / 1.5);
+
+    polynomial<int> b(a);
+    b /= 1.5;
+    BOOST_CHECK_EQUAL(b, a / 1.5);
+
+}
+
+BOOST_AUTO_TEST_CASE( test_int_remainder )
+{
+    polynomial<int> const a(d3a.begin(), d3a.end()); //3x^3 - 4x^2 - 6x + 10
+    boost::array<int, 4> const x3c = {{0,0,0,1}};
+    boost::array<int, 4> const rem3c = {{1, 0, -1, 0}};
+    polynomial<int> x3(x3c.begin(), x3c.end());
+    polynomial<int> rem3(rem3c.begin(), rem3c.end());
+
+    BOOST_CHECK_EQUAL(a, 2*(a / 2) + (a % 2));
+    BOOST_CHECK_EQUAL(a, 3*(a / 3) + (a % 3));
+
+    BOOST_CHECK_EQUAL(a % 2, x3);
+    BOOST_CHECK_EQUAL(a % 3, rem3);
+}
+
+BOOST_AUTO_TEST_CASE( test_ordering )
+{
+    boost::array<double, 3> const g2c = {{-6, 1, 9}};
+    boost::array<polynomial<double>, 7> polys = {{
+      polynomial<double>(d2.begin(), d2.end()),   // 9x^2 - 6
+      polynomial<double>(g2c.begin(), g2c.end()), // 9x^2 + x - 6
+      polynomial<double>(d3b.begin(), d3b.end()), //  x^3 + 6x^2 + 5x - 7
+      polynomial<double>(d3a.begin(), d3a.end()), // 3x^3 - 4x^2 - 6x + 10
+      polynomial<double>(d6.begin(), d6.end()),   // 3x^6 + 5x^4 - 4x^2 - 9x + 21
+      polynomial<double>(d8.begin(), d8.end()), // x^8 + x^6 - 3x^4 - 3x^3 + 8x^2 + 2x - 5
+      polynomial<double>(d8b.begin(), d8b.end()) // x^8 + x^6 - 3x^4 - 3x^3 + 8x^2 + 2x
+    }};
+
+    for (unsigned i = 0; i < polys.size(); ++i)
+    {
+        BOOST_CHECK(!(polys[i] < polys[i]));
+        BOOST_CHECK(polys[i] <= polys[i]);
+        BOOST_CHECK(polys[i] >= polys[i]);
+        BOOST_CHECK(!(polys[i] > polys[i]));
+        for (unsigned j = i + 1; j < polys.size(); ++j)
+        {
+            BOOST_CHECK(polys[i] <= polys[j]);
+            BOOST_CHECK(!(polys[i] >= polys[j]));
+            BOOST_CHECK(polys[i] < polys[j]);
+            BOOST_CHECK(polys[j] > polys[i]);
+            BOOST_CHECK(polys[j] >= polys[i]);
+            BOOST_CHECK(!(polys[j] <= polys[i]));
+        }
+    }
+}
+

From 35022a620c4494182c44ccc3370d19ebb83b855a Mon Sep 17 00:00:00 2001
From: Kolya Matteo <kundor@kundor.org>
Date: Thu, 19 May 2016 10:36:10 -0400
Subject: [PATCH 10/15] Do not template shift operators' right argument

---
 include/boost/math/tools/polynomial.hpp | 17 ++++++++---------
 1 file changed, 8 insertions(+), 9 deletions(-)

diff --git a/include/boost/math/tools/polynomial.hpp b/include/boost/math/tools/polynomial.hpp
index 68dd078a3..b31295a8c 100644
--- a/include/boost/math/tools/polynomial.hpp
+++ b/include/boost/math/tools/polynomial.hpp
@@ -453,17 +453,16 @@ class polynomial : ordered_euclidean_ring_operators< polynomial<T> >
        return *this;
    }
 
-   template <typename U>
-   polynomial& operator >>=(U const &n)
+   polynomial& operator >>=(int n)
    {
-       BOOST_ASSERT(n <= m_data.size());
+       BOOST_ASSERT(n >= 0 && n <= m_data.size());
        m_data.erase(m_data.begin(), m_data.begin() + n);
        return *this;
    }
 
-   template <typename U>
-   polynomial& operator <<=(U const &n)
+   polynomial& operator <<=(int n)
    {
+       BOOST_ASSERT(n >= 0);
        m_data.insert(m_data.begin(), n, static_cast<T>(0));
        normalize();
        return *this;
@@ -624,15 +623,15 @@ inline bool operator < (const polynomial<T> &a, const polynomial<T> &b)
                                         b.data().rbegin(), b.data().rend());
 }
 
-template <typename T, typename U>
-inline polynomial<T> operator >> (polynomial<T> a, const U& b)
+template <typename T>
+inline polynomial<T> operator >> (polynomial<T> a, int b)
 {
     a >>= b;
     return a;
 }
 
-template <typename T, typename U>
-inline polynomial<T> operator << (polynomial<T> a, const U& b)
+template <typename T>
+inline polynomial<T> operator << (polynomial<T> a, int b)
 {
     a <<= b;
     return a;

From 8646ae4b3179fc1928abe63c3ce3d89a37c28c4d Mon Sep 17 00:00:00 2001
From: Kolya Matteo <kundor@kundor.org>
Date: Thu, 19 May 2016 10:58:30 -0400
Subject: [PATCH 11/15] Silence comparison warnings

---
 include/boost/math/tools/polynomial.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/include/boost/math/tools/polynomial.hpp b/include/boost/math/tools/polynomial.hpp
index b31295a8c..3cd99de42 100644
--- a/include/boost/math/tools/polynomial.hpp
+++ b/include/boost/math/tools/polynomial.hpp
@@ -455,7 +455,7 @@ class polynomial : ordered_euclidean_ring_operators< polynomial<T> >
 
    polynomial& operator >>=(int n)
    {
-       BOOST_ASSERT(n >= 0 && n <= m_data.size());
+       BOOST_ASSERT(n >= 0 && static_cast<unsigned>(n) <= m_data.size());
        m_data.erase(m_data.begin(), m_data.begin() + n);
        return *this;
    }

From 910a62a5b386a46347d68674985c924a483d5493 Mon Sep 17 00:00:00 2001
From: Kolya Matteo <kundor@kundor.org>
Date: Thu, 19 May 2016 11:02:44 -0400
Subject: [PATCH 12/15] Remove polynomial<T> x polynomial<U> arithmetic

---
 include/boost/math/tools/polynomial.hpp | 12 +++++-------
 1 file changed, 5 insertions(+), 7 deletions(-)

diff --git a/include/boost/math/tools/polynomial.hpp b/include/boost/math/tools/polynomial.hpp
index 3cd99de42..3d617d22d 100644
--- a/include/boost/math/tools/polynomial.hpp
+++ b/include/boost/math/tools/polynomial.hpp
@@ -295,7 +295,7 @@ class polynomial : ordered_euclidean_ring_operators< polynomial<T> >
       : m_data(p.m_data) { }
 
    template <class U>
-   polynomial(const polynomial<U>& p)
+   explicit polynomial(const polynomial<U>& p)
    {
       for(unsigned i = 0; i < p.size(); ++i)
       {
@@ -408,23 +408,21 @@ class polynomial : ordered_euclidean_ring_operators< polynomial<T> >
        return *this;
    }
 
-   template <class U>
-   polynomial& operator +=(const polynomial<U>& value)
+   polynomial& operator +=(const polynomial& value)
    {
       addition(value);
       normalize();
       return *this;
    }
 
-   template <class U>
-   polynomial& operator -=(const polynomial<U>& value)
+   polynomial& operator -=(const polynomial& value)
    {
        subtraction(value);
        normalize();
        return *this;
    }
-   template <class U>
-   polynomial& operator *=(const polynomial<U>& value)
+
+   polynomial& operator *=(const polynomial& value)
    {
       // TODO: FIXME: use O(N log(N)) algorithm!!!
       polynomial const zero = zero_element(std::multiplies<polynomial>());

From 2798e2389534f3977766425093c1d0eb5aba57e9 Mon Sep 17 00:00:00 2001
From: Kolya Matteo <kundor@kundor.org>
Date: Sat, 21 May 2016 21:48:43 -0400
Subject: [PATCH 13/15] No more need for enable_if guards

We no longer need polynomial<U> arguments to prefer a different overload,
since mixed-type polynomial operations were nixed.
---
 include/boost/math/tools/polynomial.hpp | 27 +++++++++++----------------
 1 file changed, 11 insertions(+), 16 deletions(-)

diff --git a/include/boost/math/tools/polynomial.hpp b/include/boost/math/tools/polynomial.hpp
index 3d617d22d..1c1fe8a5d 100644
--- a/include/boost/math/tools/polynomial.hpp
+++ b/include/boost/math/tools/polynomial.hpp
@@ -379,8 +379,7 @@ class polynomial : ordered_euclidean_ring_operators< polynomial<T> >
    }
 
    template <class U>
-   BOOST_DEDUCED_TYPENAME enable_if<boost::is_convertible<U,T>, polynomial >::type&
-   operator /=(const U& value)
+   polynomial& operator /=(const U& value)
    {
        division(value);
        normalize();
@@ -388,23 +387,19 @@ class polynomial : ordered_euclidean_ring_operators< polynomial<T> >
    }
 
    template <class U>
-   BOOST_DEDUCED_TYPENAME enable_if_c<boost::is_convertible<U,T>::value &&
-                           std::numeric_limits<T>::is_integer, polynomial >::type&
-   operator %=(const U& value)
+   polynomial& operator %=(const U& value)
    {
        // In the case that T is integral, this preserves the semantics
        // p == r*(p/r) + (p % r), for polynomial<T> p and U r.
-       modulus(value);
-       normalize();
-       return *this;
-   }
-
-   template <class U>
-   BOOST_DEDUCED_TYPENAME enable_if_c<boost::is_convertible<U,T>::value &&
-                           !std::numeric_limits<T>::is_integer, polynomial >::type&
-   operator %=(const U& /*value*/)
-   {
-       m_data.clear();
+       if (std::numeric_limits<T>::is_integer)
+       {
+           modulus(value);
+           normalize();
+       }
+       else
+       {
+           m_data.clear();
+       }
        return *this;
    }
 

From 885e13f90c6b18fa3c9810ecb4093b5c5464a2ba Mon Sep 17 00:00:00 2001
From: Kolya Matteo <kundor@kundor.org>
Date: Sun, 22 May 2016 09:25:30 -0400
Subject: [PATCH 14/15] Update polynomial docs

---
 doc/internals/polynomial.qbk | 25 +++++++++++--------------
 1 file changed, 11 insertions(+), 14 deletions(-)

diff --git a/doc/internals/polynomial.qbk b/doc/internals/polynomial.qbk
index 1491c8033..b9f890336 100644
--- a/doc/internals/polynomial.qbk
+++ b/doc/internals/polynomial.qbk
@@ -25,6 +25,8 @@
       polynomial(I first, I last);
       template <class U>
       explicit polynomial(const U& point);
+      template <class U>
+      explicit polynomial(const polynomial<U>& p);
       polynomial(std::initializer_list<T> l); // C++11
 
       // access:
@@ -46,18 +48,13 @@
       polynomial& operator /=(const U& value);
       template <class U>
       polynomial& operator %=(const U& value);
-      template <class U>
-      polynomial& operator >>=(const U& a);
-      template <class U>
-      polynomial& operator <<=(const U& a);
-      template <class U>
-      polynomial& operator +=(const polynomial<U>& value);
-      template <class U>
-      polynomial& operator -=(const polynomial<U>& value);
-      template <class U>
-      polynomial& operator *=(const polynomial<U>& value);
-      polynomial& operator /=(const polynomial<T>& value);
-      polynomial& operator %=(const polynomial<T>& value);
+      polynomial& operator >>=(int n);
+      polynomial& operator <<=(int n);
+      polynomial& operator +=(const polynomial& value);
+      polynomial& operator -=(const polynomial& value);
+      polynomial& operator *=(const polynomial& value);
+      polynomial& operator /=(const polynomial& value);
+      polynomial& operator %=(const polynomial& value);
    };
 
    template <class T>
@@ -93,9 +90,9 @@
    polynomial<T> operator - (const polynomial<T>& a); // Unary minus
 
    template <class T>
-   polynomial<T> operator >> (polynomial<T> const &a, const U& b);
+   polynomial<T> operator >> (polynomial<T> const &a, int b);
    template <class T>
-   polynomial<T> operator << (polynomial<T> const &a, const U& b);   
+   polynomial<T> operator << (polynomial<T> const &a, int b);   
    
    template <class T>
    bool operator == (const polynomial<T> &a, const polynomial<T> &b);

From 1115101f31353cfccd21763b4995c651421ed38b Mon Sep 17 00:00:00 2001
From: Kolya Matteo <kundor@kundor.org>
Date: Sun, 22 May 2016 21:52:52 -0400
Subject: [PATCH 15/15] Unnecessary include

---
 include/boost/math/tools/polynomial.hpp | 1 -
 1 file changed, 1 deletion(-)

diff --git a/include/boost/math/tools/polynomial.hpp b/include/boost/math/tools/polynomial.hpp
index 1c1fe8a5d..b82f1dae3 100644
--- a/include/boost/math/tools/polynomial.hpp
+++ b/include/boost/math/tools/polynomial.hpp
@@ -22,7 +22,6 @@
 #include <boost/math/policies/error_handling.hpp>
 #include <boost/math/special_functions/binomial.hpp>
 #include <boost/operators.hpp>
-#include <boost/type_traits/is_convertible.hpp>
 
 #include <vector>
 #include <ostream>