+template<typename T> static void
+canonicalize(T& p, const typename T::size_type hint = std::numeric_limits<typename T::size_type>::max())
+{
+ if (p.empty())
+ return;
+
+ std::size_t i = p.size() - 1;
+ // Be fast if the polynomial is already canonicalized
+ if (!zerop(p[i]))
+ return;
+
+ if (hint < p.size())
+ i = hint;
+
+ bool is_zero = false;
+ do {
+ if (!zerop(p[i])) {
+ ++i;
+ break;
+ }
+ if (i == 0) {
+ is_zero = true;
+ break;
+ }
+ --i;
+ } while (true);
+
+ if (is_zero) {
+ p.clear();
+ return;
+ }
+
+ p.erase(p.begin() + i, p.end());
+}
+
+// END COPY FROM UPOLY.HPP
+
+static void expt_pos(umodpoly& a, unsigned int q)
+{
+ if ( a.empty() ) return;
+ cl_MI zero = a[0].ring()->zero();
+ int deg = degree(a);
+ a.resize(degree(a)*q+1, zero);
+ for ( int i=deg; i>0; --i ) {
+ a[i*q] = a[i];
+ a[i] = zero;
+ }
+}
+
+template<typename T>
+static T operator+(const T& a, const T& b)
+{
+ int sa = a.size();
+ int sb = b.size();
+ if ( sa >= sb ) {
+ T r(sa);
+ int i = 0;
+ for ( ; i<sb; ++i ) {
+ r[i] = a[i] + b[i];
+ }
+ for ( ; i<sa; ++i ) {
+ r[i] = a[i];
+ }
+ canonicalize(r);
+ return r;
+ }
+ else {
+ T r(sb);
+ int i = 0;
+ for ( ; i<sa; ++i ) {
+ r[i] = a[i] + b[i];
+ }
+ for ( ; i<sb; ++i ) {
+ r[i] = b[i];
+ }
+ canonicalize(r);
+ return r;
+ }
+}
+
+template<typename T>
+static T operator-(const T& a, const T& b)
+{
+ int sa = a.size();
+ int sb = b.size();
+ if ( sa >= sb ) {
+ T r(sa);
+ int i = 0;
+ for ( ; i<sb; ++i ) {
+ r[i] = a[i] - b[i];
+ }
+ for ( ; i<sa; ++i ) {
+ r[i] = a[i];
+ }
+ canonicalize(r);
+ return r;
+ }
+ else {
+ T r(sb);
+ int i = 0;
+ for ( ; i<sa; ++i ) {
+ r[i] = a[i] - b[i];
+ }
+ for ( ; i<sb; ++i ) {
+ r[i] = -b[i];
+ }
+ canonicalize(r);
+ return r;
+ }
+}
+
+static upoly operator*(const upoly& a, const upoly& b)
+{
+ upoly c;
+ if ( a.empty() || b.empty() ) return c;
+
+ int n = degree(a) + degree(b);
+ c.resize(n+1, 0);
+ for ( int i=0 ; i<=n; ++i ) {
+ for ( int j=0 ; j<=i; ++j ) {
+ if ( j > degree(a) || (i-j) > degree(b) ) continue;
+ c[i] = c[i] + a[j] * b[i-j];
+ }
+ }
+ canonicalize(c);
+ return c;
+}
+
+static umodpoly operator*(const umodpoly& a, const umodpoly& b)
+{
+ umodpoly c;
+ if ( a.empty() || b.empty() ) return c;
+
+ int n = degree(a) + degree(b);
+ c.resize(n+1, a[0].ring()->zero());
+ for ( int i=0 ; i<=n; ++i ) {
+ for ( int j=0 ; j<=i; ++j ) {
+ if ( j > degree(a) || (i-j) > degree(b) ) continue;
+ c[i] = c[i] + a[j] * b[i-j];
+ }
+ }
+ canonicalize(c);
+ return c;
+}
+
+static upoly operator*(const upoly& a, const cl_I& x)
+{
+ if ( zerop(x) ) {
+ upoly r;
+ return r;
+ }
+ upoly r(a.size());
+ for ( size_t i=0; i<a.size(); ++i ) {
+ r[i] = a[i] * x;
+ }
+ return r;
+}
+
+static upoly operator/(const upoly& a, const cl_I& x)
+{
+ if ( zerop(x) ) {
+ upoly r;
+ return r;
+ }
+ upoly r(a.size());
+ for ( size_t i=0; i<a.size(); ++i ) {
+ r[i] = exquo(a[i],x);
+ }
+ return r;
+}
+
+static umodpoly operator*(const umodpoly& a, const cl_MI& x)
+{
+ umodpoly r(a.size());
+ for ( size_t i=0; i<a.size(); ++i ) {
+ r[i] = a[i] * x;
+ }
+ canonicalize(r);
+ return r;
+}
+
+static void upoly_from_ex(upoly& up, const ex& e, const ex& x)