From 47b7f9c9b4a5bd91c7c83b12bac8481b55bf4b92 Mon Sep 17 00:00:00 2001
From: Alexei Sheplyakov <varg@theor.jinr.ru>
Date: Tue, 14 Oct 2008 10:44:07 +0400
Subject: [PATCH 1/1] polynomial: introduce a helper function to divide
 polynomial by ring element.

---
 ginac/polynomial/ring_traits.hpp | 13 +++++++++++++
 ginac/polynomial/upoly.hpp       | 31 +++++++++++++++++++++++++++++++
 2 files changed, 44 insertions(+)

diff --git a/ginac/polynomial/ring_traits.hpp b/ginac/polynomial/ring_traits.hpp
index 10fb172a..540a6bd7 100644
--- a/ginac/polynomial/ring_traits.hpp
+++ b/ginac/polynomial/ring_traits.hpp
@@ -10,6 +10,19 @@ static inline cln::cl_I div(const cln::cl_I& x, const cln::cl_I& y)
 	return cln::exquo(x, y);
 }
 
+/// Exact integer division.
+/// Check if y divides x, if yes put the quotient into q, otherwise don't
+/// touch q. Returns true if y divides x and false if not.
+static inline bool div(cln::cl_I& q, const cln::cl_I& x, const cln::cl_I& y)
+{
+	const cln::cl_I_div_t qr = cln::truncate2(x, y);
+	if (zerop(qr.remainder)) {
+		q = qr.quotient;
+		return true;
+	}
+	return false;
+}
+
 static inline cln::cl_I get_ring_elt(const cln::cl_I& sample, const int val)
 {
 	return cln::cl_I(val);
diff --git a/ginac/polynomial/upoly.hpp b/ginac/polynomial/upoly.hpp
index 3564ae02..2e9092f1 100644
--- a/ginac/polynomial/upoly.hpp
+++ b/ginac/polynomial/upoly.hpp
@@ -8,6 +8,7 @@
 #include <stdexcept>
 #include <iterator>
 #include <limits>
+#include "ring_traits.hpp"
 #include "debug.hpp"
 #include "compiler.h"
 
@@ -113,6 +114,36 @@ operator*=(T& p, const typename T::value_type& c)
 	return p;
 }
 
+/// Divide the polynomial @a p by the ring element @a c, put the result
+/// into @a r. If @a p is not divisible by @a c @a r is undefined.
+template<typename T> bool divide(T& r, const T& p, const typename T::value_type& c)
+{
+	if (p.empty()) {
+		r.clear();
+		return true;
+	}
+	if (r.size() != p.size())
+		r.resize(p.size());
+	bool divisible = true;
+	for (std::size_t i = p.size(); i-- != 0; ) {
+		divisible = div(r[i], p[i], c);
+		if (!divisible)
+			break;
+	}
+	return divisible;
+}
+
+template<typename T> bool divide(T& p, const typename T::value_type& c)
+{
+	if (p.empty())
+		return true;
+	T r(p.size());
+	bool divisible = divide(r, p, c);
+	if (divisible)
+		swap(p, r);
+	return divisible;
+}
+
 // Convert Z[x] -> Z/p[x]
 
 static void
-- 
2.44.0