include/algebramix/polynomial_ring_naive.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : polynomial_ring_naive.hpp
* DESCRIPTION: naive polynomial arithmetic over rings
* COPYRIGHT  : (C) 2008  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_NAIVE__HPP
#define __MMX__POLYNOMIAL_RING_NAIVE__HPP
#include <algebramix/matrix.hpp>
#include <algebramix/matrix_ring_naive.hpp>
#include <algebramix/polynomial_naive.hpp>

namespace mmx {
#define TMPL    template<typename Polynomial> 
#define C       typename scalar_type_helper<Polynomial>::val
#define Matrix  matrix<C,matrix_ring_naive<matrix_naive> >

/******************************************************************************
* Naive variant for polynomials over rings without division
******************************************************************************/

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

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

/******************************************************************************
* Naive variant for polynomials over rings without division, but with a gcd
******************************************************************************/

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

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

DEFINE_VARIANT_1 (typename V, V,
                  polynomial_gcd_ring_naive,
                  polynomial_gcd_ring_naive_inc<
                    polynomial_ring_naive<V> >)

/******************************************************************************
* Ducos optimization of naive variant over rings with a gcd
******************************************************************************/

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

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

DEFINE_VARIANT_1 (typename V, V,
                  polynomial_gcd_ring_ducos,
                  polynomial_gcd_ring_ducos_inc<
                    polynomial_gcd_ring_naive<V> >)

/******************************************************************************
* Naive but general implementation of subresultants
******************************************************************************/

template<typename V, typename W>
struct implementation<polynomial_subresultant_base,V,
                      polynomial_ring_naive<W> >:
  public implementation<polynomial_divide,V>
{
private:

TMPL static Matrix
sylvester_matrix (const Polynomial& f, const Polynomial& g, int k, int l) {
  // If l = k then it returns the matrix of the following map:
  //   P_{dp - k} x P_{dq - k} --> P_{dp + dq - 2k}
  //                   (u, v) :--> (u * f + v * g) quo x^k,
  // in the canonical basis 1, x, x^2, etc.
  // In general, the determinant of the matrix returned is the coefficient
  // of degree l of the kth subresultant.
  int n = max (deg (f), 0), m = max (deg (g), 0), k2= 2*k;
  ASSERT (k >= 0 && k <= min (n, m), "out of range");
  ASSERT (l >= 0 && l <= k, "out of range");
  Matrix S (C (0), n + m - k2, n + m - k2);
  for (int j= 0; j < m - k; j++)
    S (0, j) = f[l - j]; 
  for (int j= 0; j < n - k; j++)
    S (0, j + m - k) = g[l - j]; 
  for (int j= 0; j < m - k; j++) 
    for (int i= 1; i < n + m - k2; i++)
      S (i, j) = f[k + i - j]; 
  for (int j= 0; j < n - k; j++)
    for (int i= 1; i < n + m - k2; i++)
      S (i, j + m - k) = g[k + i - j]; 
  return S;
}

TMPL static inline C
subresultant (const Polynomial& f, const Polynomial& g, int k, int l) {
  int m= min (deg (f), deg (g));
  if (m < 0 ||
      k < 0 || k >= m ||
      l < 0 || l > k) return C(0);
  return det (sylvester_matrix (f, g, k, l));
}

TMPL static Polynomial
subresultant_co1 (const Polynomial& f, const Polynomial& g, int k) {
  int m= min (deg (f), deg (g));
  if (m < 0 ||
      k < 0 || k >= m) return Polynomial (C(0));
  Matrix S= sylvester_matrix (f, g, k, k);
  Polynomial cof(C(0));
  for (int l= 0; l < deg (g) - k; l++)
    cof += Polynomial (cofactor (S, 0, l), l);
  return cof;
}

TMPL static Polynomial
subresultant_co2 (const Polynomial& f, const Polynomial& g, int k) {
  int m= min (deg (f), deg (g));
  if (m < 0 ||
      k < 0 || k >= m) return Polynomial (C(0));
  Matrix S= sylvester_matrix (f, g, k, k);
  Polynomial cog(C(0));
  m= deg (g);
  for (int l= 0; l < deg (f) - k; l++)
    cog += Polynomial (cofactor (S, 0, m - k + l), l);
  return cog;
}

public:

TMPL static Polynomial
subresultant (const Polynomial& f, const Polynomial& g, int k) {
  Polynomial r(C(0));
  for (int l= 0; l <= k; l++)
    r += Polynomial (subresultant (f, g, k, l), l);
  return r;
}

public:

TMPL 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) {
  nat k, m= min (N (f), N (g)) - 1;
  if (m + 1 == 0) {
    d1= 0; u1= 0; v1= 0; 
    d2= 0; u2= 0; v2= 0; 
    return;
  }
  for (k= 0; k < m; k++) {
    if (k < N(res) && res[k] != 0) res[k]= subresultant     (f, g, k);
    if (k < N(cof) && cof[k] != 0) cof[k]= subresultant_co1 (f, g, k);
    if (k < N(cog) && cog[k] != 0) cog[k]= subresultant_co2 (f, g, k);
  }
  if (d1 != 0 || u1 != 0 || v1 != 0) {
    Polynomial d= subresultant (f, g, dest_deg);
    for (k= dest_deg; k < m;) {
      Polynomial tmp= C(0);
      if (d == 0) { k++; d= subresultant (f, g, k); continue; }
      if (k == m-1 || (tmp= subresultant (f, g, k+1)) != 0) break;
      else { k++; d= tmp; }
    }
    if (d1 != 0) d1= d;
    if (u1 != 0) u1= k < m ? subresultant_co1 (f, g, k) : 0;
    if (v1 != 0) v1= k < m ? subresultant_co2 (f, g, k) : 0;
  }
  if (d2 != 0 || u2 != 0 || v2 != 0) {
    Polynomial d= C(0);
    for (k= dest_deg - 1; k + 1 != 0 && d == 0; k--)
      d= subresultant (f, g, k);
    if (d2 != 0) d2= d;
    if (u2 != 0) u2= subresultant_co1 (f, g, k+1);
    if (v2 != 0) v2= subresultant_co2 (f, g, k+1); } }

}; // implementation<polynomial_subresultant_base,V,polynomial_ring_naive<W> >

