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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : polynomial_dicho.hpp
* DESCRIPTION: dichotomic algorithms including Karatsuba multiplication
* 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_DICHO__HPP
#define __MMX__POLYNOMIAL_DICHO__HPP
#include <basix/vector_sort.hpp>
#include <algebramix/polynomial_naive.hpp>
#include <algebramix/crt_polynomial.hpp>

namespace mmx {
#define TMPL  template<typename C>
#define TMPLP template <typename Polynomial>
#define Vector vector<C>

/******************************************************************************
* Variants for dichotomic algorithms on polynomials
******************************************************************************/

template<typename V>
struct polynomial_karatsuba: public V {
  typedef typename V::Vec Vec;
  typedef typename V::Naive Naive;
  typedef polynomial_karatsuba<typename V::Positive> Positive;
  typedef polynomial_karatsuba<typename V::No_simd> No_simd;
  typedef polynomial_karatsuba<typename V::No_thread> No_thread;
  typedef polynomial_karatsuba<typename V::No_scaled> No_scaled;
};

template<typename F, typename V, typename W>
struct implementation<F,V,polynomial_karatsuba<W> >:
    public implementation<F,V,W> {};

template<typename V>
struct polynomial_dicho: public V {
  typedef typename V::Vec Vec;
  typedef typename V::Naive Naive;
  typedef polynomial_dicho<typename V::Positive> Positive;
  typedef polynomial_dicho<typename V::No_simd> No_simd;
  typedef polynomial_dicho<typename V::No_thread> No_thread;
  typedef polynomial_dicho<typename V::No_scaled> No_scaled;
};

template<typename F, typename V, typename W>
struct implementation<F,V,polynomial_dicho<W> >:
  public implementation<F,V,W> {};

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

template<typename V>
struct polynomial_multiply_threshold {};

template<typename V, typename W>
struct implementation<polynomial_multiply,V,polynomial_karatsuba<W> >:
  public implementation<polynomial_linear,V>
{
  typedef polynomial_multiply_threshold<polynomial_karatsuba<W> > Th;
  typedef implementation<vector_linear,V> Vec;
  typedef implementation<polynomial_linear,W> Pol;
  typedef implementation<polynomial_multiply,W> Fallback;

TMPL static void
multiply (C* dest, const C* s1, const C* s2, nat n1, nat n2) {
  if (n1 < Threshold(C,Th) || n2 < Threshold(C,Th))
    Fallback::mul (dest, s1, s2, n1, n2);
  else {
    nat p1= n1 >> 1, p2= n2 >> 1, P= p1+p2-1;
    nat spc= aligned_size<C,V> (3*(p1+p2));
    C* buf= mmx_new<C> (spc);
    C* low1= dest;
    C* low2= low1 + p1;
    C* mid1= buf;
    C* mid2= buf + p1;
    C* hi1 = mid2 + p2;
    C* hi2 = hi1 + p1;
    C* Low = mid1;
    C* Mid = hi1;
    C* Hi  = Mid + p1 + p2;
    Vec::half_copy (low1 , s1  , p1);
    Vec::half_copy (hi1  , s1+1, p1);
    Pol::add       (mid1 , low1, hi1 , p1);
    Vec::half_copy (low2 , s2  , p2);
    Vec::half_copy (hi2  , s2+1, p2);
    Pol::add       (mid2 , low2, hi2 , p2);
    multiply       (Hi   , hi1 , hi2 , p1, p2);
    multiply       (Mid  , mid1, mid2, p1, p2);
    multiply       (Low  , low1, low2, p1, p2);
    Pol::sub       (Mid  , Low , P);
    Pol::sub       (Mid  , Hi  , P);
    Pol::add       (Low+1, Hi  , P-1);
    Low[P]= Hi[P-1];
    Vec::double_copy (dest  , Low, P+1);
    Vec::double_copy (dest+1, Mid, P);
    mmx_delete<C> (buf, spc);

    if ((n1 & 1) != 0) {
      dest[(P<<1)+1]= C(0);
      Pol::mul_add (dest + (p1<<1), s2, s1[p1<<1], p2<<1);
    }
    if ((n2 & 1) != 0) {
      dest[n1+n2-2]= C(0);
      Pol::mul_add (dest + (p2<<1), s1, s2[p2<<1], n1);
    }
  }
}

TMPL static inline void
mul (C* dest, const C* s1, const C* s2, nat n1, nat n2) {
  if (n1 == 0 && n2 == 0) return;
  if (n1 == 0 || n2 == 0)
    Pol::clear (dest, n1 + n2 - 1);
  else
    multiply (dest, s1, s2, n1, n2);
}

TMPL static void
square (C* dest, const C* s, nat n) {
  if (n == 0) return;
  if (n < Threshold(C,Th))
    Fallback::square (dest, s, n);
  else {
    nat p= n >> 1, P= 2*p-1, spc= aligned_size<C,V> (6*p);
    C* buf= mmx_new<C> (spc);
    C* low= dest;
    C* mid= buf;
    C* hi = buf + 2*p;
    C* Low= buf;
    C* Mid= hi;
    C* Hi = Mid + 2*p;
    Vec::half_copy (low  , s  , p);
    Vec::half_copy (hi   , s+1, p);
    Pol::add       (mid  , low, hi, p);
    square         (Hi   , hi , p);
    square         (Mid  , mid, p);
    square         (Low  , low, p);
    Pol::sub       (Mid  , Low, P);
    Pol::sub       (Mid  , Hi , P);
    Pol::add       (Low+1, Hi , P-1);
    Low[P]= Hi[P-1];
    Vec::double_copy (dest  , Low, P+1);
    Vec::double_copy (dest+1, Mid, P);
    mmx_delete<C> (buf, spc);

    if ((n & 1) != 0) {
      dest[(P<<1)+1]= dest[2*n-2]= C(0);
      Pol::mul_add (dest + (p<<1), s, s[p<<1], p<<1);
      Pol::mul_add (dest + (p<<1), s, s[p<<1], n);
    }
  }
}

TMPL static void
tmultiply (C* dest, const C* s1, const C* s2, nat n1, nat n2) {
  if (n1 < Threshold(C,Th) || n2 < Threshold(C,Th))
    Fallback::tmul (dest, s1, s2, n1, n2);
  else {
    nat p1= n1 >> 1, p2= n2 >> 1, P= p1+p2-1;
    nat spc= aligned_size<C,V> (3*(p1+P)+1);
    C* buf= mmx_new<C> (spc);
    C* low1= buf;
    C* low2= low1 + p1;
    C* mid1= low2 + P+1;
    C* mid2= mid1 + p1;
    C* hi1 = mid2 + P;
    C* hi2 = hi1  + p1;
    C* Low = mid1;
    C* Mid = hi1;
    C* Hi  = dest;
    Vec::half_copy (low1, s1  , p1);
    Vec::half_copy (hi1 , s1+1, p1);
    Pol::add       (mid1, low1, hi1, p1);
    Vec::half_copy (low2, s2  , P+1);
    Vec::half_copy (hi2 , s2+1, P);
    Pol::add       (mid2, low2+1, hi2, P);
    tmultiply (Hi , hi1 , hi2 , p1, p2);
    tmultiply (Mid, mid1, mid2, p1, p2);
    tmultiply (Low, low1, low2, p1, p2+1);
    Pol::sub (Mid, Hi   , p2);
    Pol::sub (Mid, Low+1, p2);
    Pol::add (Low, Hi   , p2);
    Vec::double_copy (dest  , Low, p2);
    Vec::double_copy (dest+1, Mid, p2);
    mmx_delete<C> (buf, spc);

    if ((n1 & 1) != 0)
      Pol::mul_add (dest, s2+n1-1, s1[n1-1], n2);

    if ((n2 & 1) != 0)
      dest[n2-1]= Vec::inn_prod (s1, s2+n2-1, n1);
  }
}

TMPL static inline 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);
  if (n1 == 0 || n2 == 0) return;
  tmultiply (dest, s1, s2, n1, n2);
}

}; // implementation<polynomial_multiply,V,polynomial_karatsuba<W> >

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

