From 6080afee77ffcf81702bfc759437d80f36668692 Mon Sep 17 00:00:00 2001
From: "Jeremy W. Murphy" <jeremy.william.murphy@gmail.com>
Date: Sun, 1 Nov 2015 03:28:59 +1100
Subject: [PATCH] Fix quotient_remainder to use the correct ordering of
 coefficients.

Add equality operator and zero_element() for multiplication.
---
 include/boost/math/tools/polynomial.hpp | 73 ++++++++++++++++---------
 1 file changed, 46 insertions(+), 27 deletions(-)

diff --git a/include/boost/math/tools/polynomial.hpp b/include/boost/math/tools/polynomial.hpp
index afe3f0a8a..25c89fb7a 100644
--- a/include/boost/math/tools/polynomial.hpp
+++ b/include/boost/math/tools/polynomial.hpp
@@ -11,6 +11,7 @@
 #endif
 
 #include <boost/assert.hpp>
+#include <boost/config.hpp>
 #include <boost/math/tools/rational.hpp>
 #include <boost/math/tools/real_cast.hpp>
 #include <boost/math/special_functions/binomial.hpp>
@@ -109,37 +110,51 @@ template <typename T>
 std::pair< polynomial<T>, polynomial<T> >
 unchecked_synthetic_division(const polynomial<T>& dividend, const polynomial<T>& divisor)
 {
-    BOOST_ASSERT(dividend.degree() >= divisor.degree());
-    BOOST_ASSERT(divisor.degree() != 0 || divisor[0] != T(0));
+    BOOST_ASSERT(divisor.degree() <= dividend.degree());
     
     std::vector<T> intermediate_result(dividend.data());
+    if (divisor.degree() == T(0))
+        intermediate_result.insert(begin(intermediate_result), T(0));
     
     {
-        T const normalizer = divisor[0];
-        for (std::size_t i = 0; i < dividend.degree() - divisor.degree() + 1; i++)
+        typedef BOOST_DEDUCED_TYPENAME std::vector<T>::reverse_iterator reverse_iterator;
+        T const normalizer = divisor[divisor.size() - 1];
+        for (reverse_iterator i = intermediate_result.rbegin(); 
+             i != intermediate_result.rbegin() + dividend.degree() - divisor.degree() + 1; 
+             i++)
         {
-            if (intermediate_result[i] != T(0))
+            if (*i != T(0))
             {
-                T const coefficient = intermediate_result[i] /= normalizer;
+                T const coefficient = *i /= normalizer;
                 for (std::size_t j = 1; j < divisor.degree() + 1; j++)
-                {
-                    intermediate_result[i + j] += -divisor[j] * coefficient;
-                }
+                    *(i + j) += -*(divisor.data().rbegin() + j) * coefficient;
             }
         }
     }
     
     {
         typedef BOOST_DEDUCED_TYPENAME std::vector<T>::iterator iterator;
-        iterator const f = intermediate_result.begin();
-        iterator const m = f + dividend.degree() - divisor.degree() + 1;
-        iterator const l = m + (*m == T(0) ? 1 : divisor.degree());
+        iterator const f = intermediate_result.begin(); // remainder
+        iterator const m = f + std::max(divisor.degree(), 1ul); // quotient
+        iterator const l = intermediate_result.end();
         BOOST_ASSERT(m - f > 0);
         BOOST_ASSERT(l - m > 0);
-        return std::make_pair(polynomial<T>(f, m), polynomial<T>(m, l));
+        polynomial<T> const quotient(m, l);
+        polynomial<T> const remainder = *f != T(0) ? polynomial<T>(f, m) : polynomial<T>();
+        return std::make_pair(quotient, remainder);
     }
 }
 
+
+/**
+ * Returns the zero element for multiplication of polynomials.
+ */
+template <class T>
+polynomial<T> zero_element(std::multiplies< polynomial<T> >)
+{
+    return polynomial<T>();
+}
+
 /* Calculates a / b and a % b, returning the pair (quotient, remainder) together
  * because the same amount of computation yields both.
  */
@@ -147,22 +162,20 @@ template <typename T>
 std::pair< polynomial<T>, polynomial<T> >
 quotient_remainder(const polynomial<T>& dividend, const polynomial<T>& divisor)
 {
-    if (divisor.degree() == 0 && divisor[0] == T(0))
+    polynomial<T> const zero = zero_element(std::multiplies< polynomial<T> >());
+    if (divisor == zero)
         throw std::domain_error("Divide by zero.");
-    
     if (dividend.degree() < divisor.degree())
-    {
-        return std::make_pair(polynomial<T>(T(0)), dividend);
-    }
-    
+        return std::make_pair(zero, dividend); // TODO: Is this right?
     return unchecked_synthetic_division(dividend, divisor);
 }
 
 
 template <class T>
 class polynomial : 
-    dividable< polynomial<int>
-    , modable< polynomial<int> > >
+    equality_comparable< polynomial<T>, 
+    dividable< polynomial<T>, 
+    modable< polynomial<T> > > >
 {
 public:
    // typedefs:
@@ -178,6 +191,12 @@ public:
    {
    }
    
+   template <class I>
+   polynomial(I first, I last)
+   : m_data(first, last)
+   {
+   }
+   
    template <class U>
    polynomial(const U& point)
    {
@@ -197,12 +216,6 @@ public:
       }
    }
 
-   template <typename I>
-   polynomial(I first, I last)
-      : m_data(first, last)
-   {
-   }
-   
    // access:
    size_type size()const { return m_data.size(); }
    size_type degree()const { return m_data.size() - 1; }
@@ -390,6 +403,12 @@ inline polynomial<T> operator * (const U& a, const polynomial<T>& b)
    return result;
 }
 
+template <class T>
+bool operator == (const polynomial<T> &x, const polynomial<T> &y)
+{
+    return x.data() == y.data();
+}
+
 template <class charT, class traits, class T>
 inline std::basic_ostream<charT, traits>& operator << (std::basic_ostream<charT, traits>& os, const polynomial<T>& poly)
 {