From 240e4e6e0fa73353d4fbc2e44fa91084a6dce790 Mon Sep 17 00:00:00 2001
From: "Jeremy W. Murphy" <jeremy.william.murphy@gmail.com>
Date: Sun, 25 Oct 2015 12:05:32 +1100
Subject: [PATCH] Separate out into unchecked_synthetic_division().

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

diff --git a/include/boost/math/tools/polynomial.hpp b/include/boost/math/tools/polynomial.hpp
index 4f1d277d2..27d59f2dc 100644
--- a/include/boost/math/tools/polynomial.hpp
+++ b/include/boost/math/tools/polynomial.hpp
@@ -103,15 +103,14 @@ template <typename T>
 class polynomial;
 
 
-/* Calculates a / b and a % b, returning the pair (quotient, remainder) together
- * because the same amount of computation yields both.
- */
+
 template <typename T>
 std::pair< polynomial<T>, polynomial<T> >
-quotient_remainder(const polynomial<T>& dividend, const polynomial<T>& divisor)
+unchecked_synthetic_division(const polynomial<T>& dividend, const polynomial<T>& divisor)
 {
-    if (divisor.degree() == 0 && divisor[0] == T(0))
-        throw std::domain_error("Divide by zero.");
+    BOOST_ASSERT(dividend.degree() >= divisor.degree());
+    BOOST_ASSERT(divisor.degree() != 0 || divisor[0] != T(0));
+    
     std::vector<T> intermediate_result(dividend.data());
     
     {
@@ -131,13 +130,31 @@ quotient_remainder(const polynomial<T>& dividend, const polynomial<T>& divisor)
     
     {
         typedef BOOST_DEDUCED_TYPENAME std::vector<T>::iterator iterator;
-        iterator f = intermediate_result.begin();
-        iterator m = f + dividend.degree() - divisor.degree() + 1;
-        iterator l = m + divisor.degree();
+        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());
         return std::make_pair(polynomial<T>(f, m), polynomial<T>(m, l));
     }
 }
 
+/* Calculates a / b and a % b, returning the pair (quotient, remainder) together
+ * because the same amount of computation yields both.
+ */
+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))
+        throw std::domain_error("Divide by zero.");
+    
+    if (dividend.degree() < divisor.degree())
+    {
+        return std::make_pair(polynomial<T>(T(0)), dividend);
+    }
+    
+    return unchecked_synthetic_division(dividend, divisor);
+}
+
 
 template <class T>
 class polynomial