template<typename V>
struct polynomial_divide_threshold {};
  
template<typename V, typename BV>
struct implementation<polynomial_divide,V,polynomial_dicho<BV> >:
  public implementation<polynomial_multiply,V>
{
  typedef polynomial_divide_threshold<polynomial_dicho<BV> > Th;
  typedef implementation<polynomial_multiply,V> Pol;
  typedef implementation<polynomial_divide,V,BV> Fallback;

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

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);
  }
}

TMPL static void
quo_rem (C* dest, C* s1, const C* s2, nat n1, nat n2) {
  // NOTE: requires C to be a field
  if (n1 < n2);
  else if (n1 < Threshold(C,Th))
    Fallback::quo_rem (dest, s1, s2, n1, n2);
  else {
    nat tot= aligned_size<C,V> ((n2 << 1) + n2);
    C* buf= mmx_new<C> (tot);
    C* inv= buf + (n2 << 1);
    nat nq = n1 + 1 - n2;
    nat l  = min (n2, nq);
    invert_hi (inv + n2 - l, s2 + n2 - l, l);
    while (n1 >= n2) {
      nat nq= n1 + 1 - n2;
      nat l = min (n2, nq);
      C* dest_hi= dest + nq - l;
      Pol::mul (buf, s1 + n1 - l, inv + n2 - l, l, l);
      Pol::copy (dest_hi, buf + l - 1, l);
      Pol::mul (buf, dest_hi, s2, l, n2);
      Pol::sub (s1 + n1 - (n2 + l - 1), buf, n2 + l - 1);
      n1 -= l;
    }
    mmx_delete<C> (buf, tot);
  }
}

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) occupies the one of the quotient.
  // We assume n2>0 and s2[n2-1] != 0.
  // (dest, n1) contains the output.
  // NOTE: requires C to be a field
  if (n1 < n2)
    Pol::copy (dest, s1, n1);
  else if (n1 < Threshold(C,Th))
    Fallback::tquo_rem (dest, s1, s2, n1, n2);
  else {
    nat tot= aligned_size<C,V> (3 * n2);
    C* buf= mmx_new<C> (tot);
    C*  inv= buf + 2*n2;
    nat nq = n1 + 1 - n2;
    nat l  = min (n2, nq);
    invert_hi (inv + n2 - l, s2 + n2 - l, l);
    Pol::neg (inv + n2 - l, inv + n2 - l, l);
    Pol::copy (dest, s1, n2-1);
    Pol::clear (dest+n2-1, n1-n2+1);
    nat m = n2-1;
    while (m < n1) {
      nat nq= n1 - m;
      nat l = min (n2, nq);
      Pol::clear (buf     , l-1);
      Pol::tmul  (buf+l-1 , s2          , dest + m - (n2 - 1), n2, l);
      Pol::sub   (buf+l-1 , s1 + m      , l);      
      Pol::tmul  (dest + m, inv + n2 - l, buf, l , l);
      m += l;
    }
    mmx_delete<C> (buf, tot);
  }
}

