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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : polynomial_ring_dicho.hpp
* DESCRIPTION: dichotomic algorithms over rings
* COPYRIGHT  : (C) 2009  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_RING_DICHO__HPP
#define __MMX__POLYNOMIAL_RING_DICHO__HPP
#include <algebramix/polynomial_ring_naive.hpp>
#include <algebramix/polynomial_dicho.hpp>
#include <algebramix/matrix.hpp>

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

/******************************************************************************
* Dichotomic variants for polynomials over rings without division
******************************************************************************/

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

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

DEFINE_VARIANT_1 (typename V, V,
                  polynomial_ring_dicho,
                  polynomial_ring_dicho_inc<
                    polynomial_ring_naive<
                      polynomial_dicho<V> > >)

/******************************************************************************
* Dichotomic variants for polynomials over rings with gcd without division
******************************************************************************/

template<typename V>
struct polynomial_gcd_ring_dicho_inc:
  public V {};

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

DEFINE_VARIANT_1 (typename V, V,
                  polynomial_gcd_ring_dicho,
                  polynomial_gcd_ring_dicho_inc<
                    polynomial_gcd_ring_naive<
                  polynomial_ring_dicho<V> > >)

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

template<typename V, typename W>
struct implementation<polynomial_divide,V,polynomial_ring_dicho_inc<W> >:
  public implementation<polynomial_multiply,V>
{
  typedef polynomial_divide_threshold<polynomial_ring_dicho<W> > Th;
  typedef implementation<polynomial_multiply,V> Pol;
  typedef implementation<polynomial_divide,polynomial_naive> Fallback;

TMPL static void
quo_rem (C* dest, C* s1, const C* s2, nat n1, nat n2) {
  // NOTE: C must support exact division
  if (n1 < n2);
  else if (n1 < Threshold(C,Th))
    Fallback::quo_rem (dest, s1, s2, n1, n2);
  else {
    nat l= n1-n2+1;
    if (n2 > l) {
      nat d= n2-l;
      quo_rem (dest, s1+d, s2+d, n1-d, l);
      nat buf_size= aligned_size<C,V> (l+d-1);
      C* buf= mmx_new<C> (buf_size);
      Pol::mul (buf, dest, s2, l, d);
      Pol::sub (s1, buf, l+d-1);
      mmx_delete<C> (buf, buf_size);      
    }
    else {
      nat h= (n2+1) >> 1;
      nat d= l-h;
      quo_rem (dest+d, s1+d, s2, n1-d, n2);
      quo_rem (dest, s1, s2, n1-h, n2);
    }
  }
}

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: C can be an euclidean ring (not necessarily 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 l= n1-n2+1;
    if (n2 > l) {
      nat d= n2-l;
      nat buf_size= aligned_size<C,V> (l+d-1);
      C* buf= mmx_new<C> (buf_size);
      nat tmp_size= aligned_size<C,V> (n1);
      C* tmp= mmx_new<C> (tmp_size);
      Pol::neg      (buf     , s1     , l+d-1);
      Pol::tmul     (tmp+n2-1, s2     , buf, d, l);
      Pol::add      (tmp+n2-1, s1+n2-1, l);
      Pol::copy     (tmp     , s1     , n2-1);
      Pol::tquo_rem (dest+d  , tmp+d  , s2+d, n1-d, l);
      Pol::copy     (dest    , s1     , d);
      mmx_delete<C> (buf, buf_size);
      mmx_delete<C> (tmp, tmp_size);
    }
    else {
      nat h= (n2+1) >> 1;
      nat d= l-h;
      nat tmp_size= aligned_size<C,V> (n1-d);
      C* tmp= mmx_new<C> (tmp_size);
      tquo_rem (dest    , s1      , s2, n1-h, n2);
      Pol::copy  (tmp     , dest + d, n2-1);
      Pol::copy  (tmp+n2-1, s1+n1-h , h);
      tquo_rem (dest+d  , tmp     , s2, n1-d, n2);
      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 l= n1-n2+1;
    if (n2 > l) {
      nat d= n2-l;
      pquo_rem (dest, s1+d, s2+d, n1-d, l);
      nat buf_size= aligned_size<C,V> (l+d-1);
      C* buf= mmx_new<C> (buf_size);
      Pol::mul (buf, dest, s2, l, d);
      Pol::mul_sc (s1, binpow (s2[n2-1], l), d);
      Pol::sub (s1, buf, l+d-1);
      mmx_delete<C> (buf, buf_size);      
    }
    else {
      nat h= (n2+1) >> 1;
      nat d= l-h;
      pquo_rem (dest+d, s1+d, s2, n1-d, n2);
      Pol::mul_sc (s1, binpow (s2[n2-1], h), d);
      pquo_rem (dest, s1, s2, n1-h, n2);
      Pol::mul_sc (dest+d, binpow (s2[n2-1], d), h);
    }
  }
}

}; // implementation<polynomial_divide,V,polynomial_ring_dicho_inc<W> >