/******************************************************************************
* General implementation of subresultants via pseudo-remaindering
******************************************************************************/

template<typename V, typename W>
struct implementation<polynomial_subresultant_base,V,
                      polynomial_gcd_ring_naive_inc<W> > {
// Lionel Ducos, "Optimizations of the subresultant algorithm",
// JPAA, 1998 (p. 150).

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

TMPL static inline Polynomial
defected_prem (const Polynomial& s_n, const Polynomial& s_n_1,
               const Polynomial& s_m, const C& x,
               Polynomial& u, const Polynomial& cof_n,
               const Polynomial& cof_n_1, const Polynomial& cof_m,
               bool opt_cof,
               Polynomial& v, const Polynomial& cog_n,
               const Polynomial& cog_n_1, const Polynomial& cog_m,
               bool opt_cog) {
  (void) cof_m; (void) cog_m;
  int n= degree (s_n), m= degree (s_n_1);
  C c_m= s_m[degree (s_m)], y= s_n_1[m];
  Polynomial r= s_n;
  if (opt_cof) u = y * cof_n - lshiftz (r[n] * cof_n_1, n-m);
  if (opt_cog) v = y * cog_n - lshiftz (r[n] * cog_n_1, n-m);
               r = y * r     - lshiftz (r[n] *   s_n_1, n-m);
  for (int k= n-1; k >= m; k--) {
    if (opt_cof) u = (y * u - lshiftz (r[k] * cof_n_1, k-m)) / x;
    if (opt_cog) v = (y * v - lshiftz (r[k] * cog_n_1, k-m)) / x;
                 r = (y * r - lshiftz (r[k] *   s_n_1, k-m)) / x;
  }
  y = ((n-m+1) & 1) ? -s_n[n] : s_n[n];
               r /= y;
  if (opt_cof) u /= y;
  if (opt_cog) v /= y;
  return r;
}

TMPL 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) {
  VERIFY (degree (f) >= degree (g) && degree (g) > 0, "incorrect degrees");
  nat n= degree (f), m= degree (g), k;
  bool b= false;
  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= min (min_deg_res, 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;
  if (m < dest_deg) { d1= 0; u1= 0; v1= 0; d2= 0; u2= 0; v2= 0; b= true; }
  C d= binpow (g[m], n-m);
  Polynomial q, s_n= g, s_n_1= prem (f, g, q);
  Polynomial cof_n= 0, cog_n= 1, cof_n_1= g[m] * d, cog_n_1= -q;
  n= degree (s_n); m= degree (s_n_1);
  if (n-1 < N(res) && res[n-1] != 0) res[n-1]=   s_n_1;
  if (n-1 < N(cof) && cof[n-1] != 0) cof[n-1]= cof_n_1;
  if (n-1 < N(cog) && cog[n-1] != 0) cog[n-1]= cog_n_1;
  if (m+1 < dest_deg+1 && !b) {
    if (d1 != 0) d1= 0;
    if (u1 != 0) u1= 0;
    if (v1 != 0) v1= 0;
    if (d2 != 0) d2=   s_n_1;
    if (u2 != 0) u2= cof_n_1;
    if (v2 != 0) v2= cog_n_1; b= true;
  } 
  C t0, t1;
  Polynomial s_m  , cof_m  , cog_m;
  Polynomial s_m_1, cof_m_1, cog_m_1;
  while (true) {
    if (m+1 < dest_deg+1 && !b) {
      if (d1 != 0) d1=   s_n;
      if (u1 != 0) u1= cof_n;
      if (v1 != 0) v1= cog_n;
      if (d2 != 0) d2=   s_n_1;
      if (u2 != 0) u2= cof_n_1;
      if (v2 != 0) v2= cog_n_1; b= true;
    } 
    for (k= m+1; k+1 < n; k++) {
      if (k < N(res) && res[k] != 0) res[k] = 0;
      if (k < N(cof) && cof[k] != 0) cof[k] = 0;
      if (k < N(cog) && cog[k] != 0) cog[k] = 0;
    }
    t0= binpow (s_n_1[m], n-m-1); t1= binpow (d, n-m-1);
                   s_m= (t0 *   s_n_1) / t1;
    if (opt_cof) cof_m= (t0 * cof_n_1) / t1;
    if (opt_cog) cog_m= (t0 * cog_n_1) / t1;
    if (m < N(res) && res[m] != 0) res[m] =   s_m;
    if (m < N(cof) && cof[m] != 0) cof[m] = cof_m;
    if (m < N(cog) && cog[m] != 0) cog[m] = cog_m;
    if (m+1 < min_deg+1) break;
    t0 *= square_op::op (s_n_1[m]); t1 *= ((n-m) & 1) ? s_n[n] * d : -s_n[n] * d;
                   s_m_1= prem (s_n, s_n_1, q)       / t1;
    if (opt_cof) cof_m_1= (t0 * cof_n - q * cof_n_1) / t1;
    if (opt_cog) cog_m_1= (t0 * cog_n - q * cog_n_1) / t1;
    if (m-1 < N(res) && res[m-1] != 0) res[m-1]=   s_m_1;
    if (m-1 < N(cof) && cof[m-1] != 0) cof[m-1]= cof_m_1;
    if (m-1 < N(cog) && cog[m-1] != 0) cog[m-1]= cog_m_1;
    s_n= s_m; s_n_1= s_m_1; d= s_n[m];
    n= degree (s_n); m= degree (s_n_1);
    if (opt_cof) { cof_n= cof_m; cof_n_1= cof_m_1; }
    if (opt_cog) { cog_n= cog_m; cog_n_1= cog_m_1; }
    opt_cof= m >= min_deg_cof;
    opt_cog= m >= min_deg_cog;
  }
}

}; // implementation<polynomial_subresultant_base,V,
   //                polynomial_gcd_ring_naive_inc<W> >