// Pseudo division. C can be any ring.
TMPL static void
pinvert_hi (C* dest, const C* src, nat n) {
  if (n == 1) *dest= C(1);
  else {
    nat h= (n+1) >> 1;
    nat l= n - h;
    nat tmp_size= aligned_size<C,V> (l);
    C* tmp= mmx_new<C> (tmp_size);
    pinvert_hi (dest + l, src + l, h);
    pinvert_hi (tmp, src + h, l); // FIXME: increase from l to 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_sc (dest + l, binpow (src[n-1], l), h);
    Pol::mul (buf, tmp, aux + h - 1, l, l);
      // FIXME: use truncated product
    Pol::neg (dest, buf + l - 1, l);
    mmx_delete<C> (buf, buf_size);
    mmx_delete<C> (tmp, tmp_size);
  }
}

TMPL static void
pquo_rem (C* dest, C* s1, const C* s2, nat n1, nat n2) {
  if (n1 < n2);
  else if (n1 < Threshold(C,Th))
    Fallback::pquo_rem (dest, s1, s2, n1, n2);
  else {
    nat tot= aligned_size<C,V> ((n2 << 1) + n2);
    C* buf= mmx_new<C> (tot);
    C* inv= buf + (n2 << 1);
    nat tmp_size= aligned_size<C,V> (n2);
    C* tmp= mmx_new<C> (tmp_size);
    nat nq= n1 + 1 - n2;
    nat l = min (n2, nq);
    nat l_end= nq % l;
    pinvert_hi (inv + n2 - l, s2 + n2 - l, l);
    if (l_end != 0)
      pinvert_hi (tmp + n2 - l_end, s2 + n2 - l_end, l_end);
    while (n1 >= n2) {
      nat nq= n1 + 1 - n2;
      nat l = min (n2, nq);
      C* dest_hi= dest + nq - l;
      if (l == l_end)
        Pol::mul (buf, s1 + n1 - l, tmp + n2 - l, l, l);
      else
        Pol::mul (buf, s1 + n1 - l, inv + n2 - l, l, l);
      Pol::copy (dest_hi, buf + l - 1, l);
      Pol::mul (buf, dest_hi, s2, l, n2);
      Pol::mul_sc (s1, binpow (s2[n2-1], l), n1);
      Pol::sub (s1 + n1 - (n2 + l - 1), buf, n2 + l - 1);
      Pol::mul_sc (dest_hi, binpow (s2[n2-1], nq - l), l);
      n1 -= l;
    }
    mmx_delete<C> (buf, tot);
    mmx_delete<C> (tmp, tmp_size);
  }
}

}; // implementation<polynomial_divide,V,polynomial_dicho<BV> >

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

// "Modern Computer Algebra", Algorithm 11.4

template<typename V>
struct polynomial_euclidean_threshold {};