/******************************************************************************
* Dichotomic subresultant
******************************************************************************/

template<typename V, typename BV>
struct implementation<polynomial_subresultant_base,V,
                      polynomial_gcd_ring_dicho_inc<BV> > {

template<typename T> static inline T
defected_pow (const T& x, const T& y, int n) {
  // return x^n / y^(n-1)
  VERIFY (n >= 0, "bug");
  if (n == 0) return y;
  if (n == 1) return x;
  T z= square_op::op (defected_pow (x, y, n >> 1)) / y;
  return ((n & 1) == 0) ? z : ((x * z) / y);
}

#define C typename scalar_type_helper<Polynomial>::val

TMPLP static void
subresultant_compose (matrix<Polynomial>& num_R,
                      const matrix<Polynomial>& num_S,
                      const C& den_S, nat k= 0) {
  if (k == 0)
    num_R= num_S * num_R;
  else {
    num_R(0,0)= range (range (num_R(0,0), 0, k) * range (num_S(0,0), 0, k) +
                       range (num_R(0,1), 0, k) * range (num_S(1,0), 0, k),
                       0, k);
    num_R(0,1)= range (range (num_R(0,0), 0, k) * range (num_S(0,1), 0, k) +
                       range (num_R(0,1), 0, k) * range (num_S(1,1), 0, k),
                       0, k);
    num_R(1,0)= range (range (num_R(1,0), 0, k) * range (num_S(0,0), 0, k) +
                       range (num_R(1,1), 0, k) * range (num_S(1,0), 0, k),
                       0, k);
    num_R(1,1)= range (range (num_R(1,0), 0, k) * range (num_S(0,1), 0, k) +
                       range (num_R(1,1), 0, k) * range (num_S(1,1), 0, k),
                       0, k);
  }
  num_R(0,0) /= den_S; num_R(0,1) /= den_S; 
  num_R(1,0) /= den_S; num_R(1,1) /= den_S; }

TMPLP static void
subresultant_compose (matrix<Polynomial>& num_R,
                      const vector<matrix<Polynomial> >& num_S,
                      const vector<C>& den_S, nat k=0) {
  // dichotomic computation of (big_mul (num_S) * num_R) / den_S
  // truncated to precision k (if k > 0).
  nat n= N(num_S);
  if (n == 0) return;
  if (n == 1) {
    subresultant_compose (num_R, num_S[0], den_S[0], k);
    return;
  }
  nat m= n >> 1;
  matrix<Polynomial> num_T (Polynomial(C(0)), 2, 2);
  C det_num_R= num_R(0,0)[0] * num_R(1,1)[0] - num_R(1,0)[0] * num_R(0,1)[0];
  num_T(0,0)= num_T(1,1)= det_num_R;
  subresultant_compose (num_T, range (num_S, 0, m), range (den_S, 0, m));
  subresultant_compose (num_R, num_T, det_num_R);
  det_num_R= num_R(0,0)[0] * num_R(1,1)[0] - num_R(1,0)[0] * num_R(0,1)[0];
  num_T(0,0)= num_T(1,1)= det_num_R;
  num_T(0,1)= num_T(1,0)= 0;
  subresultant_compose (num_T, range (num_S, m, n), range (den_S, m, n));
  subresultant_compose (num_R, num_T, det_num_R); }

TMPLP static void
half_subresultant_rec (const Polynomial& r0, const Polynomial& r1,
                       bool defected, const C& num,
                       matrix<Polynomial>& num_R, C& den_R, nat& k_out,
                       vector<matrix<Polynomial> >& Q, vector<C>& rho,
                       nat k) {
  // Get a maximum total of at most k terms of the quotients.
  // k = n0 performs a complete sequence.
  // defected means that r0 and r1 are strict truncations.
  nat n0= N(r0), n1= N(r1), l, l0, l1, h1, k2= (k << 1) - 1;
  VERIFY (n0 > n1, "bad input sizes");
  VERIFY (k <= n0, "index k out of range");
  if (k < n0 - n1 + 1) {
    num_R= matrix<Polynomial> (Polynomial (C(0)), 2, 2);
    num_R(0,0)= num; num_R(1,1)= num; den_R= num;
    k_out= 1;
    return;
  }
  half_subresultant_rec (r0, r1, defected, num,
                         num_R, den_R, k_out, Q, rho, k >> 1);
  VERIFY (k_out <= (k >> 1), "bug");

  Polynomial s0, s1, t0, t1, q;
  if (k2 >= n0) {
    t0= r0; t1= r1;
  }
  else {
    t0= range (r0, n0 - k2, n0); t1= range (r1, n0 - k2, n1);
    defected= true;
  }
  s1= num_R(1,0) * t0 + num_R(1,1) * t1;
  if (s1 == 0) { return; }
  s0= num_R(0,0) * t0 + num_R(0,1) * t1;
  nat m0= N(s0), m1= N(s1);
  if (defected) {
    s1= range (s1, k_out - 1, m1);
    if (s1 == 0) { return; }
    s0= range (s0, k_out - 1, m0);
    m0 -= k_out - 1; m1 -= k_out - 1;
  }
  s0 /= den_R; s1 /= den_R;
  if (k < k_out + m0 - m1) return;
  C c0= s0[m0-1], c1= s1[m1-1], den_S= square_op::op (c0);
  C a= defected_pow (c1, c0, m0 - m1 - 1);
  t0= (m0 == m1 + 1) ? s1 : ((a * s1) / c0);
  C b0= t0[m1-1];
  t1= rem (b0 * c1 * s0, s1, q);
  l0= N(t0), l1= N(t1), l= N(q);
  matrix<Polynomial> num_S (Polynomial (C(0)), 2, 2);
                       num_S(0,0)= 0         ; num_S(0,1)= a * c0;
  if ((m0-m1-1) & 1) { num_S(1,0)= (-b0) * c1; num_S(1,1)=  q; }
  else               { num_S(1,0)=   b0  * c1; num_S(1,1)= -q; }
  Q << num_S; rho << den_S;
  subresultant_compose (num_R, num_S, den_S); k_out += l - 1;
  if (defected && l0 >= l-1 && l1 >= l-1) {
    t0= range (t0, l - 1, l0); t1= range (t1, l - 1, l1);
    l0 -= l - 1; l1 -= l - 1;
  }
  t1 /= ((m0-m1-1) & 1) ? (-den_S) : den_S;
  C det_num_R= num_R(0,0)[0] * num_R(1,1)[0] - num_R(1,0)[0] * num_R(0,1)[0];
  if (t1 == 0) { return; }
  half_subresultant_rec (t0, t1, defected, det_num_R,
                         num_S, den_S, h1, Q, rho, k - k_out + 1);
  VERIFY (h1 <= k - k_out + 1, "bug");
  subresultant_compose (num_R, num_S, den_S); k_out += h1 - 1; }

TMPLP static void
half_subresultant (const Polynomial& f, const Polynomial& g,
                   matrix<Polynomial>& R,
                   vector<matrix<Polynomial> >& Q, vector<C>& rho, nat k) {
  // Handle the two first (particular) rounds before entering the recurrence.
  nat n= N(f), m= N(g), k_out, k_rem;
  VERIFY (m > 1 && n >= m, "bad input sizes");
  VERIFY (k <= n, "index k out of range");
  C d= binpow (g[m-1], n-m);
  Polynomial q;
  Polynomial s0= g, s1= prem (f, g, q); k_rem = k - N(q) + 1;
  nat m0= N(s0), m1= N(s1); 
  R= matrix<Polynomial> (Polynomial (C(0)), 2, 2);
  R(0,0)= 0         ; R(0,1)= 1;
  R(1,0)= g[m-1] * d; R(1,1)= - q;
  Q= vec (R); rho= vec (C(1)); 
  if (m1 == 0 || k < n - m1 + 1) return;
  C c0= s0[m0-1], c1= s1[m1-1];
  C a= defected_pow (c1, d, m0-m1-1);
  Polynomial t0= (m0 == m1+1) ? s1 : ((a * s1) / d);
  C b0= t0[m1-1];
  Polynomial t1= rem (b0 * c1 * s0, s1, q); k_rem -= N(q) - 1;
  matrix<Polynomial> num_S (Polynomial (C(0)), 2, 2); C den_S (d * c0);
  t1 /= ((m0-m1-1) & 1) ? -den_S : den_S;
                       num_S(0,0)= 0         ; num_S(0,1)= a * c0;
  if ((m0-m1-1) & 1) { num_S(1,0)= -(b0 * c1); num_S(1,1)=  q; }
  else               { num_S(1,0)=   b0 * c1 ; num_S(1,1)= -q; }
  Q << num_S; rho << den_S;
  subresultant_compose (R, num_S, den_S);
  if (N(t1) == 0 || k < n - N(t1) + 1) return;
  C det_R= R(0,0)[0] * R(1,1)[0] - R(1,0)[0] * R(0,1)[0];
  half_subresultant_rec (t0, t1, false, det_R,
                         num_S, den_S, k_out, Q, rho, k_rem);
  subresultant_compose (R, num_S, den_S); }

TMPLP static Polynomial
_half_gcd (const Polynomial& f, const Polynomial& g,
           const Polynomial& first_prem) {
  // Special support for multi_gcd below: the first pseudo
  // remainder is given as an argument.
  // The gcd of the contents is not taken into account within this function.
  nat n= N(f), m= N(g);
  VERIFY (m > 1 && n >= m, "bad input sizes");

  C d= binpow (g[m-1], n-m);
  Polynomial q, s0= g, s1= first_prem;
  nat m0= N(s0), m1= N(s1); 
  if (m1 == 0) return g;
  if (m1 == 1) return s1;

  C c0= s0[m0-1], c1= s1[m1-1];
  C a= defected_pow (c1, d, m0-m1-1);
  Polynomial t0= (m0 == m1+1) ? s1 : ((a * s1) / d);
  C b0= t0[m1-1];
  Polynomial t1= rem (b0 * c1 * s0, s1, q);
  if (N(t1) == 0) return s1;
  if (N(t1) == 1) return t1;
  matrix<Polynomial> num_S (Polynomial (C(0)), 2, 2); C den_S (d * c0);
  t1 /= ((m0-m1-1) & 1) ? -den_S : den_S;
                       num_S(0,0)= 0         ; num_S(0,1)= a * c0;
  if ((m0-m1-1) & 1) { num_S(1,0)= -(b0 * c1); num_S(1,1)=  q; }
  else               { num_S(1,0)=   b0 * c1 ; num_S(1,1)= -q; }

  C den= (- g[m-1] * d * (num_S(0,0)[0] * num_S(1,1)[0]
                          - num_S(1,0)[0] * num_S(0,1)[0])) / den_S;
  vector<matrix<Polynomial> > Q; vector<C> rho;
  matrix<Polynomial> num_T (Polynomial (C(0)), 2, 2); C den_T; nat k_out;
  half_subresultant_rec (t0, t1, false, den,
                         num_T, den_T, k_out, Q, rho, N(t1));
  Polynomial h;
  nat k= N(t1) - degree (num_T(1,0));
  h= range (range (num_T(1,0), 0, k) * range (t0, 0, k) +
            range (num_T(1,1), 0, k) * range (t1, 0, k), 0, k);
  if (h != 0) return (den * h) / den_T;
  k= N(t1) - degree (num_T(0,0));
  h= range (range (num_T(0,0), 0, k) * range (t0, 0, k) +
            range (num_T(0,1), 0, k) * range (t1, 0, k), 0, k);
  VERIFY (h != 0, "bug");
  return (den * h) / den_T; }
  
public:

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) {
  // See documentation in subresultant_naive.hpp.
  VERIFY (degree (f) >= degree (g) && degree (g) > 0, "incorrect degrees");
  nat n= degree (f), m= degree (g), k;
  nat min_deg_res= (d1 != 0 || d2 != 0) ? dest_deg : m;
  for (k= min_deg_res; k+1 > 0; k--)
    if (k < N(res) && res[k] != 0) min_deg_res= k;
  nat min_deg_cof= (u1 != 0 || u2 != 0) ? dest_deg : m;
  for (k= min_deg_cof; k+1 > 0; k--)
    if (k < N(cof) && cof[k] != 0) min_deg_cof= k;
  nat min_deg_cog= (v1 != 0 || v2 != 0) ? dest_deg : m;
  for (k= min_deg_cog; k+1 > 0; k--)
    if (k < N(cog) && cog[k] != 0) min_deg_cog= k;
  nat min_deg_co= min (min_deg_cof, min_deg_cog);
  bool opt_res= res != 0 || d1 != 0 || d2 != 0;
  bool opt_cof= cof != 0 || u1 != 0 || u2 != 0;
  bool opt_cog= cog != 0 || v1 != 0 || v2 != 0;
  if (!(opt_res || opt_cof || opt_cog)) return;
  vector<matrix<Polynomial> > Q; vector<C> rho; matrix<Polynomial> R;
  half_subresultant (f, g, R, Q, rho, n + 1 - dest_deg);
  if (u1 != 0) u1= R(0,0); if (v1 != 0) v1= R(0,1);
  if (u2 != 0) u2= R(1,0); if (v2 != 0) v2= R(1,1);
  if (d1 != 0) {
    k= N(g) - degree (R(0,0));
    d1= range (range (R(0,0), 0, k) * range (f, 0, k) +
               range (R(0,1), 0, k) * range (g, 0, k), 0, k);
  }
  if (d2 != 0) {
    k= N(g) - degree (R(1,0));
    d2= range (range (R(1,0), 0, k) * range (f, 0, k) +
               range (R(1,1), 0, k) * range (g, 0, k), 0, k);
  }
  R= identity_matrix<Polynomial> (2);
  vector<Polynomial> r= vec<Polynomial> (f, g);
  vector<nat> nn= vec<nat> (N(g));
  // nn is the sequence of the successive sizes.
  for (nat i= 1; i < N(Q); i++) nn << (nn[i-1] - degree (Q[i](1,1)));
  nn << 0; nat l0= 0, l1; k= m;
  while (k > 0) {
    while ((k-1 >= N(res) || res[k-1] == 0) &&
           (k-1 >= N(cof) || cof[k-1] == 0) &&
           (k-1 >= N(cog) || cog[k-1] == 0) && k > 0) k--;
    if (k == 0) break;
    l1= find (nn, k);
    // If k is found then res[k-1] is regular.
    if (l1 >= N(nn)) {
      l1= find (nn, k+1);
      // If k+1 is found then res[k] is regular and res[k-1]
      // is defective, otherwise res[k-1] is zero.
      if (l1 >= N(nn)) {
        if (k-1 < N(res)) res[k-1] = 0;
        if (k-1 < N(cof)) cof[k-1] = 0;
        if (k-1 < N(cog)) cog[k-1] = 0;
        k--; continue;
      }
      k++;
    }
    subresultant_compose (R, range (Q, l0, l1+1), range (rho, l0, l1+1),
                          k-1 >= min_deg_co ? 0 : k);
    if (k-1 < N(res) && res[k-1] != 0)
      res[k-1]= range (range (R(0,0), 0, k) * range (f, 0, k) +
                       range (R(0,1), 0, k) * range (g, 0, k), 0, k);
    if (k-1 < N(cof) && cof[k-1] != 0) cof[k-1]= R(0,0);
    if (k-1 < N(cog) && cog[k-1] != 0) cog[k-1]= R(0,1);
    k--; if (k == 0) break;
    if (k-1 < N(res) && res[k-1] != 0)
      res[k-1]= range (range (R(1,0), 0, k) * range (f, 0, k) +
                       range (R(1,1), 0, k) * range (g, 0, k), 0, k);
    if (k-1 < N(cof) && cof[k-1] != 0) cof[k-1]= R(1,0);
    if (k-1 < N(cog) && cog[k-1] != 0) cog[k-1]= R(1,1);
    k--; l0= l1+1;
  }
}

