diff --git a/doc/internals/polynomial.qbk b/doc/internals/polynomial.qbk
index 3888d5420..1e668fab8 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:
@@ -50,16 +52,13 @@
       polynomial& operator /=(const U& value);
       template <class U>
       polynomial& operator %=(const 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);
-      template <class U>
-      polynomial& operator /=(const polynomial<U>& value);
-      template <class U>
-      polynomial& operator %=(const polynomial<U>& 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>
@@ -92,16 +91,12 @@
    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);
+   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);
@@ -109,12 +104,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);
diff --git a/include/boost/math/tools/polynomial.hpp b/include/boost/math/tools/polynomial.hpp
index df1bde374..d373f17d7 100644
--- a/include/boost/math/tools/polynomial.hpp
+++ b/include/boost/math/tools/polynomial.hpp
@@ -255,12 +255,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:
@@ -296,7 +295,7 @@ class polynomial :
       : 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)
       {
@@ -388,30 +387,37 @@ class polynomial :
    }
 
    template <class U>
-   polynomial& operator %=(const U& /*value*/)
-   {
-       // We can always divide by a scalar, so there is no remainder:
-       this->set_zero();
+   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.
+       if (std::numeric_limits<T>::is_integer)
+       {
+           modulus(value);
+           normalize();
+       }
+       else
+       {
+           set_zero();
+       }
        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!!!
       if (!value)
@@ -427,31 +433,28 @@ class polynomial :
       return *this;
    }
 
-   template <typename U>
-   polynomial& operator /=(const polynomial<U>& value)
+   polynomial& operator /=(const polynomial& value)
    {
        *this = quotient_remainder(*this, value).first;
        return *this;
    }
 
-   template <typename U>
-   polynomial& operator %=(const polynomial<U>& value)
+   polynomial& operator %=(const polynomial& value)
    {
        *this = quotient_remainder(*this, value).second;
        return *this;
    }
 
-   template <typename U>
-   polynomial& operator >>=(U const &n)
+   polynomial& operator >>=(int n)
    {
-       BOOST_ASSERT(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;
    }
 
-   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;
@@ -554,64 +557,59 @@ class polynomial :
         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;
 };
 
 
-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>
-inline polynomial<T> operator - (const polynomial<T>& a, const polynomial<T>& b)
+template <class T, class U>
+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>
-inline polynomial<T> operator * (const polynomial<T>& a, const polynomial<T>& b)
+template <class T, class U>
+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>
@@ -623,51 +621,57 @@ 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>
-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 <typename T, typename U>
-polynomial<T> operator >> (const polynomial<T>& a, const U& b)
+template <class T>
+inline bool operator < (const polynomial<T> &a, const polynomial<T> &b)
 {
-    polynomial<T> result(a);
-    result >>= b;
-    return result;
+    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)
+template <typename T>
+inline polynomial<T> operator >> (polynomial<T> a, int b)
 {
-    polynomial<T> result(a);
-    result <<= b;
-    return result;
+    a >>= b;
+    return a;
+}
+
+template <typename T>
+inline polynomial<T> operator << (polynomial<T> a, int b)
+{
+    a <<= b;
+    return a;
 }
 
 // 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);
 }
diff --git a/test/test_polynomial.cpp b/test/test_polynomial.cpp
index 0ce45354a..9c37fb815 100644
--- a/test/test_polynomial.cpp
+++ b/test/test_polynomial.cpp
@@ -323,6 +323,63 @@ BOOST_AUTO_TEST_CASE_TEMPLATE( test_pow, T, all_test_types )
     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]));
+        }
+    }
+}
 
 BOOST_AUTO_TEST_CASE_TEMPLATE(test_bool, T, all_test_types)
 {