template<typename V, typename BV>
struct implementation<polynomial_euclidean,V,polynomial_dicho<BV> >:
  public implementation<polynomial_divide,V>
{
  typedef polynomial_euclidean_threshold<polynomial_dicho<BV> > Th;
  typedef implementation<vector_linear,V> Vec;
  typedef implementation<polynomial_divide,V> Pol;
  typedef implementation<polynomial_euclidean,V,BV> Fallback;

private:

TMPL static void
dot_product (C* d, nat& nd,
             const C* r0, const C* r1, nat nr0, nat nr1,
             const C* s0, const C* s1, nat ns0, nat ns1,
             C* t) {
  // d = r0 s0 + r1 s1
  if ((nr0 == 0 || ns0 == 0) && (nr1 == 0 || ns1 == 0)) {
    nd = 0; return;
  }
  if (nr0 == 0 || ns0 == 0) {
    nd = nr1 + ns1 - 1;
    Pol::mul (d, r1, s1, nr1, ns1);
    Pol::trim (d, nd);
    return;
  }
  if (nr1 == 0 || ns1 == 0) {
    nd = nr0 + ns0 - 1;
    Pol::mul (d, r0, s0, nr0, ns0);
    Pol::trim (d, nd);
  }
  if (nr0 + ns0 - 1 < nr1 + ns1 - 1) {
    nd = nr1 + ns1 - 1;
    Pol::mul (d, r1, s1, nr1, ns1);
    Pol::mul (t, r0, s0, nr0, ns0);
    Pol::add (d, t , nr0 + ns0 - 1);
  }
  else {
    nd = nr0 + ns0 - 1;
    Pol::mul (d, r0, s0, nr0, ns0);
    Pol::mul (t, r1, s1, nr1, ns1);
    Pol::add (d, t , nr1 + ns1 - 1);
  }
  Pol::trim (d, nd);
}
  
TMPL static void
matrix_vector_product (C* r0, C* r1, nat& nr0, nat& nr1,
                       const C* R00, const C* R01, const C* R10, const C* R11,
                       nat nR00, nat nR01, nat nR10, nat nR11,
                       const C* s0, const C* s1, nat n0, nat n1,
                       C* tp) {
  // (r0,r1) = R . (s0,s1)
  dot_product (r0, nr0, R00, R01, nR00, nR01, s0, s1, n0, n1, tp);
  dot_product (r1, nr1, R10, R11, nR10, nR11, s0, s1, n0, n1, tp);
}

TMPL static void
matrix_product (C* Q00, C* Q01, C* Q10, C* Q11,
                nat& nQ00, nat& nQ01, nat& nQ10, nat& nQ11,
                const C* S00, const C* S01, const C* S10, const C* S11,
                nat nS00, nat nS01, nat nS10, nat nS11,
                const C* R00, const C* R01, const C* R10, const C* R11,
                nat nR00, nat nR01, nat nR10, nat nR11, C* tp) {
  // Q = S . R
  matrix_vector_product (Q00, Q10, nQ00, nQ10,
                         S00, S01, S10, S11, nS00, nS01, nS10, nS11,
                         R00, R10, nR00, nR10, tp);
  matrix_vector_product (Q01, Q11, nQ01, nQ11, 
                         S00, S01, S10, S11, nS00, nS01, nS10, nS11,
                         R01, R11, nR01, nR11, tp);
}

TMPL static void
new_matrix (C*& M00, C*& M01, C*& M10, C*& M11, nat l) {
  M00= mmx_new<C> (l); M01= mmx_new<C> (l);
  M10= mmx_new<C> (l); M11= mmx_new<C> (l);
}

TMPL static void
delete_matrix (C* M00, C* M01, C* M10, C* M11, nat l) {
  mmx_delete<C> (M00, l); mmx_delete<C> (M01, l);
  mmx_delete<C> (M10, l); mmx_delete<C> (M11, l);
}

TMPL static void
half_gcd (C* Q00, C* Q01, C* Q10, C* Q11,
          nat& nQ00, nat& nQ01, nat& nQ10, nat& nQ11,
          const C* r0, const C* r1, nat n0, nat n1, nat k,
          C* rho, C* tp) {
  // Performs the maximum numbers of reductions less than k.
  // tp must have size n0 + n1.
  // rho is the vector of the leading coefficients of to
  // degree min (n0, n1) - 1.
  // s1 and s2 are supposed to be monic.
  VERIFY (n0 >= n1, "bad input sizes");
  VERIFY (k  <= n0, "index k out of range");
  if (n1 == 0 || k < n0 - n1 + 1) {
    Q00[0] = 1; nQ00 = 1; nQ01 = 0;
    Q11[0] = 1; nQ10 = 0; nQ11 = 1;
    return;
  }
  if (k == 1) {
    Q00[0] = 1; nQ00 = 1; nQ01 = 0;
    Q10[0] = 1; Q11[0] = - r0[n0-1] / r1[n1-1]; nQ10 = 1; nQ11 = 1;
    return;
  }
  nat h = k >> 1, h2 = (h << 1) - 1;
  nat len_R= aligned_size<C,V> (h2);
  nat nR00, nR01, nR10, nR11;
  C* R00, * R01, * R10, * R11;
  new_matrix (R00, R01, R10, R11, len_R);
  if (h2 < n0 - n1)
    half_gcd (R00, R01, R10, R11, nR00, nR01, nR10, nR11,
              r0, r1, n0, n1, h, rho, tp);
  else
    half_gcd (R00, R01, R10, R11, nR00, nR01, nR10, nR11,
              r0 + n0 - h2, r1 + n0 - h2,
              h2, h2 - (n0 - n1), h,
              rho == NULL ? rho : (rho + (n0 - h2)), tp);    
  nat len_r= aligned_size<C,V> (n0 + h);
  C* rjm1= mmx_new<C> (len_r), * rj  = mmx_new<C> (len_r);
  nat nj, njm1;
  // TODO << use truncated product
  matrix_vector_product (rjm1, rj, njm1, nj,
                         R00, R01, R10, R11, nR00, nR01, nR10, nR11,
                         r0, r1, n0, n1, tp);
  if (nj == 0 || k < n0 - nj + 1) {
    Pol::copy (Q00, R00, nR00); nQ00= nR00;
    Pol::copy (Q01, R01, nR01); nQ01= nR01;
    Pol::copy (Q10, R10, nR10); nQ10= nR10;
    Pol::copy (Q11, R11, nR11); nQ11= nR11;
    return;
  }
  nat len_S= aligned_size<C,V> (max (k - (n0 - nj), njm1 - nj + 1));
  nat nS00, nS01, nS10, nS11;
  C* S00, * S01, * S10, * S11;
  new_matrix (S00, S01, S10, S11, len_S);

  nat len_T= aligned_size<C,V> (njm1 - nj + h);
  nat nT00, nT01, nT10, nT11;
  C* T00, * T01, * T10, * T11;
  new_matrix (T00, T01, T10, T11, len_T);

  S01[0]= 1; S10[0]= 1;  
  nS00= 0; nS01= 1; nS10= 1; nS11= njm1 - nj + 1;
  quo_rem (S11, rjm1, rj, njm1, nj);
  Vec::neg (S11, nS11);
  Pol::trim (rj, nj); Pol::trim (rjm1, njm1);

  if (njm1 != 0) {
    C c= 1 / rjm1[njm1-1];
    if (rho != NULL) rho[njm1-1] = rjm1[njm1-1];
    Pol::mul_sc (rjm1, c, njm1);
    Pol::mul_sc (S11, c, nS11);
    Pol::mul_sc (S10, c, nS10);
  }
  matrix_product (T00, T01, T10, T11, nT00, nT01, nT10, nT11,
                  S00, S01, S10, S11, nS00, nS01, nS10, nS11,
                  R00, R01, R10, R11, nR00, nR01, nR10, nR11, tp);

  if (njm1 == 0 || k < n0 - njm1 + 1) {
    Pol::copy (Q00, T00, nT00); nQ00= nT00;
    Pol::copy (Q01, T01, nT01); nQ01= nT01;
    Pol::copy (Q10, T10, nT10); nQ10= nT10;
    Pol::copy (Q11, T11, nT11); nQ11= nT11;
    return;
  }
  h = k - (n0 - nj); h2 = (h << 1) - 1;
  if (h2 < nj - njm1 || h2 > nj)
    half_gcd (S00, S01, S10, S11, nS00, nS01, nS10, nS11,
              rj, rjm1, nj, njm1, h, rho, tp);
  else
    half_gcd (S00, S01, S10, S11, nS00, nS01, nS10, nS11,
              rj + nj - h2, rjm1 + nj - h2, h2, h2 - (nj - njm1), h,
              rho == NULL ? rho : (rho + (nj - h2)), tp);
  mmx_delete<C> (rjm1, len_r); mmx_delete<C> (rj, len_r);
  matrix_product (Q00, Q01, Q10, Q11, nQ00, nQ01, nQ10, nQ11,
                  S00, S01, S10, S11, nS00, nS01, nS10, nS11,
                  T00, T01, T10, T11, nT00, nT01, nT10, nT11, tp);
  delete_matrix (R00, R01, R10, R11, len_R);
  delete_matrix (S00, S01, S10, S11, len_S);
  delete_matrix (T00, T01, T10, T11, len_T);
}
  
public:

TMPL static inline 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) {
  Fallback::euclidean_sequence (s1, s2, n1, n2, d1, d2, m1, m2,
                                u1, u2, nu1, nu2, v1, v2, nv1, nv2,
                                n, rho, q, r, co1, co2, k);
}

