/Users/mourrain/Devel/mmx/algebramix/include/algebramix/polynomial.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : polynomial.hpp
* DESCRIPTION: Dense univariate polynomials
* COPYRIGHT  : (C) 2004  Joris van der Hoeven and Gregoire Lecerf
*******************************************************************************
* This software falls under the GNU general public license and comes WITHOUT
* ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
* If you don't have this file, write to the Free Software Foundation, Inc.,
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
******************************************************************************/

#ifndef __MMX_POLYNOMIAL_HPP
#define __MMX_POLYNOMIAL_HPP
#include <algebramix/polynomial_naive.hpp>

namespace mmx {
#define Polynomial_variant(C) polynomial_variant_helper<C>::PV
#define Polynomial_multiply(C) polynomial_variant_helper<C>::PV
#define TMPL_DEF \
  template<typename C, typename V=typename Polynomial_variant(C) >
#define TMPL template<typename C, typename V>
#define TMPLK template<typename C, typename V, typename K>
#define Format format<C>
#define Vector vector<C>
#define Polynomial polynomial<C,V>
#define Polynomial_rep polynomial_rep<C,V>
#define Abs_polynomial polynomial<Abs_type(C),V>
#define Real_polynomial polynomial<Real_type(C),V>
#define Center_polynomial polynomial<Center_type(C),V>
#define Radius_polynomial polynomial<Radius_type(C),V>
TMPL class polynomial_rep;
TMPL class polynomial;
TMPL inline int deg (const Polynomial& P);
TMPL inline int degree (const Polynomial& P);
TMPL inline nat N (const Polynomial& P);
TMPL inline C* seg (Polynomial& P);
TMPL inline const C* seg (const Polynomial& P);

/******************************************************************************
* Dense polynomials
******************************************************************************/

TMPL_DEF
class polynomial_rep REP_STRUCT_1(C) {
public:
  typedef implementation<vector_linear,V> Vec;
  typedef implementation<polynomial_linear,V> Pol;
private:
  C*  a;
  nat n;
  nat l;  // allocated number of coefficients
public:
  inline void normalize () {
    while ((n != 0) && is_exact_zero (a[n-1])) n--; }
  inline polynomial_rep (C* a2, nat n2, const Format& fm):
    Format (fm), a(a2), n(n2), l(n2) { normalize (); }
  inline polynomial_rep (C* a2, nat n2, nat l2, const Format& fm):
    Format (fm), a(a2), n(n2), l(l2) { normalize (); }
  inline ~polynomial_rep () {
    mmx_delete<C> (a, l); }
  void extend (nat n2) {
    nat l2= aligned_size<C,V> (n2);
    C* b= mmx_formatted_new<C> (l2, this->tfm ());
    Pol::copy (b, a, n);
    Pol::clear (b+n, l2-n);
    mmx_delete<C> (a, l);
    a= b;
    n= n2;
    l= l2;
  }
  void resize (nat n2) {
    if (n < n2) extend (n2);
    n= n2;
  }

  friend class Polynomial;
  friend int deg LESSGTR (const Polynomial& P);
  friend int degree LESSGTR (const Polynomial& P);
  friend nat N LESSGTR (const Polynomial& P);
  friend C* seg LESSGTR (Polynomial& P);
  friend const C* seg LESSGTR (const Polynomial& P);
};
  
TMPL_DEF
class polynomial {
INDIRECT_PROTO_2 (polynomial, polynomial_rep, C, V)
public:
  typedef implementation<vector_linear,V> Vec;
  typedef implementation<polynomial_linear,V> Pol;
  typedef typename Pol::template global_variables<Polynomial> S;
public:
  static inline generic get_variable_name () {
    return S::get_variable_name (); }
  static inline void set_variable_name (const generic& x) {
    S::set_variable_name (x); }