/******************************************************************************
* Optimized subresultants following Ducos
******************************************************************************/

template<typename V, typename BV>
struct implementation<polynomial_subresultant_base,V,
                      polynomial_gcd_ring_ducos_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);
}

TMPL static inline Polynomial
defected_prem (const Polynomial& s_n, const Polynomial& s_n_1,
               const Polynomial& s_m, const C& x,
               Polynomial& u, const Polynomial& cof_n,
               const Polynomial& cof_n_1, const Polynomial& cof_m,
               bool opt_cof,
               Polynomial& v, const Polynomial& cog_n,
               const Polynomial& cog_n_1, const Polynomial& cog_m, 
               bool opt_cog) {
  // Ibid. p. 153.
  int n= degree (s_n), m= degree (s_n_1);
  C c_n= s_n[n], c_n_1= s_n_1[m], c_m= s_m[m];
  // The code below corresponds to (discarding the cofactors):
  // r= rem (c_m * c_n_1 * s_n, s_n_1) / x;
  // return ((n-m+1) & 1) ? r / (-s_n[n]) : r / s_n[n];
  Polynomial t, r, h= range (s_m, 0, m), hf= cof_m, hg= cog_m;
  Polynomial d= c_m * range (s_n, 0, m) - s_n[m] * h,
    df= c_m * cof_n - s_n[m] * hf, dg= c_m * cog_n - s_n[m] * hg;
  for (int i= m + 1; i < n; i++) {
    if (opt_cof) {
      hf= lshiftz (hf, 1) - (h[m-1] * cof_n_1) / c_n_1;
      df -= s_n[i] * hf;
    }
    if (opt_cog) {
      hg= lshiftz (hg, 1) - (h[m-1] * cog_n_1) / c_n_1;
      dg -= s_n[i] * hg;
    }
    h= lshiftz (h, 1) - (h[m-1] * s_n_1) / c_n_1;
    d -= s_n[i] * h;
  }
  if (opt_cof) u= c_n_1 * (df / c_n)
                   - (c_n_1 * lshiftz (hf, 1) - h[m-1] * cof_n_1);
  if (opt_cog) v= c_n_1 * (dg / c_n)
                   - (c_n_1 * lshiftz (hg, 1) - h[m-1] * cog_n_1);
  r= c_n_1 * (d / c_n) - (c_n_1 * lshiftz (h, 1) - h[m-1] * s_n_1);
  C y= ((n-m+1) & 1) ? -x : x;
               r /= y;
  if (opt_cof) u /= y;
  if (opt_cog) v /= y;
  return r;
}