TMPL static void
gcd (C* g, nat& n, const C* s1, const C* s2, nat n1, nat n2,
     C* uu1, C* uu2, nat& nuu1, nat& nuu2) {
  // g must have allocated length min (n1, n2)
  if (n1 < Threshold(C,Th) || n2 < Threshold(C,Th)) {
    Fallback::gcd (g, n, s1, s2, n1, n2, uu1, uu2, nuu1, nuu2); return; }
  VERIFY (n1>0 && n2>0 && s1[n1-1] != 0 && s2[n2-1] != 0,
          "invalid hypothesis for gcd computation");
  nat nu1, nu2, nv1, nv2;
  C c1= 1 / s1[n1-1], c2= 1 / s2[n2-1];
  nat l1= aligned_size<C,V> (n1), l2= aligned_size<C,V> (n2);
  C* z1= mmx_new<C> (l1), * z2= mmx_new<C> (l2);
  Pol::mul_sc (z1, s1, c1, n1); Pol::mul_sc (z2, s2, c2, n2); 
  C* u1= mmx_new<C> (l2), * u2= mmx_new<C> (l1);
  C* v1= mmx_new<C> (l2), * v2= mmx_new<C> (l1);
  nat len_tp= aligned_size<C,V> (n1 + n2);
  C* tp1= mmx_new<C> (len_tp), * tp2= mmx_new<C> (len_tp);
  C* rho= NULL;
  if (n1 >= n2)
    half_gcd (u1, u2, v1, v2, nu1, nu2, nv1, nv2,
              z1, z2, n1, n2, n1, rho, tp1);
  else
    half_gcd (u2, u1, v2, v1, nu2, nu1, nv2, nv1,
              z2, z1, n2, n1, n2, rho, tp1);
  VERIFY (nu1 <= n2 && nv1 <= n2 && nu2 <= n1 && nv2 <= n1, "bug");
  dot_product (tp2, n, u1, u2, nu1, nu2, z1, z2, n1, n2, tp1);
  VERIFY (n <= min (n1, n2), "bug");
  C c= n == 0 ? 0 : (1 / tp2[n-1]);
  Pol::mul_sc (g, tp2, c, n); Pol::clear (g + n, min (n1, n2) - n);
  if (uu1 != NULL) {
    Pol::mul_sc (uu1, u1, c1 * c, nu1);
    Pol::clear (uu1 + nu1, n2 - nu1); nuu1= nu1; }
  if (uu2 != NULL) {
    Pol::mul_sc (uu2, u2, c2 * c, nu2);
    Pol::clear (uu2 + nu2, n1 - nu2); nuu2= nu2; }
  mmx_delete<C> (tp1, len_tp); mmx_delete<C> (tp2, len_tp);
  mmx_delete<C> (u1, l2); mmx_delete<C> (u2, l1);
  mmx_delete<C> (v1, l2); mmx_delete<C> (v2, l1);
}  