  //#ifdef EMPTY_FORMATS
  inline polynomial () {
    nat l= aligned_size<C,V> (0);
    rep= new Polynomial_rep (mmx_new<C> (l), 0, l, Format (no_format ()));
  }
  //#endif
  inline polynomial (const Format& fm) {
    nat l= aligned_size<C,V> (0);
    rep= new Polynomial_rep (mmx_new<C> (l), 0, l, fm);
  }
  template<typename T> inline polynomial (const T& x) {
    nat l= aligned_size<C,V> (1);
    C c= as<C> (x);
    C* a= mmx_formatted_new<C> (l, get_format (c));
    a[0]= c;
    rep= new Polynomial_rep (a, 1, l, get_format (c));
    rep->normalize ();
  }
  //#ifdef EMPTY_FORMATS
  template<typename T> inline
  polynomial (const polynomial<T,V>& P) {
    Format fm= as<Format > (CF(P));
    nat l= aligned_size<C,V> (N(P));
    C* a= mmx_formatted_new<C> (l, fm);
    Pol::cast (a, seg (P), N(P));
    rep= new Polynomial_rep (a, N(P), l, fm);
    rep->normalize ();
  }
  //#endif
  template<typename T> inline
  polynomial (const polynomial<T,V>& P, const Format& fm) {
    nat l= aligned_size<C,V> (N(P));
    C* a= mmx_formatted_new<C> (l, fm);
    Pol::cast (a, seg (P), N(P));
    rep= new Polynomial_rep (a, N(P), l, fm);
    rep->normalize ();
  }
  polynomial (const C& x, nat n) {
    nat l= aligned_size<C,V> (n+1);
    C* a= mmx_formatted_new<C> (l, get_format (x));
    Pol::clear (a, n);
    a[n]= x;
    rep= new Polynomial_rep (a, n+1, l, get_format (x));
    rep->normalize ();
  }
  inline polynomial (C* a, nat n, const Format& fm) {
    rep= new Polynomial_rep (a, n, fm);
  }
  inline polynomial (C* a, nat n, nat l, const Format& fm) {
    rep= new Polynomial_rep (a, n, l, fm);
  }
  polynomial (const vector<C>& v) {
    nat l= aligned_size<C,V> (N(v));
    C* a= mmx_formatted_new<C> (l, CF(v));
    Pol::copy (a, seg (v), N(v));
    rep= new Polynomial_rep (a, N(v), l, CF(v));
    rep->normalize ();
  }
  polynomial (const iterator<C>& it) {
    nat i, l=0;
    nat s= aligned_size<C,V> (0);
    C* a= mmx_formatted_new<C> (s, CF(it));
    rep= new Polynomial_rep (a, 0, s, CF(it));
    for (i=0; busy (it); i++, ++it) {
      if (i >= l) {
        l= max ((nat) 1, l<<1);
        rep->extend (l);
      }
      rep->a[i]= *it;
    }
    rep->normalize ();
  }
  inline friend Polynomial copy (const Polynomial& P) {
    nat n= N(P), l= aligned_size<C,V> (n);
    C* r= mmx_formatted_new<C> (l, CF(P));
    Pol::copy (r, seg (P), n);
    return Polynomial (r, n, l, CF(P));
  }
  inline C operator[] (nat i) const {
    return i<rep->n? rep->a[i]: promote (0, CF(*this)); }
  Polynomial& operator += (const Polynomial& P);
  Polynomial& operator -= (const Polynomial& P);
  Polynomial& operator *= (const C& c);
  Polynomial& operator /= (const C& c);
};
INDIRECT_IMPL_2 (polynomial, polynomial_rep, typename C, C, typename V, V)
DEFINE_UNARY_FORMAT_2 (polynomial)

TMPL inline nat N (const Polynomial& P) { return P->n; }
TMPL inline int deg (const Polynomial& P) { return ((int) P->n) - 1; }
TMPL inline int degree (const Polynomial& P) { return ((int) P->n) - 1; }
TMPL inline C* seg (Polynomial& P) { return P->a; }
TMPL inline const C* seg (const Polynomial& P) { return P->a; }
TMPL inline Format CF (const Polynomial& P) { return P->tfm (); }

TMPL inline vector<C>
coefficients (const Polynomial& P) {
  vector<C> v (get_sample (CF(P)), N(P));
  for (nat i= 0; i < N(P); i++) v[i]= P[i];
  return v;
}

TMPL inline vector<C>
coefficients (const Polynomial& P, nat b, nat e) {
  if (b >= e) return vector<C> ();
  vector<C> v (get_sample (CF(P)), e - b);
  for (nat i= 0; i < e - b; i++) v[i]= P[b+i];
  return v;
}

STYPE_TO_TYPE(TMPL,scalar_type,Polynomial,C);
STYPE_TO_TYPE(TMPL,monomial_type,Polynomial,nat);
UNARY_RETURN_TYPE(TMPL,Re_op,Polynomial,Real_polynomial);
UNARY_RETURN_TYPE(TMPL,abs_op,Polynomial,Abs_polynomial);
UNARY_RETURN_TYPE(TMPL,center_op,Polynomial,Center_polynomial);
UNARY_RETURN_TYPE(TMPL,radius_op,Polynomial,Radius_polynomial);

TMPL
struct fast_helper<Polynomial > {
  typedef Fast_type(C) FC;
  typedef polynomial<FC> fast_type;
  typedef typename Polynomial_variant(FC) FV;
  static inline fast_type dd (const Polynomial& p) {
    format<FC> fm= fast (CF(p));
    nat l= aligned_size<FC,FV> (N(p));
    FC* r= mmx_formatted_new<FC> (l, fm);
    for (nat i=0; i<N(p); i++) r[i]= fast (p[i]);
    return fast_type (r, N(p), l, fm); }
  static inline Polynomial uu (const fast_type& p) {
    Format fm= slow<C> (CF(p));
    nat l= aligned_size<C,V> (N(p));
    C* r= mmx_formatted_new<C> (l, fm);
    for (nat i=0; i<N(p); i++) r[i]= slow<C> (p[i]);
    return Polynomial (r, N(p), l, fm); }
};

template<typename T, typename F, typename TV, typename FV>
struct as_helper<polynomial<T,TV>,polynomial<F,FV> > {
  static inline polynomial<T,TV>
  cv (const polynomial<F,FV>& p) {
    format<T> fm= as<format<T> > (CF(p));
    nat l= aligned_size<T,TV> (N(p));
    T* c= mmx_formatted_new<T> (l, fm);
    for (nat i=0; i<N(p); i++) c[i]= as<T> (p[i]);
    //return polynomial<T,TV> (c, N(p), l, get_format (c[0]));
    // FIXME: a direct return does not seem to work; a compiler bug?
    polynomial<T,TV> r= polynomial<T,TV> (c, N(p), l, fm);
    return r;
  }
};

template<typename T, typename F, typename TV, typename FV> inline void
set_as (polynomial<T,TV>& r, const polynomial<F,FV>& p) {
  r.secure ();
  inside (r) -> resize (N(p));
  T* c= seg (r);
  for (nat i=0; i<N(p); i++) set_as (c[i], p[i]);
  inside (r) -> normalize ();
}

template<typename C, typename V, typename T> inline void
set_as (Polynomial& r, const T& x) {
  r.secure ();
  inside (r) -> resize (1);
  set_as (seg(r)[0], x);
  inside (r) -> normalize ();
}

template<typename Op, typename C, typename V> nat
unary_hash (const Polynomial& p) {
  register nat i, h= 642531, n= N(p);
  for (i=0; i<n; i++)
    h= (h<<1) ^ (h<<5) ^ (h>>27) ^ Op::op (p[i]);
  return h;
}

/******************************************************************************
* Polynomial iteration
******************************************************************************/

TMPL_DEF
class polynomial_iterator_rep: public iterator_rep<C> {
  Polynomial P;  // the polynomial to traverse
  nat i;         // current position
protected:
  bool is_busy () { return ((int) i) <= deg (P); }
  void advance () { i++; }
  C current () { return P[i]; }
public:
  polynomial_iterator_rep (const Polynomial& P2):
    iterator_rep<C> (CF(P2)), P(P2), i(0) {}
};

TMPL iterator<C>
iterate (const Polynomial& P) {
  return iterator<C> (new polynomial_iterator_rep<C,V> (P));
}

/******************************************************************************
* Output routines for dense polynomials
******************************************************************************/

TMPL inline generic var (const Polynomial& P) {
  (void) P; return Polynomial::get_variable_name (); }

TMPL inline void set_variable_name (const Polynomial& P, const generic& x) {
  (void) P; return Polynomial::set_variable_name (x); }

TMPL syntactic
flatten (const Polynomial& P, const syntactic& v) {
  int i, len= deg(P);
  syntactic s= 0;
  for (i=len; i>=0; i--)
    s= s + flatten (P[i]) * pow (v, i);
  return s;
}

TMPL syntactic
flatten (const Polynomial& P) {
  return flatten (P, as_syntactic (var (P)));
}

TMPL
struct binary_polynomial_helper: public void_binary_helper<Polynomial > {
  static inline string short_type_name () {
    return "P" * Short_type_name (C); }
  static inline generic full_type_name () {
    return gen ("Polynomial", Full_type_name (C)); }
};

TMPL
struct binary_helper<Polynomial >: public binary_polynomial_helper<C,V> {
  static inline generic disassemble (const Polynomial& p) {
    return binary_disassemble<vector<C> > (coefficients (p)); }
  static inline Polynomial assemble (const generic& v) {
    return Polynomial (binary_assemble<vector<C> > (v)); }
  static inline void write (const port& out, const Polynomial& p) {
    binary_write<Format > (out, CF(p));
    binary_write<nat> (out, N(p));
    for (nat i=0; i<N(p); i++)
      binary_write<C> (out, p[i]); }
  static inline Polynomial read (const port& in) {
    Format fm= binary_read<Format > (in);
    nat n= binary_read<nat> (in);
    nat l= aligned_size<C,V> (n);
    C* r= mmx_formatted_new<C> (l, fm);
    for (nat i=0; i<n; i++)
      r[i]= binary_read<C> (in);
    return Polynomial (r, n, l, fm); }
};

/******************************************************************************
* Valuation and comparisons
******************************************************************************/

TMPL int
val (const Polynomial& P) {
  for (nat i=0; i<N(P); i++)
    if (P[i] != 0) return (int) i;
  return (int) (((nat) (-1)) >> 1);
}

TMPL int
sign (const Polynomial& P) {
  for (nat i=0; i<N(P); i++) {
    int s= sign (P[i]);
    if (s != 0) return s;
  }
  return 0;
}

TMPL int
compare (const Polynomial& P1, const Polynomial& P2) {
  nat n1= N(P1), n2= N(P2), n= max (n1, n2);
  for (nat i=0; i<n; i++) {
    int s= sign (P1[i] - P2[i]);
    if (s != 0) return s;
  }
  return 0;
}

template<typename Op,typename C,typename V> inline bool
binary_test (const Polynomial& P1, const Polynomial& P2) {
  nat n1= N(P1), n2= N(P2), n= max (n1, n2);
  for (nat i=0; i<n; i++)
    if (Op::not_op (P1[i], P2[i])) return false;
  return true;
}

TRUE_IDENTITY_OP_SUGAR(TMPL,Polynomial)
EXACT_IDENTITY_OP_SUGAR(TMPL,Polynomial)
COMPARE_SUGAR(TMPL,Polynomial)
EQUAL_SCALAR_SUGAR(TMPL,Polynomial,C)
COMPARE_SCALAR_SUGAR(TMPL,Polynomial,C)
EQUAL_INT_SUGAR(TMPL,Polynomial)
COMPARE_INT_SUGAR(TMPL,Polynomial)

/******************************************************************************
* In place arithmetic
******************************************************************************/

TMPL Polynomial&
Polynomial::operator += (const Polynomial& P) {
  typedef implementation<polynomial_linear,V> Pol;
  secure ();
  nat n1= rep->n, n2= N(P);
  if (n1 < n2) rep -> extend (n2);
  Pol::add (rep->a, seg (P), n2);
  if (n1 == n2) rep -> normalize ();
  return *this;
}

TMPL Polynomial&
Polynomial::operator -= (const Polynomial& P) {
  typedef implementation<polynomial_linear,V> Pol;
  secure ();
  nat n1= rep->n, n2= N(P);
  if (n1 < n2) rep->extend (n2);
  Pol::sub (rep->a, seg (P), n2);
  if (n1 == n2) rep->normalize ();
  return *this;
}

TMPL Polynomial&
Polynomial::operator *= (const C& c) {
  typedef implementation<polynomial_linear,V> Pol;
  secure ();
  Pol::mul_sc (rep->a, c, rep->n);
  rep->normalize ();
  return *this;
}

TMPL Polynomial&
Polynomial::operator /= (const C& c) {
  typedef implementation<polynomial_linear,V> Pol;
  secure ();
  Pol::div_sc (rep->a, c, rep->n);
  rep->normalize ();
  return *this;
}

/******************************************************************************
* Basic arithmetic operations
******************************************************************************/

TMPL Polynomial
operator - (const Polynomial& P) {
  typedef implementation<polynomial_linear,V> Pol;
  nat n= N(P);
  nat l= aligned_size<C,V> (n);
  C* r= mmx_formatted_new<C> (l, CF(P));
  Pol::neg (r, seg (P), n);
  return Polynomial (r, n, l, CF(P));
}

TMPL Polynomial
operator + (const Polynomial& P1, const Polynomial& P2) {
  typedef implementation<polynomial_linear,V> Pol;
  nat m= min (N(P1), N(P2));
  nat n= max (N(P1), N(P2));
  nat l= aligned_size<C,V> (n);
  C* r= mmx_formatted_new<C> (l, CF(P1));
  Pol::add (r, seg (P1), seg (P2), m);
  if (N(P1) > m) Pol::copy (r+m, seg(P1)+m, n-m);
  if (N(P2) > m) Pol::copy (r+m, seg(P2)+m, n-m);
  return Polynomial (r, n, l, CF(P1));
}

TMPL Polynomial operator + (const Polynomial& P, const C& c) {
  return P + Polynomial (c); }
TMPL Polynomial operator + (const C& c, const Polynomial& P) {
  return Polynomial (c) + P; }

TMPL Polynomial
big_add (const vector<Polynomial >& a) {
  typedef implementation<polynomial_linear,V> Pol;
  nat i, k= N(a), n=0;
  nat l= aligned_size<C,V> (n);
  for (i=0; i<k; i++) n= max (N(a[i]), n);
  C* r= mmx_formatted_new<C> (l, CF (get_sample (CF (a))));
  Pol::set (r, 0, n);
  for (i=0; i<k; i++) Pol::add (r, seg (a[i]), N(a[i]));
  return Polynomial (r, n, l, Format (CF(a)));
}

TMPL Polynomial
operator - (const Polynomial& P1, const Polynomial& P2) {
  typedef implementation<polynomial_linear,V> Pol;
  nat m= min (N(P1), N(P2));
  nat n= max (N(P1), N(P2));
  nat l= aligned_size<C,V> (n);
  C* r= mmx_formatted_new<C> (l, CF(P1));
  Pol::sub (r, seg (P1), seg (P2), m);
  if (N(P1) > m) Pol::copy (r+m, seg(P1)+m, n-m);
  if (N(P2) > m) Pol::neg  (r+m, seg(P2)+m, n-m);
  return Polynomial (r, n, l, CF(P1));
}

TMPL Polynomial operator - (const Polynomial& P, const C& c) {
  return P - Polynomial (c); }
TMPL Polynomial operator - (const C& c, const Polynomial& P) {
  return Polynomial (c) - P; }

TMPL Polynomial
operator * (const Polynomial& P1, const Polynomial& P2) {
  typedef implementation<polynomial_multiply,V> Pol;
  nat n1= N(P1), n2= N(P2);
  if (n1 == 0 || n2 == 0) return Polynomial (CF(P1));
  nat l= aligned_size<C,V> (n1+n2-1);
  C* r= mmx_formatted_new<C> (l, CF(P1));
  Pol::mul (r, seg (P1), seg (P2), n1, n2);
  return Polynomial (r, n1+n2-1, l, CF(P1));
}

TMPL Polynomial
tmul (int d2, const Polynomial& P1, const Polynomial& P2) {
  typedef implementation<polynomial_multiply,V> Pol;
  nat n2 = max (0, d2 + 1), n1= N(P1), n= N(P2);
  if (n1 == 0 || n2 == 0) return Polynomial (CF(P1));
  ASSERT (n < n1 + n2, "bad dimension in tmul");
  nat l2= aligned_size<C,V> (n2);
  C* r= mmx_formatted_new<C> (l2, CF(P1));
  if (n != n1 + n2 - 1) {
    nat l= aligned_size<C,V> (n1+n2-1);
    C* s2= mmx_formatted_new<C> (l, CF(P1));
    Pol::copy (s2, seg (P2), n);
    Pol::clear (s2 + n, n1+n2-n-1);
    Pol::tmul (r, seg (P1), s2, n1, n2);
    mmx_delete<C> (s2, l);
  }
  else
    Pol::tmul (r, seg (P1), seg (P2), n1, n2);
  return Polynomial (r, n2, l2, CF(P1));
}

TMPL Polynomial
square (const Polynomial& P) {
  typedef implementation<polynomial_multiply,V> Pol;
  nat n= N(P);
  if (n == 0) return P;
  nat l= aligned_size<C,V> (2*n-1);
  C* r= mmx_formatted_new<C> (l, CF(P));
  Pol::square (r, seg (P), n);
  return Polynomial (r, 2*n-1, l, CF(P));
}

TMPL Polynomial
operator * (const Polynomial& P, const C& c) {
  typedef implementation<polynomial_linear,V> Pol;
  nat n= N(P);
  nat l= aligned_size<C,V> (n);
  C* r= mmx_formatted_new<C> (l, CF(P));
  Pol::mul_sc (r, seg (P), c, n);
  return Polynomial (r, n, l, CF(P));
}

TMPL Polynomial
operator * (const C& c, const Polynomial& P) {
  typedef implementation<polynomial_linear,V> Pol;
  nat n= N(P);
  nat l= aligned_size<C,V> (n);
  C* r= mmx_formatted_new<C> (l, CF(P));
  Pol::mul_sc (r, seg (P), c, n);
  return Polynomial (r, n, l, CF(P));
}

TMPL Polynomial
big_mul (const vector<Polynomial >& a) {
  typedef implementation<vector_linear,V> Vec;
  return Vec::template vec_unary_big_dicho<mul_op> (seg (a), N(a));
}

TMPL Polynomial
operator / (const Polynomial& P, const C& c) {
  typedef implementation<polynomial_linear,V> Pol;
  nat n= N(P);
  nat l= aligned_size<C,V> (n);
  C* r= mmx_formatted_new<C> (l, CF(P));
  Pol::div_sc (r, seg (P), c, n);
  return Polynomial (r, n, l, CF(P));
}

TMPL Polynomial
operator / (const Polynomial& P1, const Polynomial& P2) {
  typedef implementation<polynomial_exact_divide,V> Pol;
  nat n1= N(P1), n2= N(P2);
  ASSERT (n2 != 0, "division by zero");
  if (n1 < n2) return Polynomial (CF(P1));
  nat lq= aligned_size<C,V> (n1-n2+1);
  C* q= mmx_formatted_new<C> (lq, CF(P1));
  Pol::div (q, seg (P1), seg (P2), n1, n2);
  return Polynomial (q, n1-n2+1, lq, CF(P1));
}

TMPL Polynomial
skew_div (const Polynomial& P, const C& c, bool left) {
  typedef implementation<polynomial_linear,V> Pol;
  if (left) {
    nat n= N(P);
    nat l= aligned_size<C,V> (n);
    C* r= mmx_formatted_new<C> (n, CF(P));
    Pol::div_sc (r, seg (P), c, n);
    return Polynomial (r, n, l, CF(P));
  }
  else return P / c;
}

TMPL Polynomial
quo (const Polynomial& P, const C& c) {
  typedef implementation<polynomial_linear,V> Pol;
  nat n= N(P);
  nat l= aligned_size<C,V> (n);
  C* r= mmx_formatted_new<C> (l, CF(P));
  Pol::quo_sc (r, seg (P), c, n);
  return Polynomial (r, n, l, CF(P));
}

TMPL Polynomial
rem (const Polynomial& P, const C& c) {
  typedef implementation<polynomial_linear,V> Pol;
  nat n= N(P);
  nat l= aligned_size<C,V> (n);
  C* r= mmx_formatted_new<C> (l, CF(P));
  Pol::rem_sc (r, seg (P), c, n);
  return Polynomial (r, n, l, CF(P));
}

TMPL Polynomial
quo (const Polynomial& P1, const Polynomial& P2) {
  typedef implementation<polynomial_divide,V> Pol;
  nat n1= N(P1), n2= N(P2);
  if (n1 < n2 || n2 == 0) return Polynomial (CF(P1));
  nat lq= aligned_size<C,V> (n1-n2+1);
  nat lr= aligned_size<C,V> (n1);
  C* q= mmx_formatted_new<C> (lq, CF(P1));
  C* r= mmx_formatted_new<C> (lr, CF(P1));
  Pol::copy (r, seg (P1), n1);
  Pol::quo_rem (q, r, seg (P2), n1, n2);
  mmx_delete<C> (r, lr);
  return Polynomial (q, n1-n2+1, lq, CF(P1));
}

TMPL Polynomial
tquo (int d1, const Polynomial& P1, const Polynomial& P2) {
  typedef implementation<polynomial_divide,V> Pol;
  nat n1= max (0, d1+1), n= N(P1), n2= N(P2);
  ASSERT (n <= n1-n2+1, "bad dimension in tquo");
  if (n1 < n2 || n2 == 0) return Polynomial (CF(P1));
  nat l= aligned_size<C,V> (n1);
  C* q= mmx_formatted_new<C> (l, CF(P1));
  C* r= mmx_formatted_new<C> (l, CF(P1));
  Pol::clear (r, n1);
  Pol::copy  (r+n2-1, seg (P1), n);
  Pol::tquo_rem (q, r, seg (P2), n1, n2);
  mmx_delete<C> (r, l);
  return Polynomial (q, n1, l, CF(P1));
}

TMPL Polynomial
rem (const Polynomial& P1, const Polynomial& P2) {
  typedef implementation<polynomial_divide,V> Pol;
  nat n1= N(P1), n2= N(P2);
  if (n1 < n2 || n2 == 0) return P1;
  nat lq= aligned_size<C,V> (n1-n2+1);
  nat lr= aligned_size<C,V> (n1);
  C* q= mmx_formatted_new<C> (lq, CF(P1));
  C* r= mmx_formatted_new<C> (lr, CF(P1));
  Pol::copy (r, seg (P1), n1);
  Pol::quo_rem (q, r, seg (P2), n1, n2);
  mmx_delete<C> (q, lq);
  return Polynomial (r, n2 - 1, lr, CF(P1));
}

TMPL bool
divides (const Polynomial& P1, const Polynomial& P2) {
  return rem (P2, P1) == 0;
}

TMPL Polynomial
rem (const Polynomial& P1, const Polynomial& P2, Polynomial& Q) {
  typedef implementation<polynomial_divide,V> Pol;
  nat n1= N(P1), n2= N(P2);
  if (n1 < n2 || n2 == 0) return Polynomial (CF(P1));
  nat lq= aligned_size<C,V> (n1-n2+1);
  nat lr= aligned_size<C,V> (n1);
  C* q= mmx_formatted_new<C> (lq, CF(P1));
  C* r= mmx_formatted_new<C> (lr, CF(P1));
  Pol::copy (r, seg (P1), n1);
  Pol::quo_rem (q, r, seg (P2), n1, n2);
  Q= Polynomial (q, n1-n2+1, lq, CF(P1));
  return Polynomial (r, n2 - 1, lr, CF(P1));
}

TMPL Polynomial
trem (int d1, const Polynomial& P1, const Polynomial& P2) {
  typedef implementation<polynomial_divide,V> Pol;
  nat n1= max (0, d1+1), n= N(P1), n2= N(P2);
  ASSERT (n <= n2-1, "bad dimension in trem");
  if (n1 < n2 || n2 == 0) return P1;
  nat l= aligned_size<C,V> (n1);
  C* q= mmx_formatted_new<C> (l, CF(P1));
  C* r= mmx_formatted_new<C> (l, CF(P1));
  Pol::copy  (r, seg (P1), n);
  Pol::clear (r+n, n1-n);
  Pol::tquo_rem (q, r, seg (P2), n1, n2);
  mmx_delete<C> (r, l);
  return Polynomial (q, n1, l, CF(P1));
}

TMPL inline int
pexponent (const Polynomial& P1, const Polynomial& P2) {
  return max (0, deg (P1) - deg (P2) + 1);
}

TMPL Polynomial
pquo (const Polynomial& P1, const Polynomial& P2) {
  typedef implementation<polynomial_divide,V> Pol;
  nat n1= N(P1), n2= N(P2);
  if (n1 < n2 || n2 == 0) return Polynomial (CF(P1));
  nat lq= aligned_size<C,V> (n1-n2+1);
  nat lr= aligned_size<C,V> (n1);
  C* q= mmx_formatted_new<C> (lq, CF(P1));
  C* r= mmx_formatted_new<C> (lr, CF(P1));
  Pol::copy (r, seg (P1), n1);
  Pol::pquo_rem (q, r, seg (P2), n1, n2);
  mmx_delete<C> (r, lr);
  return Polynomial (q, n1-n2+1, lq, CF(P1));
}

TMPL Polynomial
prem (const Polynomial& P1, const Polynomial& P2) {
  typedef implementation<polynomial_divide,V> Pol;
  nat n1= N(P1), n2= N(P2);
  if (n1 < n2 || n2 == 0) return P1;
  nat lq= aligned_size<C,V> (n1-n2+1);
  nat lr= aligned_size<C,V> (n1);
  C* q= mmx_formatted_new<C> (lq, CF(P1));
  C* r= mmx_formatted_new<C> (lr, CF(P1));
  Pol::copy (r, seg (P1), n1);
  Pol::pquo_rem (q, r, seg (P2), n1, n2);
  mmx_delete<C> (q, lq);
  return Polynomial (r, n2 - 1, lr, CF(P1));
}

TMPL Polynomial
prem (const Polynomial& P1, const Polynomial& P2, Polynomial& Q) {
  typedef implementation<polynomial_divide,V> Pol;
  nat n1= N(P1), n2= N(P2);
  if (n1 < n2 || n2 == 0) return Polynomial (CF(P1));
  nat lq= aligned_size<C,V> (n1-n2+1);
  nat lr= aligned_size<C,V> (n1);
  C* q= mmx_formatted_new<C> (lq, CF(P1));
  C* r= mmx_formatted_new<C> (lr, CF(P1));
  Pol::copy (r, seg (P1), n1);
  Pol::pquo_rem (q, r, seg (P2), n1, n2);
  Q= Polynomial (q, n1-n2+1, lq, CF(P1));
  return Polynomial (r, n2 - 1, lr, CF(P1));
}

TMPL inline Polynomial operator / (const C& c, const Polynomial& P) {
  return Polynomial (c) / P; }

ARITH_SCALAR_INT_SUGAR(TMPL,Polynomial)

/******************************************************************************
* Gcd computations
******************************************************************************/

// gcd (P1, P2) = U1 P1 + U2 P2, with
// deg U1 < deg P2 and deg U2 < deg P1 (for non-degenerated cases)

TMPL inline C
contents (const Polynomial& P) {
  typedef implementation<polynomial_linear,V> Pol;
  return Pol::big_gcd (seg (P), N(P));
}

TMPL inline Polynomial
primitive_part (const Polynomial& P, C& c) {
  typedef implementation<polynomial_linear,V> Pol;
  c= contents (P);
  nat n= N(P);
  nat l= aligned_size<C,V> (n);
  C* r= mmx_formatted_new<C> (l, CF(P));
  Pol::div_sc (r, seg (P), c, n);
  return Polynomial (r, n, l, CF(P));
}

TMPL inline Polynomial
primitive_part (const Polynomial& P) {
  C c;
  return primitive_part (P, c);
}

TMPL inline Polynomial
gcd (const Polynomial& P1, const Polynomial& P2,
     Polynomial& U1, Polynomial& U2) {
  typedef implementation<polynomial_gcd,V> Pol;
  C one= promote (1, CF(P1));
  U1= one; U2= one;
  return Pol::gcd (P1, P2, U1, U2);
}

TMPL inline Polynomial
gcd (const Polynomial& P1, const Polynomial& P2, Polynomial& U1) {
  Polynomial U2(0); U1= promote (1, CF(P1));
  return gcd (P1, P2, U1, U2);
}

TMPL inline Polynomial
gcd (const Polynomial& P1, const Polynomial& P2) {
  typedef implementation<polynomial_gcd,V> Pol;
  return Pol::gcd (P1, P2);
}

TMPL static vector<Polynomial>
gcd (const Polynomial& p, const vector<Polynomial>& q) {
  typedef implementation<polynomial_evaluate,V> Pol;
  return Pol::multi_gcd (p, q);
}

TMPL inline Polynomial
lcm (const Polynomial& P1, const Polynomial& P2) {
  typedef implementation<polynomial_gcd,V> Pol;
  return Pol::lcm (P1, P2);
}

TMPL inline Polynomial
invert_modulo (const Polynomial& P, const Polynomial& Q) {
  typedef implementation<polynomial_gcd,V> Pol;
  return Pol::invert_mod (P, Q);
}

TMPL inline void
pade (const Polynomial& P, nat n, nat k, Polynomial& Num, Polynomial &Den) {
  typedef implementation<polynomial_gcd,V> Pol;
  Pol::pade (P, n, k, Num, Den);
}

template<typename C,typename V,typename W> inline void
minimal_polynomial_bis (Polynomial& p, const vector<C,W>& v) {
  // p is the minimal polynomial of the sequence of the elements in v
  // Algorithm 12.19 of "Modern Computer Algebra"
  nat n= N(v); nat k= n >> 1;
  ASSERT ((n & 1) == 0, "even size expected");
  Polynomial h (v), s, t;
  pade (h, n, k, s, t);
  t= reverse (t);
  p= (deg(t) < 1 + deg(s)) ? lshiftz (t, 1 + deg (s) - deg (t)) : t;
}

template<typename C,typename W> inline
polynomial<C,typename Polynomial_variant(C) >
minimal_polynomial (const vector<C,W>& v) {
  polynomial<C,typename Polynomial_variant(C) > P;
  minimal_polynomial_bis (P, v);
  return P;
}

/******************************************************************************
* Subresultants
******************************************************************************/

// Recall that subresultants are only well-defined until
// degree min (deg (P), deg (Q)) - 1.

TMPL inline vector<Polynomial>
subresultants (const Polynomial& P, const Polynomial& Q,
               vector<Polynomial>& co_P, vector<Polynomial>& co_Q) {
  typedef implementation<polynomial_subresultant,V> Pol;
  int n= deg (P), m= deg (Q); nat l= max (min (n, m), 0);
  C one= promote (1, CF(P));
  Polynomial d; vector<Polynomial> res (Polynomial (one), l);
  co_P= vector<Polynomial> (Polynomial (one), l);
  co_Q= vector<Polynomial> (Polynomial (one), l); 
  Pol::subresultant_sequence (P, Q, res, co_P, co_Q,
                              d, d, d, d, d, d, 0);
  return res;
}

TMPL inline vector<Polynomial>
subresultants (const Polynomial& P, const Polynomial& Q,
               vector<Polynomial>& co_P) {
  typedef implementation<polynomial_subresultant,V> Pol;
  int n= deg (P), m= deg (Q); nat l= max (min (n, m), 0);
  Polynomial d;
  C zero= promote (0, CF(P)), one= promote (1, CF(P));
  vector<Polynomial> res (Polynomial (one), l), co_Q (Polynomial (zero), 0); 
  co_P= vector<Polynomial> (Polynomial (one), l);
  Pol::subresultant_sequence (P, Q, res, co_P, co_Q,
                              d, d, d, d, d, d, 0);
  return res;
}

TMPL inline vector<Polynomial>
subresultants (const Polynomial& P, const Polynomial& Q) {
  typedef implementation<polynomial_subresultant,V> Pol;
  int n= deg (P), m= deg (Q); nat l= max (min (n, m), 0);
  Polynomial d;
  C zero= promote (0, CF(P)), one= promote (1, CF(P));
  vector<Polynomial>
    res (Polynomial(one), l),
    co_P (Polynomial (zero), 0),
    co_Q (Polynomial (zero), 0); 
  Pol::subresultant_sequence (P, Q, res, co_P, co_Q,
                              d, d, d, d, d, d, 0);
  return res;
}

TMPL inline Polynomial
subresultant (const Polynomial& P, const Polynomial& Q, int k,
              Polynomial& coP, Polynomial& coQ) {
  typedef implementation<polynomial_subresultant,V> Pol;
  int n= deg (P), m= deg (Q);
  nat l= max (min (n, m), 0);
  ASSERT (k < l, "index out of range");
  C zero= promote (0, CF(P)), one= promote (1, CF(P));
  Polynomial d; 
  vector<Polynomial>
    res (Polynomial (zero), l),
    co_P (Polynomial (zero), l),
    co_Q (Polynomial (zero), l); 
  res [k]= Polynomial (one);
  co_P[k]= Polynomial (one);
  co_Q[k]= Polynomial (one);
  Pol::subresultant_sequence (P, Q, res, co_P, co_Q,
                              d, d, d, d, d, d, 0);
  coP= co_P[k]; coQ= co_Q[k];
  return res[k];
}

TMPL inline Polynomial
subresultant (const Polynomial& P, const Polynomial& Q, int k,
              Polynomial& coP) {
  typedef implementation<polynomial_subresultant,V> Pol;
  int n= deg (P), m= deg (Q);
  nat l= max (min (n, m), 0);
  ASSERT (k < l, "index out of range");
  Polynomial d;
  C zero= promote (0, CF(P)), one= promote (1, CF(P));
  vector<Polynomial>
    res (Polynomial (zero), l),
    co_P (Polynomial (zero), l),
    co_Q (Polynomial (zero), 0); 
  res [k]= Polynomial (one);
  co_P[k]= Polynomial (one);
  Pol::subresultant_sequence (P, Q, res, co_P, co_Q,
                              d, d, d, d, d, d, 0);
  coP= co_P[k];
  return res[k];
}

TMPL inline Polynomial
subresultant (const Polynomial& P, const Polynomial& Q, int k) {
  typedef implementation<polynomial_subresultant,V> Pol;
  if (k < 0) return promote (0, P);
  int n= deg (P), m= deg (Q);
  nat l= max (min (n, m), 0);
  ASSERT ((nat) k < l, "index out of range");
  Polynomial d;
  C zero= promote (0, CF(P)), one= promote (1, CF(P));
  vector<Polynomial>
    res (Polynomial (zero), l),
    co_P (Polynomial (zero), 0),
    co_Q (Polynomial (zero), 0);
  res [k]= Polynomial (one);
  Pol::subresultant_sequence (P, Q, res, co_P, co_Q,
                              d, d, d, d, d, d, 0);
  return res[k];
}

TMPL inline C
resultant (const Polynomial& P, const Polynomial& Q) {
  typedef implementation<polynomial_subresultant,V> Pol;
  int n= degree (P), m= degree (Q);
  if (n < 0 || m < 0) return promote (0, CF(P));
  if (m == 0) return binpow (Q[0], n);
  if (n == 0) {
    C r= binpow (P[0], m); 
    return (n & 1) ? -r : r;
  }
  Polynomial d;
  C zero= promote (0, CF(P)), one= promote (1, CF(P));
  vector<Polynomial>
    res (Polynomial (one), 1),
    co_P (Polynomial (zero), 0),
    co_Q (Polynomial (zero), 0);
  Pol::subresultant_sequence (P, Q, res, co_P, co_Q,
                              d, d, d, d, d, d, 0);
  return res[0][0];
}

TMPL inline C
discriminant (const Polynomial& P) {
  typedef implementation<polynomial_subresultant,V> Pol;
  return resultant (derive (P), P);
}

/******************************************************************************
* Differentiation, integration and composition
******************************************************************************/

TMPL Polynomial
derive (const Polynomial& P) {
  typedef implementation<polynomial_linear,V> Pol;
  nat n= N(P);
  if (n <= 1) return promote (0, P);
  nat l= aligned_size<C,V> (n-1);
  C* r= mmx_formatted_new<C> (l, CF(P));
  Pol::derive (r, seg (P), n);
  return Polynomial (r, n-1, l, CF(P));
}

TMPL Polynomial
derive (const Polynomial& P, const nat& order) {
  typedef implementation<polynomial_linear,V> Pol;
  nat n= N(P);
  if (n <= order) return promote (0, P);
  nat l= aligned_size<C,V> (n-order);
  C* r= mmx_formatted_new<C> (l, CF(P));
  Pol::derive (r, seg (P), n, order);
  return Polynomial (r, n-order, l, CF(P));
}

TMPL Polynomial
xderive (const Polynomial& P) {
  typedef implementation<polynomial_linear,V> Pol;
  nat n= N(P);
  nat l= aligned_size<C,V> (n);
  C* r= mmx_formatted_new<C> (l, CF(P));
  Pol::xderive (r, seg (P), n);
  return Polynomial (r, n, l, CF(P));
}

TMPL Polynomial
integrate (const Polynomial& P) {
  typedef implementation<polynomial_linear,V> Pol;
  nat n= N(P);
  nat l= aligned_size<C,V> (n+1);
  C* r= mmx_formatted_new<C> (l, CF(P));
  Pol::integrate (r, seg (P), n);
  return Polynomial (r, n+1, l, CF(P));
}

template<typename C,typename V,typename K> Polynomial
compose (const Polynomial& P1, const polynomial<K,V>& P2) {
  typedef implementation<polynomial_compose,V> Pol;
  nat n1= N(P1), n2= N(P2);
  if (n1 <= 1) return P1;
  if (n2 == 0) return P1[0];
  nat n= (n1-1) * (n2-1) + 1;
  nat l= aligned_size<C,V> (n);
  C* r= mmx_formatted_new<C> (l, CF(P1));
  Pol::compose (r, seg (P1), seg (P2), n1, n2);
  return Polynomial (r, n, l, CF(P1));
}

TMPL Polynomial
shift (const Polynomial& P, const C& sh) {
  typedef implementation<polynomial_compose,V> Pol;
  nat n= N(P);
  if (n <= 1 || sh == 0) return P;
  nat l= aligned_size<C,V> (n);
  C* r= mmx_formatted_new<C> (l, CF(P));
  Pol::shift (r, seg (P), sh, n);
  return Polynomial (r, n, l, CF(P));
}

TMPL Polynomial
shift (const Polynomial& P, int i) {
  return shift (P, promote (i, CF(P)));
}

template<typename C,typename V,typename K> Polynomial
q_difference (const Polynomial& P, const K& q) {
  typedef implementation<polynomial_linear,V> Pol;
  nat n= N(P);
  if (n <= 1) return P;
  nat l= aligned_size<C,V> (n);
  C* r= mmx_formatted_new<C> (l, CF(P));
  Pol::q_difference (r, seg (P), q, n);
  return Polynomial (r, n, l, CF(P));
}

TMPL Polynomial
dilate (const Polynomial& P, nat p) {
  typedef implementation<polynomial_linear,V> Pol;
  if (p == 1) return P;
  nat n= N(P);
  if (n <= 1) return P;
  nat k= (n-1)*p + 1;
  nat l= aligned_size<C,V> (k);
  C* r= mmx_formatted_new<C> (l, CF(P));
  Pol::dilate (r, seg (P), p, n);
  return Polynomial (r, k, l, CF(P));
}

/******************************************************************************
* Evaluation
******************************************************************************/

TMPL inline vector<Polynomial>
rem (const Polynomial& p, const vector<Polynomial>& q) {
  typedef implementation<polynomial_evaluate,V> Pol;
  return Pol::multi_rem (p, q);
}

TMPL inline vector<Polynomial>
prem (const Polynomial& p, const vector<Polynomial>& q) {
  typedef implementation<polynomial_evaluate,V> Pol;
  return Pol::multi_prem (p, q);
}

TMPL inline Polynomial
annulator_bis (const Vector& x) {
  typedef implementation<polynomial_evaluate,V> Pol;
  return Pol::template annulator<Polynomial> (x);
}

template<typename C> inline polynomial<C,typename Polynomial_variant(C) >
annulator (const Vector& x) {
  return annulator_bis<C,typename Polynomial_variant(C) > (x);
}

TMPL inline C
evaluate (const Polynomial& p, const C& x) {
  typedef implementation<polynomial_evaluate,V> Pol;
  return Pol::evaluate (seg (p), x, N(p));
}

TMPL inline C
eval (const Polynomial& p, const C& x) {
  return evaluate (p, x);
}

TMPL inline Polynomial
tevaluate_bis (const C& v, const C& x, nat n) {
  typedef implementation<polynomial_evaluate,V> Pol;
  nat l= aligned_size<C,V> (n);
  C* buf= mmx_formatted_new<C> (l, get_format (v));
  Pol::tevaluate (v, buf, x, n);
  return Polynomial (buf, n, l, get_format (v));
}

template<typename C> inline polynomial<C,typename Polynomial_variant(C) >
tevaluate (const C& v, const C& x, nat n) {
  return tevaluate_bis<C,typename Polynomial_variant(C) > (v, x, n);
}

TMPL inline Vector
evaluate (const Polynomial& p, const Vector& x) {
  typedef implementation<polynomial_evaluate,V> Pol;
  return Pol::evaluate (p, x);
}

TMPL inline Vector
eval (const Polynomial& p, const Vector& x) {
  return evaluate (p, x);
}

TMPL inline Polynomial
tevaluate_bis (const Vector& v, const Vector& x, nat l) {
  typedef implementation<polynomial_evaluate,V> Pol;
  return Pol::template tevaluate<Polynomial> (v, x, l);
}

template<typename C> inline polynomial<C,typename Polynomial_variant(C) >
tevaluate (const Vector& v, const Vector& x, nat l) {
  return tevaluate_bis<C,typename Polynomial_variant(C) > (v, x, l);
}

TMPL inline Polynomial
interpolate_bis (const Vector& v, const Vector& x) {
  typedef implementation<polynomial_evaluate,V> Pol;
  return Pol::template interpolate<Polynomial> (v, x);
}

template<typename C> inline polynomial<C,typename Polynomial_variant(C) >
interpolate (const Vector& v, const Vector& x) {
  return interpolate_bis<C,typename Polynomial_variant(C) > (v, x);
}

TMPL Vector
tinterpolate (const Polynomial& p, const Vector& x) {
  typedef implementation<polynomial_evaluate,V> Pol;
  return Pol::tinterpolate (p, x);
}

TMPL vector<Polynomial >
expand (const Polynomial& p, const Vector& v, const vector<nat>& mu) {
  typedef implementation<polynomial_evaluate,V> Pol;
  nat k= N(v);
  ASSERT (N(mu) == k, "dimensions don't match");
  C** r= mmx_new<C*> (k);
  for (nat i=0; i<k; i++)
    r[i]= mmx_formatted_new<C> (aligned_size<C,V> (mu[i]), CF(p));
  Pol::expand (r, seg (p), seg (v), seg (mu), N (p), k);
  nat l= default_aligned_size<Polynomial > (k);
  Polynomial* ret= mmx_formatted_new<Polynomial > (l, get_format (p));
  for (nat i=0; i<k; i++)
    ret[i]= Polynomial (r[i], mu[i], aligned_size<C,V> (mu[i]), CF(p));
  mmx_delete<C*> (r, k);
  return vector<Polynomial > (ret, k, l);
}

/******************************************************************************
* Miscellaneous
******************************************************************************/

TMPL Polynomial
range (const Polynomial& P, nat start, nat end) {
  typedef implementation<polynomial_linear,V> Pol;
  nat l= aligned_size<C,V> (end-start);
  C* r= mmx_formatted_new<C> (l, CF(P));
  for (nat i=start; i<end; i++) r[i-start]= P[i];
  return Polynomial (r, end-start, l, CF(P));
}

TMPL Polynomial
extract_mod (const Polynomial& P, nat k, nat p) {
  typedef implementation<polynomial_linear,V> Pol;
  nat n= (N(P) - k + p - 1) / p;
  nat l= aligned_size<C,V> (n);
  C* r= mmx_formatted_new<C> (l, CF(P));
  for (nat i=0; i<n; i++) r[i]= P[i*p+k];
  return Polynomial (r, n, l, CF(P));
}

TMPL Polynomial
lshiftz (const Polynomial& P, const int& shift) {
  typedef implementation<polynomial_linear,V> Pol;
  if (shift == 0) return P;
  else if (shift <= -((int) N(P))) return promote (0, P);
  else if (shift < 0) return range (P, (nat) (-shift), N(P));
  else {
    nat n= N(P) + shift;
    nat l= aligned_size<C,V> (n);
    C* r= mmx_formatted_new<C> (l, CF(P));
    Pol::copy (r+shift, seg (P), N(P));
    return Polynomial (r, n, l, CF(P));
  }
}

TMPL Polynomial
reverse (const Polynomial& P) {
  typedef implementation<polynomial_linear,V> Pol;
  typedef implementation<vector_linear,V> Vec;
  nat n= N(P);
  nat l= aligned_size<C,V> (n);
  C* r= mmx_formatted_new<C> (l, CF(P));
  Vec::vec_reverse (r, seg (P), n);
  return Polynomial (r, n, l, CF(P));
}

TMPL Polynomial
graeffe (const Polynomial& P) {
  typedef implementation<polynomial_graeffe,V> Pol;
  nat n= N(P);
  nat l= aligned_size<C,V> (n);
  C* r= mmx_formatted_new<C> (l, CF(P));
  Pol::graeffe (r, seg (P), n);
  return Polynomial (r, n, l, CF(P));
}

TMPL Polynomial
invert_lo (const Polynomial& P, nat m) {
  typedef implementation<polynomial_divide,V> Pol;
  nat n= N(P);
  nat l= aligned_size<C,V> (m);
  C* r= mmx_formatted_new<C> (l, CF(P));
  if (n >= m)
    Pol::invert_lo (r, seg (P), m);
  else {
    C* t= mmx_formatted_new<C> (l, CF(P));
    Pol::copy (t, seg (P), n);
    Pol::clear (t + n, m - n);
    Pol::invert_lo (r, t, m);
  }
  return Polynomial (r, m, l, CF(P));
}

TMPL inline Polynomial
invert_lo (const Polynomial& P) {
  return invert_lo (P, N(P));
}

TMPL Polynomial
invert_hi (const Polynomial& P) {
  typedef implementation<polynomial_divide,V> Pol;
  nat n= N(P);
  nat l= aligned_size<C,V> (n);
  C* r= mmx_formatted_new<C> (l, CF(P));
  Pol::invert_hi (r, seg (P), n);
  return Polynomial (r, n, l, CF(P));
}

/******************************************************************************
* Mapping operations
******************************************************************************/

template<typename Op, typename C, typename V>
polynomial<Unary_return_type(Op,C),V>
unary_map (const polynomial<C,V>& p) {
  typedef implementation<vector_linear,V> Vec;
  typedef Unary_return_type(Op,C) T;
  nat n= N(p);
  nat l= aligned_size<T,V> (n);
  format<T> fm= unary_map<Op> (CF(p));
  T* r= mmx_formatted_new<T> (l, fm);
  Vec::template vec_unary<Op> (r, seg (p), n);
  return polynomial<T,V> (r, n, l, fm);
}

template<typename Op, typename C, typename V, typename X>
polynomial<Binary_return_type(Op,C,X),V>
binary_map_scalar (const polynomial<C,V>& p, const X& x) {
  typedef implementation<vector_linear,V> Vec;
  typedef Binary_return_type(Op,C,X) T;
  nat n= N(p);
  nat l= aligned_size<T,V> (n);
  format<T> fm= binary_map_scalar<Op> (CF(p), x);
  T* r= mmx_formatted_new<T> (l, fm);
  Vec::template vec_binary_scalar<Op> (r, seg (p), x, n);
  return polynomial<T,V> (r, n, l, fm);
}

template<typename Op, typename C, typename V> Unary_return_type(Op,C)
big (const polynomial<C,V>& v) {
  typedef implementation<vector_linear,V> Vec;
  return Vec::template vec_unary_big<Op> (seg (v), N(v));
}

TMPL Polynomial copy (const Polynomial& P) {
  return unary_map<id_op> (P); }
TMPL Polynomial duplicate (const Polynomial& P) {
  return unary_map<duplicate_op> (P); }

/******************************************************************************
* Dynamic mappers
******************************************************************************/

template<typename S1, typename D> polynomial<D>
map (const function_1<D,Argument(S1) >& fun,
     const polynomial<S1>& p1, const format<D>& fm)
{
  nat n= N(p1);
  nat l= default_aligned_size<D> (n);
  D* r= mmx_formatted_new<D> (l, fm);
  for (nat i=0; i<n; i++) r[i]= fun (p1[i]);
  return polynomial<D> (r, n, l, fm);
}

template<typename S1, typename D> inline polynomial<D>
map (D (*fun) (const S1&),
     const polynomial<S1>& p1, const format<D>& fm) {
  return map (function_1<D,Argument(S1) > (fun), p1, fm);
}

/******************************************************************************
* Floating point related functions
******************************************************************************/

TMPL inline bool is_finite (const Polynomial& p) {
  return big<and_is_finite_op> (p); }
TMPL inline bool is_nan (const Polynomial& p) {
  return big<or_is_nan_op> (p); }
TMPL inline bool is_infinite (const Polynomial& p) {
  return !is_nan (p) && big<or_is_infinite_op> (p); }
TMPL inline bool is_fuzz (const Polynomial& p) {
  return !is_nan (p) && big<or_is_fuzz_op> (p); }
TMPL inline bool is_reliable (const Polynomial& p) {
  return is_reliable (promote (0, CF(p))); }

TMPL inline Polynomial change_precision (const Polynomial& P, xnat p) {
  return binary_map_scalar<change_precision_op> (P, p); }
TMPL inline xnat precision (const Polynomial& p) {
  return big<min_precision_op> (p); }

TMPL inline xint exponent (const Polynomial& p) {
  return big<max_exponent_op> (p); }
TMPL inline double magnitude (const Polynomial& p) {
  return big<max_magnitude_op> (p); }

/******************************************************************************
* Complex and ball polynomials
******************************************************************************/

TMPL Abs_polynomial abs (const Polynomial& p) { return unary_map<abs_op> (p); }
TMPL Real_polynomial Re (const Polynomial& p) { return unary_map<Re_op> (p); }
TMPL Real_polynomial Im (const Polynomial& p) { return unary_map<Im_op> (p); }
TMPL Polynomial conj (const Polynomial& p) { return unary_map<conj_op> (p); }

TMPL Center_polynomial center (const Polynomial& p) {
  return unary_map<center_op> (p); }
TMPL Radius_polynomial radius (const Polynomial& p) {
  return unary_map<radius_op> (p); }
TMPL Center_polynomial lower (const Polynomial& p) {
  return unary_map<lower_op> (p); }
TMPL Center_polynomial upper (const Polynomial& p) {
  return unary_map<upper_op> (p); }
TMPL inline Polynomial sharpen (const Polynomial& p) {
  return unary_map<sharpen_op> (p); }
TMPLK inline Polynomial blur (const Polynomial& p, const K& x) {
  return binary_map_scalar<blur_op> (p, x); }

#undef TMPL_DEF
#undef TMPL
#undef TMPLK
#undef Format
#undef Vector
#undef Polynomial
#undef Polynomial_rep
#undef Abs_polynomial
#undef Real_polynomial
#undef Center_polynomial
#undef Radius_polynomial
} // namespace mmx
#endif // __MMX_POLYNOMIAL_HPP