TMPL 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) {
  VERIFY (degree (f) >= degree (g) && degree (g) > 0, "incorrect degrees");
  nat n= degree (f), m= degree (g), k;
  bool b= false;
  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= min (min_deg_res, 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;
  if (m < dest_deg) { d1= 0; u1= 0; v1= 0; d2= 0; u2= 0; v2= 0; b= true; }
  C d= binpow (g[m], n-m);
  Polynomial q, s_n= g, s_n_1= prem (f, g, q);
  Polynomial cof_n= 0, cof_n_1= g[m] * d, cog_n= 1, cog_n_1= -q;
  n= degree (s_n); m= degree (s_n_1);
  if (n-1 < N(res) && res[n-1] != 0) res[n-1]=   s_n_1;
  if (n-1 < N(cof) && cof[n-1] != 0) cof[n-1]= cof_n_1;
  if (n-1 < N(cog) && cog[n-1] != 0) cog[n-1]= cog_n_1;
  if (m+1 < dest_deg+1 && !b) {
    if (d1 != 0) d1= 0;
    if (u1 != 0) u1= 0;
    if (v1 != 0) v1= 0;
    if (d2 != 0) d2=   s_n_1;
    if (u2 != 0) u2= cof_n_1;
    if (v2 != 0) v2= cog_n_1; b= true;
  } 
  C t0, t1;
  Polynomial s_m  , cof_m  , cog_m;
  Polynomial s_m_1, cof_m_1, cog_m_1;
  while (true) {
    if (m+1 < dest_deg+1 && !b) {
      if (d1 != 0) d1=   s_n;
      if (u1 != 0) u1= cof_n;
      if (v1 != 0) v1= cog_n;
      if (d2 != 0) d2=   s_n_1;
      if (u2 != 0) u2= cof_n_1;
      if (v2 != 0) v2= cog_n_1; b= true;
    } 
    for (k= m+1; k+1 < n; k++) {
      if (k < N(res) && res[k] != 0) res[k] = 0;
      if (k < N(cof) && cof[k] != 0) cof[k] = 0;
      if (k < N(cog) && cog[k] != 0) cog[k] = 0;
    }
    if (n == m+1) {
      s_m= s_n_1;
      if (opt_cof) cof_m= cof_n_1;
      if (opt_cog) cog_m= cog_n_1;
    }
    else {
      C a= defected_pow (s_n_1[m], d, n-m-1);
      s_m= (a * s_n_1) / d;
      if (opt_cof) cof_m= (a * cof_n_1) / d;
      if (opt_cog) cog_m= (a * cog_n_1) / d;
    }
    if (m < N (res) && res[m] != 0) res[m] = s_m;
    if (m < N (cof) && cof[m] != 0) cof[m] = cof_m;
    if (m < N (cog) && cog[m] != 0) cog[m] = cog_m;
    if (m+1 < min_deg+1) break;
    s_m_1= defected_prem (s_n, s_n_1, s_m, d,
                          cof_m_1, cof_n, cof_n_1, cof_m, opt_cof,
                          cog_m_1, cog_n, cog_n_1, cog_m, opt_cog);
    if (m-1 < N (res) && res[m-1] != 0) res[m-1]= s_m_1;
    if (m-1 < N (cof) && cof[m-1] != 0) cof[m-1]= cof_m_1;
    if (m-1 < N (cog) && cog[m-1] != 0) cog[m-1]= cog_m_1;
    s_n= s_m; s_n_1= s_m_1; d= s_n[m];
    if (opt_cof) { cof_n= cof_m; cof_n_1= cof_m_1; }
    if (opt_cog) { cog_n= cog_m; cog_n_1= cog_m_1; }
    n= degree (s_n); m= degree (s_n_1);
    opt_cof= m >= min_deg_cof;
    opt_cog= m >= min_deg_cog;
  }
}

}; // implementation<polynomial_subresultant_base,V,
   //                polynomial_gcd_ring_ducos_inc<BV> >

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