TMPL static void
gcd (C* g, nat& n, const C* s1, const C* s2, nat n1, nat n2) {
  nat nuu1, nuu2;
  C* uu1= NULL, * uu2= NULL;
  gcd (g, n, s1, s2, n1, n2, uu1, uu2, nuu1, nuu2);
}

TMPL static void
gcd (C* g, nat& n, const C* s1, const C* s2, nat n1, nat n2,
     C* uu1, nat& nuu1) {
  C* uu2= NULL; nat nuu2;
  gcd (g, n, s1, s2, n1, n2, uu1, uu2, nuu1, nuu2);
}  

TMPL static void
pade (C* r, C* t, const C* s, nat m, nat n, nat k) {
  // r (resp. t) must have allocated length at least k (resp. n-k+1)
  // s = r / t + O(x^n), deg r < k, deg t <= n-k
  // "Modern Computer Algebra", Corollary 5.21
  if (n < Threshold(C,Th) || k == n || k == 0) {
    Fallback::pade (r, t, s, m, n, k); return; }
  VERIFY (n > k && n > 0 && m > 0 && m <= n && s[m-1] != 0,
          "invalid hypothesis for gcd computation");
  nat nu1, nu2, nv1, nv2;
  nat n1= n+1, n2= m;
  nat l1= aligned_size<C,V> (n1), l2= aligned_size<C,V> (n2);
  C* u1= mmx_new<C> (l2), * u2= mmx_new<C> (l1);
  C* v1= mmx_new<C> (l2), * v2= mmx_new<C> (l1);
  C* s1= mmx_new<C> (l1); Pol::clear (s1, n); s1[n]= 1; const C* s2= s;
  nat len_tp= aligned_size<C,V> (n1 + n2);
  C* tp= mmx_new<C> (len_tp);
  C* rho= NULL;
  half_gcd (u1, u2, v1, v2, nu1, nu2, nv1, nv2,
            s1, s2, n1, n2, n-k, rho, tp);

  nat len_r= aligned_size<C,V> (2*n1-k);
  C* rjm1= mmx_new<C> (len_r), * rj  = mmx_new<C> (len_r);
  nat nj, njm1;
  matrix_vector_product (rjm1, rj, njm1, nj,
                         u1, u2, v1, v2, nu1, nu2, nv1, nv2,
                         s1, s2, n1, n2, tp);
  if (nj <= k) {
    VERIFY (nv2 <= n-k+1, "bug");
    Pol::copy (r, rj, nj); Pol::clear (r + nj, k - nj);
    Pol::copy (t, v2, nv2); Pol::clear (t + nv2, n-k+1 - nv2);
  }
  else {
    VERIFY (nj > k, "bug");
    nat nq= njm1 - nj + 1, len_q= aligned_size<C,V> (nq);
    C* q= mmx_new<C> (len_q);
    quo_rem (q, rjm1, rj, njm1, nj);
    Pol::trim (rjm1, njm1);
    Pol::mul (tp, q, v2, nq, nv2);
    Pol::clear (u2 + nu2, n1 - nu2); nu2= n1;
    Pol::sub (u2, tp, nq + nv2 - 1);
    Pol::trim (u2, nu2);
    VERIFY (njm1 <= k, "bug");
    VERIFY (nu2 <= n-k+1, "bug");
    Pol::copy (r, rjm1, njm1); Pol::clear (r + njm1, k - njm1);
    Pol::copy (t, u2, nu2); Pol::clear (t + nu2, n-k+1 - nu2);
    mmx_delete<C> (q, len_q);
  }
  mmx_delete<C> (rjm1, len_r); mmx_delete<C> (rj, len_r);
  mmx_delete<C> (s1, l1); mmx_delete<C> (tp, len_tp);
  mmx_delete<C> (u1, l2); mmx_delete<C> (u2, l1);
  mmx_delete<C> (v1, l2); mmx_delete<C> (v2, l1);
}

}; // implementation<polynomial_euclidean,V,polynomial_dicho<BV> >

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

