include/algebramix/polynomial_naive.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : polynomial_naive.hpp
* DESCRIPTION: naive polynomial arithmetic over fields
* COPYRIGHT  : (C) 2003  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_NAIVE__HPP
#define __MMX__POLYNOMIAL_NAIVE__HPP
#include <basix/port.hpp>
#include <basix/vector.hpp>

namespace mmx {
#define TMPL  template <typename C>
#define TMPLK template <typename C, typename K>
#define TMPLX template <typename C, typename X>
#define TMPLP template <typename Polynomial>
#define Vector vector<C>
template <typename C, typename V> struct polynomial;

/******************************************************************************
* Naive variant for polynomials
******************************************************************************/

struct polynomial_naive {
  typedef vector_naive Vec;            // Polynomials as vectors
  typedef polynomial_naive Naive;      // Naive variant
  typedef polynomial_naive Positive;   // Multiplication without subtractions
  typedef polynomial_naive No_simd;    // Variant without SIMD instructions
  typedef polynomial_naive No_thread;  // Variant without threads
  typedef polynomial_naive No_scaled;  // Variant without preconditioning
};

template<typename C>
struct polynomial_variant_helper {
  typedef polynomial_naive PV;
};

/******************************************************************************
* Variable name storage, default is R[x]
******************************************************************************/

struct polynomial_defaults {};

template<typename V>
struct implementation<polynomial_defaults,V,polynomial_naive> {
  template<typename P>
  class global_variables {
    static inline generic& dyn_name () {
      static generic name = "x";
      return name; }
  public:
    static inline void set_variable_name (const generic& x) { dyn_name () = x; }
    static inline generic get_variable_name () { return dyn_name (); }
  };

}; // implementation<polynomial_defaults,V,polynomial_naive>

/******************************************************************************
* Inheritance from vectors
******************************************************************************/

template<typename V>
struct implementation<vector_defaults,V,polynomial_naive>:
  public implementation<vector_defaults,V,typename V::Vec> {};

template<typename V>
struct implementation<vector_allocate,V,polynomial_naive>:
  public implementation<vector_allocate,V,typename V::Vec> {};

template<typename V>
struct implementation<vector_abstractions,V,polynomial_naive>:
  public implementation<vector_abstractions,V,typename V::Vec> {};

template<typename V>
struct implementation<vector_abstractions_stride,V,polynomial_naive>:
  public implementation<vector_abstractions_stride,V,typename V::Vec> {};

template<typename V>
struct implementation<vector_linear,V,polynomial_naive>:
  public implementation<vector_linear,V,typename V::Vec> {};

/******************************************************************************
* Vectorial operations
******************************************************************************/

struct polynomial_vectorial {};

template<typename V>
struct implementation<polynomial_vectorial,V,polynomial_naive>:
  public implementation<polynomial_defaults,V>
{
  typedef implementation<vector_linear,V> Vec;

TMPLX static inline void set (C* dest, const X& c, nat n) {
  Vec::set (dest, c, n); }
TMPL static inline void clear (C* dest, nat n) {
  Vec::clear (dest, n); }
TMPL static inline void copy (C* dest, const C* s, nat n) {
  Vec::copy (dest, s, n); }
TMPLX static inline void cast (C* dest, const X* s, nat n) {
  Vec::cast (dest, s, n); }
TMPL static inline void reverse (C* dest, const C* s, nat n) {
  Vec::vec_reverse (dest, s, n); }
TMPL static inline void neg (C* dest, const C* s, nat n) {
  Vec::neg (dest, s, n); }
TMPL static inline void neg (C* dest, nat n) {
  Vec::neg (dest, n); }
TMPL static inline void add (C* dest, const C* s1, const C* s2, nat n) {
  Vec::add (dest, s1, s2, n); }
TMPL static inline void add (C* dest, const C* s, nat n) {
  Vec::add (dest, s, n); }
TMPL static inline void sub (C* dest, const C* s1, const C* s2, nat n) {
  Vec::sub (dest, s1, s2, n); }
TMPL static inline void sub (C* dest, const C* s, nat n) {
  Vec::sub (dest, s, n); }
TMPL static inline void mul_add (C* dest, const C* s1, const C* s2, nat n) {
  Vec::mul_add (dest, s1, s2, n); }
TMPLX static inline void mul_add (C* dest, const C* s, const X& c, nat n) {
  Vec::mul_add (dest, s, c, n); }
TMPLX static inline void mul_add (C* dest, const X& c, const C* s, nat n) {
  Vec::mul_add (dest, c, s, n); }
TMPL static inline void mul_sc (C* dest, const C* s, const C& c, nat n) {
  Vec::mul (dest, s, c, n); }
TMPL static inline void mul_sc (C* dest, const C& c, nat n) {
  Vec::mul (dest, c, n); }
TMPL static inline void div_sc (C* dest, const C* s, const C& c, nat n) {
  Vec::div (dest, s, c, n); }
TMPL static inline void div_sc (C* dest, const C& c, nat n) {
  Vec::div (dest, c, n); }
TMPL static inline void quo_sc (C* dest, const C* s, const C& c, nat n) {
  Vec::quo (dest, s, c, n); }
TMPL static inline void quo_sc (C* dest, const C& c, nat n) {
  Vec::quo (dest, c, n); }
TMPL static inline void rem_sc (C* dest, const C* s, const C& c, nat n) {
  Vec::rem (dest, s, c, n); }
TMPL static inline void rem_sc (C* dest, const C& c, nat n) {
  Vec::rem (dest, c, n); }
TMPL static inline bool equal (const C* s1, const C* s2, nat n) {
  return Vec::equal (s1, s2, n); }
TMPL static inline bool equal (const C* s, const C& c, nat n) {
  return Vec::equal (s, c, n); }
TMPL static inline bool exact_eq (const C* s1, const C* s2, nat n) {
  return Vec::exact_eq (s1, s2, n); }
TMPL static inline bool exact_eq (const C* s, const C& c, nat n) {
  return Vec::exact_eq (s, c, n); }
TMPL static inline void trim (const C* s, nat& n) {
  while (n > 0 && s[n-1] == 0) n--; }

}; // implementation<polynomial_vectorial,V,polynomial_naive>

/******************************************************************************
* Other linear time operators
******************************************************************************/

struct polynomial_linear {};

template<typename V>
struct implementation<polynomial_linear,V,polynomial_naive>:
  public implementation<polynomial_vectorial,V>
{
  typedef implementation<vector_linear,V> Vec;
  typedef implementation<polynomial_vectorial,V> Pol;

TMPL static inline int
sign (C* s, nat n) {
  while (n != 0) {
    int sgn= sign (*s);
    if (sgn != 0) return sgn;
    s++; n--;
  }
  return 0;
}

TMPL static inline int
cmp (const C* s1, const C* s2, nat n) {
  while (n != 0) {
    int sgn= compare (*s1, *s2);
    if (sgn != 0) return sgn;
    s1++; s2++; n--;
  }
  return 0;
}

TMPL static inline void
derive (C* dest, nat n) {
  for (nat i=1; i<n; i++)
    dest[i-1]= promote ((int) i, dest[i]) * dest[i];
}

TMPL static inline void
derive (C* dest, nat n, nat order) {
  if (order == 0) return;
  for (nat i=order; i<n; i++) {
    dest[i-order]= promote (i, dest[i]) * dest[i];
    for (nat j=1; j<order; j++)
      dest[i-order] *= promote ((int) (i-j), dest[i]);
  }
}

TMPL static inline void
derive (C* dest, const C* s, nat n) {
  for (nat i=1; i<n; i++)
    dest[i-1]= promote ((int) i, s[i]) * s[i];
}

TMPL static inline void
derive (C* dest, const C* s, nat n, nat order) {
  for (nat i=order; i<n; i++) {
    dest[i-order]= s[i];
    for (nat j=0; j<order; j++)
      dest[i-order] *= promote ((int) (i-j), s[i]);
  }
}

TMPL static inline void
xderive (C* dest, nat n) {
  for (nat i=0; i<n; i++)
    dest[i]= promote ((int) i, dest[i]) * dest[i];
}

TMPL static inline void
xderive (C* dest, const C* s, nat n) {
  for (nat i=0; i<n; i++)
    dest[i]= promote ((int) i, s[i]) * s[i];
}

TMPL static inline void
integrate (C* dest, const C* s, nat n) {
  dest[0]= promote (0, dest[0]);
  for (nat i=0; i<n; i++)
    dest[i+1]= s[i] / promote ((int) (i+1), s[i]);
}

TMPLK static inline void
q_difference (C* dest, const C* s, const K& q, nat n) {
  K p= promote (1, q);
  while (n != 0) {
    *dest= p * (*s);
    dest++; s++; n--; p *= q;
  }
}

TMPL static inline void
dilate (C* dest, const C* s, nat p, nat n) {
  Vec::clear (dest, (n-1)*p + 1);
  while (n != 0) {
    *dest= *s;
    dest += p; s++; n--;
  }
}

TMPL inline static C
big_gcd (const C* s, nat n) {
  return Vec::template vec_unary_big<gcd_op, C> (s, n);
}

}; // implementation<polynomial_linear,V,polynomial_naive>

/******************************************************************************
* Multiplication
******************************************************************************/

struct polynomial_multiply {};

template<typename V>
struct implementation<polynomial_multiply,V,polynomial_naive>:
  public implementation<polynomial_linear,V>
{
  typedef implementation<polynomial_linear,V> Pol;
  typedef implementation<vector_linear,V> Vec;

TMPLK static void
mul (C* dest, const C* s1, const K* s2, nat n1, nat n2) {
  if (n1+n2>0)
    Pol::clear (dest, n1+n2-1);
  while (n2 != 0) {
    Pol::mul_add (dest, s1, *s2, n1);
    dest++; s2++; n2--;
  }
}

TMPL static void
square (C* dest, const C* s, nat n) {
  if (n>0)
    Pol::clear (dest, n+n-1);
  else return;
  for (nat i=0; i<n-1; i++)
    Pol::mul_add (dest+i+i+1, s+i+1, s[i], n-i-1);
  for (nat i=1; i<n+n-2; i++)
    dest[i] = dest[i] + dest[i];
  for (nat i=0; i<n; i++)
    dest[i+i] += s[i] * s[i];
}

TMPL static void
tmul (C* dest, const C* s1, const C* s2, nat n1, nat n2) {
  // transposed multiplication
  // s1 has length n1, dest has length n2 and s2 has length n1+n2-1
  Pol::clear (dest, n2);
  while (n1 != 0) {
    Pol::mul_add (dest, *s1, s2, n2);
    s1++; s2++; n1--;
  }
}

}; // implementation<polynomial_multiply,V,polynomial_naive>

/******************************************************************************
* Graeffe transforms
******************************************************************************/

struct polynomial_graeffe {};

template<typename V>
struct implementation<polynomial_graeffe,V,polynomial_naive>:
  public implementation<polynomial_multiply,V>
{
  typedef implementation<polynomial_multiply,V> Pol;
  typedef implementation<vector_linear,V> Vec;

TMPL static void
graeffe (C* dest, const C* s, nat n) {
  if (n <= 1) return;
  nat l1= (n+1) >> 1, l2= n >> 1;
  nat l= aligned_size<C,V> (n << 1);
  C* c1 = mmx_new<C> (l);
  C* c2 = c1 + l1;
  C* aux= c2 + l2;
  if (n>0) dest[n-1]= C(0);
  Vec::half_copy (c1, s, l1);
  Vec::half_copy (c2, s+1, l2);
  Pol::square (dest, c1, l1);
  Pol::square (aux , c2, l2);
  Pol::sub (dest+1, aux, (l2<<1) - 1);
  mmx_delete<C> (c1, l);
}

}; // implementation<polynomial_graeffe,V,polynomial_naive>

/******************************************************************************
* Division
******************************************************************************/

struct polynomial_divide {};

template<typename V>
struct implementation<polynomial_divide,V,polynomial_naive>:
  public implementation<polynomial_multiply,V>
{
  typedef implementation<polynomial_multiply,V> Pol;
  typedef implementation<vector_linear,V> Vec;

TMPL static void
quo_rem (C* dest, C* s1, const C* s2, nat n1, nat n2) {
  // (s1, n1) contains the numerator on input and the remainder on output
  // (s2, n2) contains the denominator. We assume n2>0 and s2[n2-1] != 0
  // (dest, n1-n2+1) contains the quotient
  while (n1 >= n2 && n1 != 0) {
    int d= n1-n2;
    C q= quo (s1[n1-1], s2[n2-1]);
    dest[d]= q;
    for (nat i=0; i<n2; i++)
      s1[i+d] -= q * s2[i];
    n1--;
  }
}

TMPL static void
tquo_rem (C* dest, const C* s1, const C* s2, nat n1, nat n2) {
  // (s1, n2-1) contains the part corresponding to the remainder, while
  // (s2+n2-1, n1-n2+1) contains the one of the quotient.
  // We assume n2>0 and s2[n2-1] != 0.
  // (dest, n1) contains the output.
  C c = s2[n2-1];
  Pol::copy (dest, s1, n2-1);
  s1 += n2-1;
  for (nat i = n2-1; i < n1; i++, dest++, s1++)
    dest[n2-1] = quo (*s1 - Vec::inn_prod (dest, s2, n2-1), c);
}

TMPL static void
pquo_rem (C* dest, C* s1, const C* s2, nat n1, nat n2) {
  // (s1, n1) contains the numerator on input
  //          and the pseudo-remainder on output
  // (s2, n2) contains the denominator.
  // We assume n2>0 and s2[n2-1] != 0
  // (dest, n1-n2+1) contains the pseudo-quotient
  nat l= n1 - n2 + 1;
  while (n1 >= n2 && n1 != 0) {
    int d= n1-n2;
    dest[d]= s1[n1-1];
    for (nat i=0; i<n1-1; i++)
      s1[i]*= s2[n2-1];
    for (nat i=0; i<n2-1; i++)
      s1[i+d]-= s1[n1-1] * s2[i];
    s1[n1-1]=0;
    n1--;
  }
  C c= s2[n2-1];
  for (nat i = 1; i < l; i++) {
    dest[i] *= c;
    c *= s2[n2-1];
  }
}

TMPL static void
invert_hi (C* dest, const C* src, nat n) {
  if (n == 1) *dest= C(1) / *src;
  else {
    nat h= (n+1) >> 1;
    nat l= n - h;
    invert_hi (dest + l, src + l, h);
    nat buf_size= aligned_size<C,V> (n << 1);
    C* buf= mmx_new<C> (buf_size);
    C* aux= buf + l;
    Pol::mul (aux, src, dest + l, n, h);
      // FIXME: use middle product
    Pol::mul (buf, dest + h, aux + h - 1, l, l);
      // FIXME: use truncated product
    Pol::neg (dest, buf + l - 1, l);
    mmx_delete<C> (buf, buf_size);
  }
}

}; // implementation<polynomial_divide,V,polynomial_naive>

/******************************************************************************
* Exact division
******************************************************************************/

struct polynomial_exact_divide {};

template<typename V>
struct polynomial_exact_divide_threshold {};

template<typename V>
struct implementation<polynomial_exact_divide,V,polynomial_naive>:
  public implementation<polynomial_divide,V>
{
  typedef implementation<polynomial_divide,V> Pol;
  typedef implementation<vector_linear,V> Vec;

TMPL static void
div (C* dest, const C* s1, const C* s2, nat n1, nat n2) {
  // (s1, n1) contains the numerator on input
  // (s2, n2) contains the denominator. We assume n2>0 and s2[n2-1] != 0
  // (dest, n1-n2+1) contains the quotient
  nat l1= aligned_size<C,V> (n1);
  C* r= mmx_new<C> (l1);
  Vec::copy (r, s1, n1);
  Pol::quo_rem (dest, r, s2, n1, n2);
  if (! Pol::equal (r, C(0), n1)) {
    mmx_delete<C> (r, l1);
    ERROR ("division must be exact");
  }
  mmx_delete<C> (r, l1);
}

}; // implementation<polynomial_exact_divide,V,polynomial_naive>:

/******************************************************************************
* Composition
******************************************************************************/

struct polynomial_compose {};

template<typename V>
struct implementation<polynomial_compose,V,polynomial_naive>:
  // FIXME: really a dichotomic algorithm; also implement naive version
  public implementation<polynomial_multiply,V>
{
  typedef implementation<polynomial_multiply,V> Pol;

TMPLK static void
compose (C* dest, K* p, const C* s1, const K* s2, nat n1, nat n2) {
  // computes compose ((s1, n1), (s2, n2)) and pow ((s2, n2), n1)
  if (n1 == 1) {
    *dest= *s1;
    Pol::copy (p, s2, n2);
  }
  else {
    nat d2= n2 - 1;
    nat h1= n1 >> 1, h2= n1 - h1;
    nat l1= (h1-1) * d2 + 1, l2= (h2-1) * d2 + 1;
    nat k1= l1 + d2, k2= l2 + d2;
    nat alloc_c1= aligned_size<C,V> (l1 + l2);
    C* c1= mmx_new<C> (alloc_c1);
    C* c2= c1 + l1;
    nat alloc_p1= aligned_size<K,V> (k1 + k2);
    C* p1= mmx_new<K> (alloc_p1);
    K* p2= p1 + k1;
    compose (c1, p1, s1, s2, h1, n2);
    compose (c2, p2, s1+h1, s2, h2, n2);
    Pol::mul (dest, c2, p1, l2, k1);
    Pol::add (dest, c1, l1);
    Pol::mul (p, p1, p2, k1, k2);
    mmx_delete<K> (p1, alloc_p1);
    mmx_delete<C> (c1, alloc_c1);
  }
}

TMPLK static void
compose (C* dest, const C* s1, const K* s2, nat n1, nat n2) {
  ASSERT (n1 > 0 && n2 > 0, "composition of empty segments");
  if (n1 == 1) *dest= *s1;
  else {
    nat d2= n2 - 1;
    nat h1= n1 >> 1, h2= n1 - h1;
    nat l1= (h1-1) * d2 + 1, l2= (h2-1) * d2 + 1;
    nat k1= l1 + d2;
    nat alloc_c1= aligned_size<C,V> (l1 + l2);
    C* c1= mmx_new<C> (alloc_c1);
    C* c2= c1 + l1;
    nat alloc_p1= aligned_size<K,V> (k1);
    K* p1= mmx_new<K> (alloc_p1);
    compose (c1, p1, s1, s2, h1, n2);
    compose (c2, s1+h1, s2, h2, n2);
    Pol::mul (dest, c2, p1, l2, k1);
    Pol::add (dest, c1, l1);
    mmx_delete<K> (p1, alloc_p1);
    mmx_delete<C> (c1, alloc_c1);
  }
}

TMPLK static void
shift (C* dest, const C* s, const K& sh, nat n) {
  if (n <= 1 || sh == 0) Pol::copy (dest, s, n);
  else {
    K q[2]; q[0]= sh; q[1]= promote (1, sh);
    compose (dest, s, q, n, 2);
  }
}

}; // implementation<polynomial_compose,V,polynomial_naive>

/******************************************************************************
* Euclidean sequence
******************************************************************************/

struct polynomial_euclidean {};

template<typename V>
struct implementation<polynomial_euclidean,V,polynomial_naive>:
  public implementation<polynomial_divide,V>
{
  typedef implementation<polynomial_divide,V> Pol;
  typedef implementation<vector_linear,V> Vec;

private:
// Cf "Modern Computer Algebra", Chap. 6.
// d1 and d2 contains the given polynomials at input, with n1 >= n2.
// At the end, we have d1 = u1 s1 + u2 s2,  d2 = v1 s1 + v2 s2, being
// the Bezout relations where the algorithm stops.

// The normalized remainder of size i is stored in r[i].
// Defected cases are also stored as for a subresultant sequence.
// The respective quotient is stored in q.
// The normalization coefficient is stored in rho[i].

// n contains the sequence of the sizes of the successive remainders.

TMPL static void
euclidean_sequence (C*& d1, C*& d2, nat& n1 , nat& n2,
                    C*& u1, C*& u2, nat& nu1, nat& nu2,
                    C*& v1, C*& v2, nat& nv1, nat& nv2,
                    nat* n, C* rho, C* q, C** r, C** co1, C** co2,
                    nat k) {
  // d1 and d2 must have size n1 >= n2.
  // n, rho, q, r, co1, co2 must have size n2.
  // u1 and v1 (resp. u2 and v2) must have size n2 (resp. n1).
  // The algorithm stops with the last remainder d1 of size bigger than k,
  // and d2 of size at most k.
  VERIFY (n1 >= n2, "bug");
  VERIFY (n1 > k, "bug");
  if (n != NULL) Vec::clear (n, n2);
  C c1= invert (d1[n1-1]), c2= invert (d2[n2-1]);
  Vec::mul (d1, c1, n1); Vec::mul (d2, c2, n2);
  if (rho != NULL) Pol::clear (rho, n2);
  if (u1 != NULL) {
    VERIFY (v1 != NULL, "bug");
    Pol::clear (u1, n2); u1[0]= c1; nu1= 1;
    Pol::clear (v1, n2);            nv1= 0; }
  if (u2 != NULL) {
    VERIFY (v2 != NULL, "bug");
    Pol::clear (u2, n1);            nu2= 0;
    Pol::clear (v2, n1); v2[0]= c2; nv2= 1; }
  nat i= 0;
  while (n2 > k)
    if (n1 >= n2) {
      C m= - d1[n1-1]; nat s= n1-n2;
      if (q != NULL) q[n1-1]= m;
      Vec::mul_add (d1 + s, m, d2, n2); Pol::trim (d1, n1);
      if (u1 != NULL && nv1 != 0) {
        Vec::mul_add (u1 + s, m, v1, nv1);
        nu1 = max (nu1, nv1 + s); Pol::trim (u1, nu1);
      }
      if (u2 != NULL && nv2 != 0 ) {
        Vec::mul_add (u2 + s, m, v2, nv2);
        nu2 = max (nu2, nv2 + s); Pol::trim (u2, nu2);
      }
    }
    else {
      if (n != NULL) n[i]= n1; i++;
      if (rho != NULL)
        rho[n1]= rho[n2-1]= n1 == 0 ? promote (1, c1) : d1[n1-1];
      c1= invert (n1 == 0 ? promote (1, c1) : d1[n1-1]);
      Vec::mul (d1, c1, n1);
      if (u1 != NULL) Vec::mul (u1, c1, nu1);
      if (u2 != NULL) Vec::mul (u2, c1, nu2);
      if (r != NULL && r[n1] != NULL)
        Pol::copy (r[n1], d1, n1);
      if (n1 + 1 != n2 && r != NULL && r[n2-1] != NULL)
        Pol::copy (r[n2-1], d1, n1);
      if (co1 != NULL && co1[n1] != NULL)
        Pol::copy (co1[n1], u1, nu1);
      if (n1 + 1 != n2 && co1 != NULL && co1[n2-1] != NULL)
        Pol::copy (co1[n2-1], u1, nu1);
      if (co2 != NULL && co2[n1] != NULL)
        Pol::copy (co2[n1], u2, nu2);
      if (n1 + 1 != n2 && co2 != NULL && co2[n2-1] != NULL)
        Pol::copy (co2[n2-1], u2, nu2);
      swap (d1, d2); swap (n1 , n2 );
      swap (u1, v1); swap (nu1, nv1);
      swap (u2, v2); swap (nu2, nv2);
    }
}

public:

TMPL static void
euclidean_sequence (const C* s1, const C* s2, nat n1, nat n2,
                    C* d1, C* d2, nat& m1 , nat& m2,
                    C* u1, C* u2, nat& nu1, nat& nu2,
                    C* v1, C* v2, nat& nv1, nat& nv2,
                    nat* n, C* rho, C* q, C** r, C** co1, C** co2, nat k= 0) {
  // Wrap the previous function, so that any argument becomes optional.
  VERIFY (n1 >= n2, "bug");
  nat l1= aligned_size<C,V> (n1), l2= aligned_size<C,V> (n2);
  nat l= max (l1, l2);
  C* t1= mmx_new<C> (l), * t2= mmx_new<C> (l);
  Pol::copy (t1, s1, n1); Pol::copy (t2, s2, n2);
  bool opt_co1= ! (u1 == NULL && v1 == NULL && co1 == NULL);
  bool opt_co2= ! (u2 == NULL && v2 == NULL && co2 == NULL);
  C* uu1= NULL, * uu2= NULL, * vv1= NULL, * vv2= NULL;
  if (opt_co1) { uu1= mmx_new<C> (l2); vv1= mmx_new<C> (l2); }
  if (opt_co2) { uu2= mmx_new<C> (l1); vv2= mmx_new<C> (l1); }
  nat mm1= n1, mm2= n2, mu1, mu2, mv1, mv2;
  // If the final remainders are not requested then we stop earlier:
  if (d1 == NULL && d2 == NULL && !opt_co1 && !opt_co2) {
    while (k < n2 &&
           (r   == NULL || r  [k] == NULL) &&
           (co1 == NULL || co1[k] == NULL) &&
           (co2 == NULL || co2[k] == NULL)) k++;
  }
  euclidean_sequence (t1 , t2 , mm1, mm2,
                      uu1, uu2, mu1, mu2,
                      vv1, vv2, mv1, mv2,
                      n, rho, q, r, co1, co2, k);
  if (d1 != NULL) { Pol::copy (d1, t1, mm1); m1= mm1; }
  if (d2 != NULL) { Pol::copy (d2, t2, mm2); m2= mm2; }
  mmx_delete<C> (t1, l); mmx_delete<C> (t2, l);
  if (u1 != NULL) { Pol::copy (u1, uu1, mu1); nu1= mu1; }
  if (u2 != NULL) { Pol::copy (u2, uu2, mu2); nu2= mu2; }
  if (v1 != NULL) { Pol::copy (v1, vv1, mv1); nv1= mv1; }
  if (v2 != NULL) { Pol::copy (v2, vv2, mv2); nv2= mv2; }
  if (opt_co1) { mmx_delete<C> (uu1, l2); mmx_delete<C> (vv1, l2); }
  if (opt_co2) { mmx_delete<C> (uu2, l1); mmx_delete<C> (vv2, l1); }
}
  
TMPL static void
gcd (C* g, nat& n, const C* s1, const C* s2, nat n1, nat n2,
     C* u1, C* u2, nat& nu1, nat& nu2) {
  VERIFY (n1>0 && n2>0 && s1[n1-1] != 0 && s2[n2-1] != 0,
          "invalid hypothesis for gcd computation");
  if (n1 < n2) return gcd (g, n, s2, s1, n2, n1, u2, u1, nu2, nu1);
  n= max (n1, n2); nat m= min (n1, n2);
  nat l1= aligned_size<C,V> (n1), l2= aligned_size<C,V> (n2);
  nat l= aligned_size<C,V> (n);
  C* d1= mmx_new<C> (l), * d2= mmx_new<C> (l);
  C* v1= u1 == NULL ? NULL : mmx_new<C> (l2);
  C* v2= u2 == NULL ? NULL : mmx_new<C> (l1);
  nat m1=0, m2=0, nv1=0, nv2=0;
  euclidean_sequence (s1, s2, n1, n2,
                      d1, d2, m1, m2,
                      u1, u2, nu1, nu2,
                      v1, v2, nv1, nv2,
                      (nat*)NULL, (C*)NULL, (C*)NULL,
                      (C**)NULL, (C**)NULL, (C**)NULL);
  Pol::copy (g, d1, m1); Pol::clear (g + m1, m - m1); n= m1;
  mmx_delete<C> (d1, l); mmx_delete<C> (d2, l);
  if (u1 != NULL) mmx_delete<C> (v1, l2);
  if (u2 != NULL) mmx_delete<C> (v2, l1);
}

TMPL static void
reconstruct (C* r, C* t, const C* s, nat m, const C* p, nat n, nat k) {
  // r (resp. t) must have allocated length at least k (resp. n-k+1)
  // s = r / t mod p, deg p = n, deg r < k, deg t <= n-k
  VERIFY (n > k && n > 0 && s[m-1]!= 0 && m <= n,
          "invalid hypothesis for reconstruction");
  nat n1= n+1, n2= m, l1= aligned_size<C,V> (n1);
  C* s1= mmx_new<C> (l1);
  Pol::copy (s1, p, n1);
  C* d1= NULL, * d2= mmx_new<C> (l1);
  C* u1= NULL, * u2= mmx_new<C> (l1);
  C* v1= NULL, * v2= mmx_new<C> (l1);
  nat nu1=0, nu2=0, nv1=0, nv2=0, nd1=0, nd2=0;
  euclidean_sequence (s1, s , n1, n2,
                      d1, d2, nd1, nd2,
                      u1, u2, nu1, nu2,
                      v1, v2, nv1, nv2,
                      (nat*)NULL, (C*)NULL, (C*)NULL,
                      (C**)NULL, (C**)NULL, (C**)NULL, k);
  VERIFY (nd2 <= k, "bug");
  VERIFY (nv2 <= n-k+1, "bug");
  Pol::copy (r, d2, nd2); Pol::clear (r + nd2, k - nd2);
  Pol::copy (t, v2, nv2); Pol::clear (t + nv2, n-k+1 - nv2);
  mmx_delete<C> (d2, l1);
  mmx_delete<C> (u2, l1); mmx_delete<C> (v2, l1);
}

}; // implementation<polynomial_euclidean,V,polynomial_naive>

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

struct polynomial_gcd {};

template<typename V>
struct implementation<polynomial_gcd,V,polynomial_naive>:
  public implementation<polynomial_divide,V>
{
  typedef implementation<polynomial_euclidean,V> Pol;
  typedef implementation<vector_linear,V> Vec;

#define C typename scalar_type_helper<Polynomial >::val
  
TMPLP static Polynomial
gcd (const Polynomial& P1, const Polynomial& P2,
     Polynomial& U1, Polynomial& U2) {
  bool opt1= U1 != 0, opt2= U2 != 0;
  nat n1= N(P1), n2= N(P2);
  C c1= P1[n1-1], c2= P2[n2-1];
  if (n1 == 0 && n2 == 0) {
    if (opt1) set_as (U1, 0);
    if (opt2) set_as (U2, 0);
    return promote (0, P1); }
  if (n2 == 0) {
    if (opt1) U1= promote (1, c1) / c1;
    if (opt2) set_as (U2, 0);
    return P1 / c1; }
  if (n1 == 0) {
    if (opt1) set_as (U1, 0);
    if (opt2) U2= promote (1, c2) / c2;
    return P2 / c2; }
  nat m = min (n1, n2), nu1, nu2;
  nat l = aligned_size<C,V> (m);
  nat l1= aligned_size<C,V> (n1), l2= aligned_size<C,V> (n2);
  C* tmp= mmx_new<C> (l);
  C* u1= opt1 ? mmx_new<C> (l2) : NULL;
  C* u2= opt2 ? mmx_new<C> (l1) : NULL;
  Pol::gcd (tmp, m, seg (P1), seg (P2), n1, n2, u1, u2, nu1, nu2);
  if (opt1) U1= Polynomial (u1, nu1, l2, CF(P1));
  if (opt2) U2= Polynomial (u2, nu2, l1, CF(P1));
  nat lg = aligned_size<C,V> (m);
  C* g= mmx_new<C> (lg);
  Pol::copy (g, tmp, m);
  mmx_delete<C> (tmp, l);
  return Polynomial (g, m, lg, CF(P1)); }

TMPLP static Polynomial
gcd (const Polynomial& P1, const Polynomial& P2) {
  Polynomial U1= promote (0, P1), U2= promote (0, P2);
  return gcd (P1, P2, U1, U2); }

TMPLP static inline Polynomial
lcm (const Polynomial& P1, const Polynomial& P2) {
  Polynomial U1= promote (0, P1), U2= promote (0, P1);
  return
    (P1 == 0 && P2 == 0) ?
    promote (0, P1) :
    P1 * (P2 / gcd (P1, P2, U1, U2)); }

TMPLP static Polynomial
invert_mod (const Polynomial& P, const Polynomial& Q) {
  // inverse of P modulo Q, return 0 if non-invertible
  Polynomial inv_P= C(1), dummy= C(0);
  Polynomial G= gcd (P, Q, inv_P, dummy);
  if (deg (G) != 0 || G[0] == 0) inv_P= 0;
  return inv_P / G[0]; }

TMPLP static void
reconstruct (const Polynomial& P, const Polynomial& Q, nat k,
             Polynomial& Num, Polynomial& Den) {
  // In return P = Num/Den mod Q, with deg (Q) = n,
  // deg (Num) < k and deg (Den) <= n-k.
  // Note that Num and Den are not ensured to be coprime.
  VERIFY (deg (Q) >= 0, "invalid argument");
  nat n= deg (Q);
  VERIFY (k <= n && k > 0, "invalid argument");
  nat m= N(P);
  if (n == 0 || m == 0) { Num= 0; Den= 1; return; }
  if (k == n) { Num= P; Den= 1; return; }
  nat lnum = aligned_size<C,V> (k), lden= aligned_size<C,V> (n-k+1);
  C* num= mmx_new<C> (lnum), * den= mmx_new<C> (lden);
  Pol::reconstruct (num, den, seg (P), m, seg (Q), n, k);
  Num= Polynomial (num, k, lnum, CF(P));
  Den= Polynomial (den, n-k+1, lden, CF(P)); }

#undef C
}; // implementation<polynomial_gcd,V,polynomial_naive>

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

struct polynomial_subresultant_base {};
struct polynomial_subresultant {};

template<typename V>
struct implementation<polynomial_subresultant_base,V,polynomial_naive>:
  public implementation<polynomial_divide,V>
{
  typedef implementation<polynomial_euclidean,V> Pol;
  typedef implementation<vector_linear,V> Vec;

TMPLP static void
subresultant_sequence (const Polynomial& f, const Polynomial& g,
                       vector<Polynomial>& res,
                       vector<Polynomial>& cof, vector<Polynomial>& cog,
                       Polynomial& d1, Polynomial& u1, Polynomial& v1,
                       Polynomial& d2, Polynomial& u2, Polynomial& v2,
                       nat dest_deg= 0) {
  /* In return res[k] contains the kth subresultant of f and g.
   * If res[k] is zero in input then this computation is skipped.
   * The same rule applies to the cofactors cof and cog, and d1, u1, etc
   * When all the following entries are actually filled we have:
   *   res[k] = cof[k] * f[k] + cog[k] * g[k]
   * d1 contains the smallest subresultant of degree at least dest_deg.
   * u1 and v1 are the corresponding cofactors: d1= u1 * f + v1 * g.
   * d2 contains the largest subresultant of degree less than dest_deg.
   * u2 and v2 are the corresponding cofactors: d2= u2 * f + v2 * g.
   */
  typedef typename scalar_type_helper<Polynomial>::val C;
  nat n1= N(f), n2= N(g);
  VERIFY (n1 >= n2, "bug");
  if (n2 <= 1) return;
  nat l1= aligned_size<C,V> (n1), l2= aligned_size<C,V> (n2);
  vector<C> rho (C(0), n2); C** sub= mmx_new<C*> (n2, (C*) NULL);
  vector<nat> ns ((nat) 0, n2);
  C* dd1= (d1 == 0 && u1 == 0 && v1 == 0) ? NULL : mmx_new<C> (l1);
  C* dd2= (d2 == 0 && u2 == 0 && v2 == 0) ? NULL : mmx_new<C> (l1);
  C* uu1= u1 == 0 ? NULL : mmx_new<C> (l2);
  C* vv1= v1 == 0 ? NULL : mmx_new<C> (l1);
  C* uu2= u2 == 0 ? NULL : mmx_new<C> (l2);
  C* vv2= v2 == 0 ? NULL : mmx_new<C> (l1);
  nat nd1= 0, nd2= 0, nu1= 0, nu2= 0, nv1= 0, nv2= 0;
  C** co1= mmx_new<C*> (n2, (C*) NULL); C** co2= mmx_new<C*> (n2, (C*) NULL);
  for (nat i = 1; i < n2; i++) {
    if (i-1 < N(res) && res[i-1] != 0)
      sub[i]= mmx_new<C> (aligned_size<C,V> (i), C(0));
    if (i-1 < N(cof) && cof[i-1] != 0)
      co1[i]= mmx_new<C> (l2, C(0));
    if (i-1 < N(cog) && cog[i-1] != 0)
      co2[i]= mmx_new<C> (l1, C(0));
  }
  Pol::euclidean_sequence (seg(f), seg(g), n1, n2,
                           dd1, dd2, nd1, nd2,
                           uu1, uu2, nu1, nu2,
                           vv1, vv2, nv1, nv2,
                           seg (ns), seg (rho), (C*) NULL,
                           sub, co1, co2, dest_deg);
  for (nat i = 1; i < n2; i++) {
    if (i-1 < N(res) && res[i-1] != 0)
      res[i-1]= Polynomial (sub[i], i, aligned_size<C,V> (i), CF(f));
    if (i-1 < N(cof) && cof[i-1] != 0)
      cof[i-1]= Polynomial (co1[i], n2, l2, CF(f));
    if (i-1 < N(cog) && cog[i-1] != 0)
      cog[i-1]= Polynomial (co2[i], n1, l1, CF(f));
  }
  if (d1 != 0) d1= Polynomial (dd1, nd1, l1, CF(f));
  if (d2 != 0) d2= Polynomial (dd2, nd2, l1, CF(f));
  if (u1 != 0) u1= Polynomial (uu1, nu1, l2, CF(f));
  if (u2 != 0) u2= Polynomial (uu2, nu2, l2, CF(f));
  if (v1 != 0) v1= Polynomial (vv1, nv1, l1, CF(f));
  if (v2 != 0) v2= Polynomial (vv2, nv2, l1, CF(f));
  mmx_delete<C*> (sub, n2);
  mmx_delete<C*> (co1, n2); mmx_delete<C*> (co2, n2);
  vector<C> c (promote (0, CF(f)), n2);
    // vector of the leading coefficients of the
    // subresultants (including the defective cases)
  nat n= n2 - 1, m= ns[0] - 1, i= 0;
  C c_n= g[n], c_n_1= f[n1-1] * binpow (g[n], n1 - n2 + 1) * rho[ns[0]];
  C d= binpow (g[n], n1 - n2);
  c[n-1]= c_n_1;
  while (m + 1 > 0) {
    C t0= binpow (c_n_1, n-m-1), t1= binpow (d, n-m-1);
    c[m]= (t0 / t1) * c_n_1;
    if (m == 0) break;
    t0 *= square_op::op (c_n_1); t1 *= ((n-m) & 1) ? c_n * d : -c_n * d;
    c[m-1]= (t0 * c_n * rho[m]) / t1;
    n= m; c_n= c[m]; d= c[m]; c_n_1= c[m-1]; i++; m= ns[i] - 1;
  }
  for (nat i= 0; i < N(res); i++) res[i] *= c[i];
  for (nat i= 0; i < N(cof); i++) cof[i] *= c[i];
  for (nat i= 0; i < N(cog); i++) cog[i] *= c[i];
  d= (nd1 > n2 || nd1 == 0) ? promote (0, d) : c[nd1 - 1];
  d1 *= d; u1 *= d; v1 *= d;
  d= (nd1 > n2 || nd1 <  2) ? promote (1, d) : c[nd1 - 2];
  d2 *= d; u2 *= d; v2 *= d; }

}; // implementation<polynomial_subresultant_base,V,polynomial_naive>

template<typename V>
struct implementation<polynomial_subresultant,V,polynomial_naive>:
  public implementation<polynomial_subresultant_base,V>
{
  // Common handle for when degree (f) < degree (g) or f or g is zero.
  typedef implementation<polynomial_subresultant_base,V> Pol;

TMPLP static void
subresultant_sequence (const Polynomial& f, const Polynomial& g,
                       vector<Polynomial>& res,
                       vector<Polynomial>& cof, vector<Polynomial>& cog,
                       Polynomial& d1, Polynomial& u1, Polynomial& v1,
                       Polynomial& d2, Polynomial& u2, Polynomial& v2,
                       nat dest_deg= 0) {
  int n= degree (f), m= degree (g);
  if (n <= 0 || m <= 0) {
    set_as (d1, 0); set_as (u1, 0); set_as (v1, 0);
    set_as (d2, 0); set_as (u2, 0); set_as (v2, 0);
    return;
  }
  if (n < m) {
    Pol::subresultant_sequence (g, f, res, cog, cof,
                                d1, v1, u1, d2, v2, u2, dest_deg);
    for (nat k= 0; k < (nat) n; k++)
      if ((m+n-1) & (n-k) & 1) {
        if (k < N(res)) res[k]= -res[k];
        if (k < N(cof)) cof[k]= -cof[k];
        if (k < N(cog)) cog[k]= -cog[k];
      }
    if ((m+n-1) & (n-degree (d1)) & 1) { d1= -d1; u1= -u1; v1= -v1; }
    else { d2= -d2; u2= -u2; v2= -v2; }
  }
  else
    Pol::subresultant_sequence (f, g, res, cof, cog,
                                d1, u1, v1, d2, u2, v2, dest_deg); }

}; // implementation<polynomial_subresultant,V,polynomial_naive>

/******************************************************************************
* Efficient evaluation of polynomials
******************************************************************************/

struct polynomial_evaluate {};

template<typename V>
struct implementation<polynomial_evaluate,V,polynomial_naive>:
  public implementation<polynomial_divide,V>
{
  typedef implementation<polynomial_divide,V> Pol;
  typedef implementation<vector_linear,V> Vec;

TMPL static void
multi_mod (C* r, const C* p, const C* q, const nat* d, nat n, nat k) {
  // The polynomial (p, n) is divided by k polynomials
  // (q, d[0]), (q + d[0], d[1]), ..., (q + d[0] + ... + d[k-2], d[k-1])
  // The remainders of these divisions are put into
  // (r, d[0]-1), (r + d[0], d[1]-1), ..., (r + d[0] + ... + d[k-2], d[k-1]-1)
  nat l= aligned_size<C,V> (n);
  C* tmp_q= mmx_new <C> (l);
  C* tmp_r= mmx_new <C> (l);
  for (nat i= 0; i < k; i++) {
    if (n < d[i]) {
      Pol::copy (r, p, n);
      Pol::clear (r+n, d[i]-1-n);
    }
    else {
      Pol::copy (tmp_r, p, n);
      Pol::quo_rem (tmp_q, tmp_r, q, n, d[i]);
      Pol::copy (r, tmp_r, d[i]-1);
    }
    r += d[i]; q += d[i];
  }
  mmx_delete<C> (tmp_q, l);
  mmx_delete<C> (tmp_r, l);
}

TMPL static void
annulator (C* d, const C* s, nat n) {
  // Given n constants s[0],...s[n-1],
  // construct the polynomial (z-s[0]) ... (z-s[n-1])
  switch (n) {
  case 0:
    d[0]= promote (1, d[0]);
    break;
  case 1:
    d[1]= promote (1, d[0]);
    d[0]= -s[0];
    break;
  case 2:
    d[2]= promote (1, d[0]);
    d[1]= -(s[0]+s[1]);
    d[0]= s[0]*s[1];
    break;
  case 3:
    d[3]= promote (1, d[0]);
    d[2]= -(s[0]+s[1]+s[2]);
    d[1]= s[0]*s[1] + s[1]*s[2]+ s[2]*s[0];
    d[0]= -s[0]*s[1]*s[2];
    break;
  default:
    {
      nat h= n>>1;
      nat l= aligned_size<C,V> (n+2);
      C* aux= mmx_new<C> (l);
      annulator (aux, s, h);
      annulator (aux+h+1, s+h, n-h);
      Pol::mul (d, aux, aux+h+1, h+1, n+1-h);
      mmx_delete<C> (aux, l);
    }
  }
}

TMPL static C
evaluate (const C* p, const C& x, nat l) {
  // Evaluate the polynomial (p, n) at x
  if (l == 0) return promote (0, x);
  C r= p[l-1];
  for (nat i=l-2; i!=((nat)(-1)); i--)
    r= (x*r) + p[i];
  return r;
}

TMPL void
tevaluate (const C& v, C* p, const C& x, nat l) {
  // Transposed evaluation of the polynomial (p, l) at x
  if (l == 0) return;
  p[0] = v;
  for (nat i= 1; i < l; i++) p[i] = x * p[i-1];
}

TMPL static void
evaluate (C* v, const C* p, const C* x, nat l, nat n) {
  // Evaluate the polynomial (p, l) at points x[0], ..., x[n-1]
  for (nat j=0; j<n; j++)
    if (l == 0) v[j]= 0;
    else {
      C r= p[l-1];
      for (nat i=l-2; i!=((nat)(-1)); i--)
        r= (x[j]*r) + p[i];
      v[j]= r;
    }
}

TMPL static void
tevaluate (const C* v, C* p, const C* x, nat l, nat n) {
  // Transposed evaluation of (p, l) at points x[0], ..., x[n-1]
  nat l_y= aligned_size<C,V> (n);
  C* y= mmx_new<C> (l_y);
  Pol::set (y, C(1), n);
  for (nat j= 0; j < l; j++) {
    p[j]= Vec::inn_prod (y, v, n);
    Vec::mul (y, x, n);
  }
  mmx_delete<C> (y, l_y);
}

TMPL static void
interpolate (C* p, const C* v, const C* x, nat n) {
  // Interpolate the polynomial (p, n) from values
  // v[0], ... , v[n-1] at points x[0], ..., x[n-1]
  if (n == 0) return;
  format<C> fm= get_format (v[0]);
  nat l= aligned_size<C,V> (n+1), i, j;
  C ** buf= mmx_new<C*> (n+1);
  for (i= 0; i < n+1; i++)
    buf[i]= mmx_formatted_new<C> (l, fm);
  C * tmp0= mmx_formatted_new<C> (l, fm);
  C * tmp1= mmx_formatted_new<C> (l, fm);
  buf[0][0]= promote (1, fm);
  for (i= 0; i < n; i++) {
    for (j= 0; j < i; j++) {
      Pol::copy (tmp0 + 1, buf[j], i); tmp0[0]= promote (0, fm);
      Pol::mul_sc (tmp1, buf[j], x[i], i); tmp1[i]= promote (0, fm);
      Pol::sub (buf[j], tmp0, tmp1, i+1);
    }
    Pol::copy (tmp0 + 1, buf[j], i + 1); tmp0[0]= promote (0, fm);
    Pol::mul_sc (tmp1, buf[j], x[i], i + 1); tmp1[i+1]= promote (0, fm);
    Pol::sub (buf[j+1], tmp0, tmp1, i + 2);
  }
  Pol::derive (tmp1, buf[n], n+1);
  evaluate (tmp0, tmp1, x, n, n);
  Vec::div (tmp1, v, tmp0, n);
  Pol::clear (p, n);
  for (i= 0; i < n; i++) {
    Pol::mul_sc (tmp0, buf[i], tmp1[i], n);
    Pol::add (p, tmp0, n);
  }
  mmx_delete<C> (tmp0, l); mmx_delete<C> (tmp1, l);
  for (i= 0; i < n+1; i++) mmx_delete<C> (buf[i], l);
  mmx_delete<C*> (buf, n+1);
}

TMPL static void
tinterpolate (const C* p, C* v, const C* x, nat n) {
  // Interpolate the polynomial (p, n) from values
  // v[0], ... , v[n-1] at points x[0], ..., x[n-1]
  if (n == 0) return;
  format<C> fm= get_format (v[0]);
  nat l= aligned_size<C,V> (n+1), i, j;
  C ** buf= mmx_new<C*> (n+1);
  for (i= 0; i < n+1; i++)
    buf[i]= mmx_formatted_new<C> (l, fm);
  C * tmp0= mmx_formatted_new<C> (l, fm);
  C * tmp1= mmx_formatted_new<C> (l, fm);
  buf[0][0]= promote (1, fm);
  for (i= 0; i < n; i++) {
    for (j= 0; j < i; j++) {
      Pol::copy (tmp0 + 1, buf[j], i); tmp0[0]= promote (0, fm);
      Pol::mul_sc (tmp1, buf[j], x[i], i); tmp1[i]= promote (0, fm);
      Pol::sub (buf[j], tmp0, tmp1, i+1);
    }
    Pol::copy (tmp0 + 1, buf[j], i + 1); tmp0[0]= promote (0, fm);
    Pol::mul_sc (tmp1, buf[j], x[i], i + 1); tmp1[i+1]= promote (0, fm);
    Pol::sub (buf[j+1], tmp0, tmp1, i + 2);
  }
  Pol::derive (tmp1, buf[n], n+1);
  evaluate (tmp0, tmp1, x, n, n);
  for (i= 0; i < n; i++)
    v[i]= Vec::inn_prod (buf[i], p, n);
  Vec::div (v, tmp0, n);
  mmx_delete<C> (tmp0, l); mmx_delete<C> (tmp1, l);
  for (i= 0; i < n+1; i++) mmx_delete<C> (buf[i], l);
  mmx_delete<C*> (buf, n+1);
}

TMPL static void
expand (C** v, const C* p, const C* x, const nat* nu, nat n, nat k) {
  (void) v; (void) p; (void) x; (void) nu; (void) n; (void) k;
  ERROR ("expand not yet implemented");
}

#define C typename scalar_type_helper<Polynomial>::val

TMPLP static vector<Polynomial>
multi_rem (const Polynomial& p, const vector<Polynomial>& q) {
  ASSERT (is_non_scalar (q), "non-scalar vector expected");
  nat k= N(q), n= N(p);
  if (k == 0) return q;
  vector<nat> d ((nat) 0, k);
  for (nat i= 0; i < k; i++) {
    d[i]= N(q[i]);
    ASSERT (d[i] > 0, "division by 0");
  }
  if (n == 0) return vector<Polynomial> (Polynomial (C(0)), k);
  nat m= big_add (d), l= aligned_size<C,V> (m);
  C* qq= mmx_new<C> (l), * r= mmx_new<C> (l), * tmp= qq;
  for (nat i= 0; i < k; i++) {
    Pol::copy (tmp, seg (q[i]), d[i]);
    if (q[i][d[i]-1] != 1)
      Pol::mul_sc (tmp, 1/q[i][d[i]-1], d[i]);
    tmp += d[i];
  }
  multi_mod (r, seg(p), qq, seg (d), n, k);

  mmx_delete<C> (qq, l);
  vector<Polynomial> ans (Polynomial (C(0)), k);  
  tmp= r;
  for (nat i= 0; i < k; i++) {
    nat l_t= aligned_size<C,V> (d[i]-1);
    C* t= mmx_new<C> (l_t);
    Pol::copy (t, tmp, d[i]-1);
    ans[i]= Polynomial (t, d[i]-1, l_t, CF(p));
    tmp += d[i];
  }
  mmx_delete<C> (r, l);
  return ans; }

TMPLP static vector<Polynomial>
multi_prem (const Polynomial& P, const vector<Polynomial>& Q) {
  return binary_map_scalar<prem_op> (Q, P); }

TMPLP static vector<Polynomial>
multi_gcd (const Polynomial& P, const vector<Polynomial>& Q) {
  return binary_map<gcd_op> (rem (P, Q), Q); }

TMPLP static Polynomial
annulator (const Vector& x) {
  ASSERT (is_non_scalar (x), "non-scalar vector expected");
  nat l= aligned_size<C,V> (N(x)+1);
  C* r= mmx_new<C> (l);
  annulator (r, seg(x), N(x));
  return Polynomial (r, N(x)+1, l, CF(x)); }

TMPLP static inline Vector
evaluate (const Polynomial& p, const Vector& x) {
  ASSERT (is_non_scalar (x), "non-scalar vector expected");
  nat l= default_aligned_size<C> (N(x));
  C* r= mmx_new<C> (l);
  evaluate (r, seg(p), seg(x), N(p), N(x));
  return Vector (r, N(x), l, CF(x)); }

TMPLP static inline Polynomial
tevaluate (const Vector& v, const Vector& x, nat l) {
  ASSERT (is_non_scalar (x), "non-scalar vector expected");
  nat l_r= aligned_size<C,V> (l), n= N(x);
  C* r= mmx_new<C> (l_r);
  tevaluate (seg(v), r, seg(x), l, n);
  return Polynomial (r, l, l_r, CF(v)); }

TMPLP static Polynomial
interpolate (const Vector& v, const Vector& x) {
  ASSERT (is_non_scalar (x), "non-scalar vector expected");
  ASSERT (N(v) == N(x), "dimensions don't match");
  nat n= N(v), l= aligned_size<C,V> (n);
  C* r= mmx_new<C> (l);
  interpolate (r, seg(v), seg(x), n);
  return Polynomial (r, n, l, CF(v)); }

TMPLP static Vector
tinterpolate (const Polynomial& p, const Vector& x) {
  ASSERT (is_non_scalar (x), "non-scalar vector expected");
  ASSERT (N(p) <= N(x), "dimensions don't match");
  nat n= N(x), lv = default_aligned_size<C> (n);
  nat lp= aligned_size<C,V> (n); 
  C* r= mmx_new<C> (lv); C* tmp= mmx_new<C> (lp);
  Pol::copy (tmp, seg(p), N(p)); Pol::clear (tmp + N(p), n - N(p)); 
  tinterpolate (tmp, r, seg(x), n);
  mmx_delete<C> (tmp, lp);
  return Vector (r, n, lv); }
#undef C

}; // implementation<polynomial_evaluate,V,polynomial_naive>

#undef TMPL
#undef TMPLK
#undef TMPLX
#undef TMPLP
} // namespace mmx
#endif //__MMX__POLYNOMIAL_NAIVE__HPP

---