#undef C
}; // implementation<polynomial_subresultant_base,V,polynomial_gcd_ring_dicho_inc<BV> >

/******************************************************************************
* Dichotomic multi-gcd
******************************************************************************/

template<typename V, typename W>
struct implementation<polynomial_evaluate,V,polynomial_gcd_ring_dicho_inc<W> >:
  public implementation<polynomial_evaluate,V,W>
{
  typedef implementation<polynomial_subresultant_base,V,
                         polynomial_gcd_ring_dicho_inc<W> > Pol;

#define C typename scalar_type_helper<Polynomial>::val
private:
TMPLP static Polynomial
gcd_with_first_prem (const Polynomial& P1, const Polynomial& P2,
                     const Polynomial& R) {
  // This is a special hack: the pseudoremainder R of P1 and P2
  // is forwarded to the subresultant computation.
  // R is only used when deg(P1) >= deg (P2).
  nat n1= N(P1), n2= N(P2);
  if (n1 < n2) return gcd_op::op (P2, P1);
  if (n2 == 0) { return P1; }
  if (n1 == 0) { return P2; }
  if (R == 0) {
    C c2; Polynomial T= primitive_part (P2, c2);
    return T * gcd_op::op (contents (P1), c2);
  }
  Polynomial G= Pol::_half_gcd (P1, P2, R);
  return primitive_part (G) * gcd_op::op (contents (P1), contents (P2));
}
#undef C

public:
TMPLP static vector<Polynomial>
multi_gcd (const Polynomial& P, const vector<Polynomial>& Q) {
  vector<Polynomial> R= prem (P, Q);
  for (nat i = 0; i < N(Q); i++)
    R[i]= gcd_with_first_prem (P, Q[i], R[i]);
  return R;
}

}; // implementation<polynomial_evaluate,V,polynomial_gcd_ring_dicho_inc<W> >

#undef TMPL
#undef TMPLP
} // namespace mmx
#endif //__MMX__POLYNOMIAL_RING_DICHO__HPP