template<typename V>
struct polynomial_evaluate_threshold {};

template<typename V, typename BV>
struct implementation<polynomial_evaluate,V,polynomial_dicho<BV> >:
  public implementation<polynomial_divide,V>
{
  typedef polynomial_evaluate_threshold<polynomial_dicho<BV> > Th;
  typedef implementation<vector_linear,V> Vec;
  typedef implementation<polynomial_divide,V> Pol;
  typedef implementation<polynomial_evaluate,V,BV> Fallback;

TMPL static inline void
factorials (C* dest, nat n) {
  if (n > 0) dest[0]= C(1);
  for (nat i=1; i<n; i++)
    dest[i]= C(i) * dest[i-1];
}

TMPL static void
shift (C* dest, const C* s, const C& sh, nat n) {
  if (n <= 1 || sh == 0) Pol::copy (dest, s, n);
  else {
    nat l= aligned_size<C,V> (5 * n);
    C* u    = mmx_new<C> (l);
    C* v    = u + n;
    C* w    = v + n;
    C* facts= w + n + n;
    factorials (facts, n);
    Vec::mul (u, s, facts, n);
    Vec::vec_reverse (u, n);
    Vec::set (v, 1, n);
    if (sh != 0) Pol::q_difference (v, v, sh, n);
    Vec::div (v, facts, n);
    Pol::mul (w, u, v, n, n); // FIXME: rather use truncated multiplication
    Vec::vec_reverse (w, n);
    Vec::div (dest, w, facts, n);
    mmx_delete<C> (u, l);
  }
}

TMPL static inline C
evaluate (const C* p, const C& x, nat l) {
  return Fallback::evaluate (p, x, l);
}

TMPL static void
q_binomial (C* dest, const C& q, nat mu) {
  dest[mu]= 1;
  for (nat i=1; i<=mu; i++)
    dest[mu-i]= (C (mu+1-i) * dest[mu+1-i] * q) / C(i);
}

TMPL static void
expand (C** v, const C* p, const C* x, const nat* nu, nat n, nat k) {
  nat tot= 0;
  for (nat i=0; i<k; i++)
    tot += nu[i] + 1;
  nat* d= mmx_new<nat> (k);
  nat l= aligned_size<C,V> (tot << 1);
  C* q= mmx_new<C> (l);
  C* r= q + tot;
  nat off= 0;
  for (nat i=0; i<k; i++) {
    q_binomial (q + off, -x[i], nu[i]);
    d[i]= nu[i] + 1;
    off += d[i];
  }
  multi_mod (r, p, q, d, n, k);
  off= 0;
  for (nat i=0; i<k; i++) {
    shift (v[i], r + off, x[i], nu[i]);
    off += d[i];
  }
  mmx_delete<C> (q, l);
  mmx_delete<nat> (d, k);
}

#define C typename scalar_type_helper<Polynomial >::val

struct _vector_sort_by_increasing_degree_op {
  TMPLP static bool
  op (const Polynomial& p, const Polynomial& q) {
    return deg (p) < deg (q); }
  TMPLP static bool
  not_op (const Polynomial& p, const Polynomial& q) {
    return deg (p) >= deg (q); }
};

template<typename Op, typename Polynomial> static inline vector<Polynomial>
_multi_rem (const Polynomial& p, const vector<Polynomial>& q) {
  if (p == 0) return vector<Polynomial> (Polynomial(C(0)), N(q));
  nat n= degree (p);
  vector<Polynomial> sorted_q (q), r (Polynomial(C(0)), N(q));
  vector<nat> sigma; // permutation
  sort_leq<_vector_sort_by_increasing_degree_op> (sorted_q, sigma);
  nat start= 0, sum= 0;
  for (nat i= 0; i < N(q); i++) {
    sum += degree (sorted_q[i]);
    if (sum > n / 2 || i+1 == N(q)) {
      Crt_polynomial_transformer(Polynomial)
        crter (range (sorted_q, start, i+1));
      vector<Polynomial> tmp; direct_crt (tmp, p, crter);
      for (nat j= start; j < i+1; j++) r[sigma[j]]= tmp[j-start];
      start= i+1; sum= 0;
    }
  }
  return r; }

TMPLP static inline vector<Polynomial>
multi_rem (const Polynomial& p, const vector<Polynomial>& q) {
  return _multi_rem<rem_op> (p, q); }

TMPLP static inline vector<Polynomial>
multi_prem (const Polynomial& p, const vector<Polynomial>& q) {
  return _multi_rem<prem_op> (p, q); }

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 xector expected");
  if (N(x) == 0) return 1;
  vector<Polynomial> q (Polynomial (), N(x));
  Polynomial z (C(1), 1);
  for (nat i= 0; i < N(x); i++) q[i]= z - x[i];
  Crt_polynomial_transformer(Polynomial) crter (q);
  return * moduli_product (crter); }