template<typename V, typename W>
struct implementation<polynomial_gcd,V,polynomial_ring_naive<W> >:
  public implementation<polynomial_divide,V>
{
  typedef implementation<polynomial_subresultant,V> Pol;

TMPL static Polynomial
gcd (const Polynomial& P1, const Polynomial& P2) {
  nat n1= N(P1), n2= N(P2);
  C c1, c2;
  if (n2 == 0) { return P1; }
  if (n1 == 0) { return P2; }
  Polynomial Q1= primitive_part (P1, c1);
  Polynomial Q2= primitive_part (P2, c2);
  C c= gcd_op::op (c1, c2);
  Polynomial zero= promote (0, CF(P1));
  vector<Polynomial> T (zero, 0);
  Polynomial G(C(1)), H(C(1)), D(C(0));
  Pol::subresultant_sequence (Q1, Q2, T, T, T, G, D, D, H, D, D, 1);
  if (H != 0) { G= H; }
  if (G == 0) { if (n1 < n2) G= Q1; else G= Q2; }
  else { C c; G= primitive_part (G, c); }
  return G * c;
}

TMPL static Polynomial
gcd (const Polynomial& P1, const Polynomial& P2,
     Polynomial& U1, Polynomial& U2) {
  ERROR ("extended gcd is not defined in " *
         string (typeid (Polynomial).name ()));
}

TMPL static vector<Polynomial>
multi_gcd (const Polynomial& P, const vector<Polynomial>& Q) {
  return binary_map_scalar<gcd_op> (Q, P); }

TMPL static inline Polynomial
lcm (const Polynomial& P1, const Polynomial& P2) {
  Polynomial dummy (0);
  return (P1 == 0 && P2 == 0) ? 0 : P1 * (P2 / gcd (P1, P2));
}

TMPL 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), G= gcd (P, Q, inv_P, dummy);
  if (deg (G) != 0) inv_P= 0;
  return inv_P / G[0];
}

TMPL 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.
  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; }
  Polynomial S= Q, D (0);
  Polynomial zero= promote (0, CF(P));
  vector<Polynomial> T (zero, 0);
  Num= 1; Den= 1;
  Pol::subresultant_sequence (S, P, T, T, T, D, D, D, Num, D, Den, k);
  if (Num == 0) { Num= P; Den= 1; }
}

}; // implementation<polynomial_gcd,V,polynomial_ring_naive<W> >

#undef Matrix
#undef TMPL
#undef C
} // namespace mmx
#endif //__MMX__POLYNOMIAL_RING_NAIVE__HPP

---