TMPLP static inline Vector
evaluate (const Polynomial& p, const Vector& x) {
  ASSERT (is_non_scalar (x), "non-scalar vector expected");
  if (N(x) == 0) return Vector (C(0), 0);
  vector<Polynomial> q (Polynomial (), N(x));
  Polynomial z (C(1), 1);
  for (nat i= 0; i < N(x); i++) q[i]= z - x[i];
  vector<Polynomial> tmp= multi_rem (p, q);
  Vector r (C(0), N(x));
  for (nat i= 0; i < N(x); i++) r[i]= tmp[i][0];
  return r;
}

TMPLP static inline Polynomial
tevaluate (const Vector& v, const Vector& x, nat l) {
  ASSERT (is_non_scalar (x), "non-scalar vector expected");
  if (l == 0) return Polynomial (0);
  nat n= N(x), ll= aligned_size<C,V> (l);
  vector<Polynomial> q (Polynomial (), n);
  Polynomial z (C(1), 1);
  for (nat i= 0; i < N(v); i++) q[i]= Polynomial (1) - x[i] * z;
  Crt_polynomial_transformer(Polynomial) crter (q);
  Polynomial num (combine_crt (v, crter)), den (* moduli_product (crter));
  // retrieve the coefficients of the series num / den
  C* tmp= mmx_new<C> (ll), * inv= mmx_new<C> (ll);
  Pol::clear (tmp, l); Pol::copy (tmp, seg(den), min (N(den), l));
  Pol::invert_lo (inv, tmp, l);
  Polynomial b (inv, l, ll), c (num * b);
  for (nat i= 0; i < l; i++) tmp[i]= c[i];
  return Polynomial (tmp, l, ll);
}

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");
  if (N(x) == 0) return Polynomial (0);
  vector<Polynomial> q (Polynomial (), N(x));
  Polynomial z (C(1), 1);
  for (nat i= 0; i < N(x); i++) q[i]= z - x[i];
  Crt_polynomial_transformer(Polynomial) crter (q);
  Polynomial ans; inverse_crt (ans, as<vector<Polynomial> > (v), crter);
  return ans;
}

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);
  if (n == 0) return Vector (C(0), 0);
  vector<Polynomial> q (Polynomial (), n);
  Polynomial z (C(1), 1);
  for (nat i= 0; i < n; i++) q[i]= z - x[i];
  Crt_polynomial_transformer(Polynomial) crter (q);
  Polynomial rp (lshiftz (reverse (p), (int) (n - N(p))));
  Polynomial den (* moduli_product (crter));
  vector<Polynomial> a; direct_crt (a, range (den * rp, n, 2*n), crter);
  vector<Polynomial> b; direct_crt (b, derive (den), crter);
  Vector v (C(0), n); for (nat i= 0; i < n; i++) v[i]= a[i][0] / b[i][0];
  return v;
}

#undef C
}; // implementation<polynomial_evaluate,V,polynomial_dicho<BV> >

#undef TMPL
#undef TMPLP
#undef Vector
} // namespace mmx
#endif //__MMX__POLYNOMIAL_DICHO__HPP