include/algebramix/series_carry_linear_algebra.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : series_carry_linear_algebra.hpp
* DESCRIPTION: linear system solving for lazy p-adic numbers
* COPYRIGHT  : (C) 2011 Jeremy Berthomieu and Romain Lebreton
*******************************************************************************
* 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_SERIES_CARRY_LINEAR_ALGEBRA_HPP
#define __MMX_SERIES_CARRY_LINEAR_ALGEBRA_HPP
#include <algebramix/series_carry.hpp>
#include <algebramix/series_matrix.hpp>

namespace mmx {
#define TMPL template<typename M, typename V>
#define C typename M::C
#define Series series<M,V>
#define Vector vector<M>
#define Matrix matrix<M>
#define Series_rep series_rep<M,V>
#define Series_vector series<Vector,V>
#define Series_vector_rep series_rep<Vector,V>
#define Vector_series vector<Series>
#define Series_matrix series<Matrix,V>
#define Series_matrix_rep series_rep<Matrix,V>
#define Matrix_series matrix<Series>

/******************************************************************************
* Specific converters
******************************************************************************/

//as_vector already exists in series_vector

// TMPL Vector_series
// asvector (const Series_vector& f, nat n) {
//   Vector_series v (Series (get_format1 (CF(f))), n);
//   for (nat i=0; i<n; i++)
//     v[i]= access (f, i);
//   return v;
// }

// TMPL inline Vector_series
// asvector (const Series_vector& f) {
//   return asvector (f, N(f[0]));
// }

TMPL Matrix_series
asmatrix (const Series_matrix& f) {
  nat nr= rows (f[0]), nc= cols (f[0]);
  Matrix_series m (Series (0), nr, nc);
  for (nat i=0; i<nr; i++)
    for (nat j=0; j<nc; j++)
      m (i, j)= access (f, i, j);
  return m;
}

TMPL Series_matrix
asseries (const Matrix_series& m) {
  return (series_rep<Matrix,V>*) 
    new matrix_series_rep<M,V,typename Matrix_variant(M)> (m);
}

TMPL Series_vector
asseries (const Vector_series& v) {
  return (series_rep<Vector,V>*) 
    new vector_series_rep<M,V,typename Vector_variant(M)> (v);
}

/******************************************************************************
* Specific shifts and accessors
******************************************************************************/

// TMPL
// class lshiftz_series_vector_rep: public series_rep<Vector,V> {
//   const Series_vector f;
//   const nat r;
//   const int shift;
// public:
//   lshiftz_series_vector_rep (const Series_vector& f2, nat r2, int shift2):
//     series_rep<Vector,V> (CF(f2)),
//     f (f2),  r (r2), shift (shift2) {}
//   syntactic expression (const syntactic& z) const {
//     return lshiftz (flatten (f, z), flatten (shift)); }
//   void Increase_order (nat l) {
//     series_rep<Vector,V>::Increase_order (l);
//     increase_order (f, max (0, ((int) l) - shift)); }
//   Vector next () { return shift > ((int) this->n)? 
//       Vector (M (), r): f[this->n - shift]; }
// };

// TMPL inline Series_vector
// lshiftz_series_vector (const Series_vector& f, const nat& r,
//                     const int& shift=1) {
//   if (is_exact_zero (f)) {
//     Series zero (get_format1 (CF(f)));
//     return as_series (Vector_series (zero, r));
//   }
//   return Series_vector ((series_rep<Vector,V>*) 
//                      new lshiftz_series_vector_rep<M,V> (f, r, shift));
// }

TMPL
class lshiftz_series_matrix_rep: public series_rep<Matrix,V> {
  const Series_matrix f;
  const nat r,c;
  const int shift;
public:
  lshiftz_series_matrix_rep (const Series_matrix& f2, nat r2, nat c2,
                             int shift2):
    series_rep<Matrix,V> (CF(f2)),
    f (f2),  r (r2), c (c2), shift (shift2) {}
  syntactic expression (const syntactic& z) const {
    return lshiftz (flatten (f, z), flatten (shift)); }
  void Increase_order (nat l) {
    series_rep<Matrix,V>::Increase_order (l);
    increase_order (f, max (0, ((int) l) - shift)); }
  Matrix next () { return shift > ((int) this->n)? 
      Matrix (M (), r, c): f[this->n - shift]; }
};

TMPL inline Series_matrix
lshiftz_series_matrix (const Series_matrix& f, const nat& r, const nat& c,
                       const int& shift=1) {
  if (is_exact_zero (f)) {
    Series zero (get_format1 (CF(f)));
    return as_series (Matrix_series (zero, r, c));
  }
  return Series_matrix ((series_rep<Matrix,V>*) 
                        new lshiftz_series_matrix_rep<M,V> (f, r, c, shift));
}

TMPL inline Vector_series
lshiftz (const Vector_series& v, const int& shift=1) {
  nat n= N (v);
  Vector_series lv (Series(), n);
  for (nat i=0; i<n; i++)
    lv[i] = lshiftz (v[i], shift);
  return lv;
}

TMPL inline Vector_series
rshiftz (const Vector_series& v, const int& shift=1) {
  return lshiftz (v, -shift);
}

TMPL inline Matrix_series
lshiftz (const Matrix_series& m, const int& shift=1) {
  nat c= cols (m), r= rows (m);
  Matrix_series lm (Series(), r, c);
  for (nat i=0; i<r; i++)
    for (nat j=0; j<c; j++)
      lm(i,j) = lshiftz (m(i,j), shift);
  return lm;
}

TMPL inline Matrix_series
rshiftz (const Matrix_series& m, const int& shift=1) {
  return lshiftz (m, -shift);
}

TMPL inline Vector
access (const Vector_series& v, nat n) {
  nat nv= N (v);
  Vector av (M(), nv);
  for (nat i=0; i<nv; i++)
    av[i] = v[i][n];
  return av;
}

TMPL inline Matrix
access (const Matrix_series& m, nat n) {
  nat c= cols (m), r= rows (m);
  Matrix am (M(), r, c);
  for (nat i=0; i<r; i++)
    for (nat j=0; j<c; j++)
      am(i,j) = m(i,j)[n];
  return am;
}

/******************************************************************************
* Multiplication quotient and remainder between a matrix and a vector
******************************************************************************/

struct series_vector_divide {};

template<typename U>
struct implementation<series_vector_divide,U,series_naive> {

template<typename M, typename V, typename X>
class vector_carry_mul_quo_series_rep: public series_rep<Vector,V> {
  typedef Lift_type(M) L;
protected:
  Series_vector f;
  matrix<L> x;
  vector<L> fn, c;
  Vector q;
  L p;
  nat cx, rx, nf;
public:
  inline vector_carry_mul_quo_series_rep (const matrix<X>& x2, 
                                          const Series_vector& f2):
    series_rep<Vector,V> (get_format (vec<Vector> 
                                      (get_sample (CF(f2))))),
    f(f2), cx(cols(x2)), rx (rows (x2)), nf (N (f2[0]))
  {
    ASSERT ((cx == rx) && (cx == nf), "Matrix and vector sizes do not match");
    p = L (* M::get_modulus());
    c = vector<L> (L (0), nf);
    fn= vector<L> (L (0), nf);
    x = matrix<L> (L (0), rx, cx);
    for (nat i=0; i < rx; i++)
      for (nat j=0; j < cx; j++)
        x(i,j)= L (* x2(i,j));
  }
  syntactic expression (const syntactic& z) const {
    return apply (syntactic ("vector_carry_mul_quo"), flatten (f, z),
                  flatten (x)); }
  virtual void Increase_order (nat l) {
    series_rep<Vector,V>::Increase_order (l);
    increase_order (f, l); }
  virtual Vector next () {
    for (nat i=0; i < nf; i++)
      fn[i]= L (* f[this->n][i]);
    fn= quo (x * fn, p) + c; // TODO Use Modulus as in series_carry_naive
                             // so that rem/quo are faster when modulus_preinv
    q = as<Vector> (rem (fn, p));
    c = quo (fn, p);
    return q; }
};

TMPL
class vector_carry_mul_rem_series_rep: public series_rep<Vector, V> {
protected:
  Matrix m;
  Series_vector sv;
public:
  inline vector_carry_mul_rem_series_rep (const Matrix& m2,
                                          const Series_vector& sv2):
    series_rep<Vector, V> (CF(m2)), m (m2), sv(sv2) {}
  syntactic expression (const syntactic& z) const {
    return apply (apply (syntactic ("vector_carry_mul_rem"), flatten (m),
                         flatten (sv,z))); }
  virtual void Increase_order (nat l) {
    series_rep<Vector, V>::Increase_order (l);
    increase_order (sv, l); }
  virtual Vector next () {
    return m*sv[this->n]; }
};

TMPL
class ldiv_sc_vec_series_rep: public recursive_series_rep<Vector,V> {
protected:
  Vector_series f;
  Matrix c;
  nat nf, rc, cc;
  Matrix inv;
  bool init;
  Vector vinit;
public:
  ldiv_sc_vec_series_rep (const Matrix& c2, const Vector_series& f2):
    recursive_series_rep<Vector,V> (get_vector_format (f2)),
    f(f2), c(c2), nf(N(f2)), rc(rows(c2)), cc(cols(c2)), inv(invert(c)),
    init(false) {}
  ldiv_sc_vec_series_rep (const Matrix& c2, const Vector_series& f2,
                          const Matrix& inv2):
    recursive_series_rep<Vector,V> (get_vector_format (f2)),
    f(f2), c(c2), nf(N(f2)), rc(rows (c2)),cc(cols(c2)), inv(inv2),
    init(false) {}
  ldiv_sc_vec_series_rep (const Matrix& c2, const Vector_series& f2,
                          const Matrix& inv2, const Vector& vinit2):
    recursive_series_rep<Vector,V> (get_vector_format (f2)),
    f(f2), c(c2), nf(N(f2)), rc(rows(c2)), cc(cols(c2)), inv(inv2), init(true),
    vinit(vinit2) {}
  syntactic expression (const syntactic& z) const {
    return invert (flatten (c)) * flatten (f,z); }
  virtual void Increase_order (nat l) {
    recursive_series_rep<Vector,V>::Increase_order (l);
    for (nat i = 0; i < nf; i++)
      increase_order (f[i],l);
  }
  Series_vector initialize () {
    ASSERT ((cc == rc) && (cc == nf), "Matrix and vector sizes do not match");
    if (init) this->initial (0)= vinit;
    else this->initial (0)= inv * access (f, 0);
    Series_vector ls_me = lshiftz_series_vector (this -> me (), nf);
    Vector_series tmp= f - as_vector (lmul_quo (c, ls_me));
    return lmul_rem (inv, asseries (tmp));
  }
};

TMPL
class ldiv_vec_series_rep: public recursive_series_rep<Vector,V> {
protected:
  Vector_series f;
  Matrix_series m;
  nat nf, rm, cm;
  bool init;
  Vector vinit;
public:
  ldiv_vec_series_rep (const Matrix_series& m2, const Vector_series& f2):
    recursive_series_rep<Vector,V> (get_vector_format (f2)), f(f2), m(m2),
    nf(N(f2)), rm(rows(m2)), cm(cols(m2)), init(false) {}
  ldiv_vec_series_rep (const Matrix_series& m2, const Vector_series& f2,
                       const Vector& vinit2):
    recursive_series_rep<Vector,V> (get_vector_format (f2)), f(f2), m(m2),
    nf(N(f2)), rm(rows(m2)), cm(cols(m2)), init(true), vinit (vinit2) {}
  syntactic expression (const syntactic& z) const {
    return invert (flatten (m)) * flatten (f,z); }
  virtual void Increase_order (nat l) {
    recursive_series_rep<Vector,V>::Increase_order (l);
    for (nat i = 0; i < nf; i++) {
      increase_order (f[i],l);
      for (nat j = 0; j < cm; j++)
        increase_order (m(i,j),l);
    }
  }
  Series_vector initialize () {
    ASSERT ((cm == rm) && (cm == nf), "Matrix and vector sizes do not match");
    Matrix inv= invert (access (m, 0));
    if(init) this->initial (0)= vinit;
    else this->initial (0)= inv * access (f, 0);
    Vector_series va = f - lshiftz (rshiftz (m) * as_vector (this-> me ()));
    return ldiv (access (m, 0), va, inv);
  }
};

TMPL static inline Series_vector
ser_lmul_quo_sc (const Matrix& c, const Series_vector& f) {
  return Series_vector ((series_rep<Vector,V> *)
                        new vector_carry_mul_quo_series_rep<M,V,M> (c, f));
}

TMPL static inline Series_vector
ser_lmul_rem_sc (const Matrix& m, const Series_vector& sv) {
  return Series_vector ((series_rep<Vector, V> *) 
                        new vector_carry_mul_rem_series_rep<M,V> (m, sv));
}

TMPL static inline Series_vector
ser_ldiv_sc (const Matrix& c, const Vector_series& f) {
  typedef ldiv_sc_vec_series_rep<M,V> Div_sc_rep;
  return recursive (Series_vector ((series_rep<Vector,V>*)
                                   new Div_sc_rep (c,f)));
}

TMPL static inline Series_vector
ser_ldiv_sc (const Matrix& c, const Vector_series& f, const Matrix& inv) {
  typedef ldiv_sc_vec_series_rep<M,V> Div_sc_rep;
  return recursive (Series_vector ((series_rep<Vector,V>*)
                                   new Div_sc_rep (c, f, inv)));
}

TMPL static inline Series_vector
ser_ldiv_sc (const Matrix& c, const Vector_series& f, const Matrix& inv,
             const Vector& init) {
  typedef ldiv_sc_vec_series_rep<M,V> Div_sc_rep;
  return recursive (Series_vector ((series_rep<Vector,V>*)
                                   new Div_sc_rep (c, f, inv, init)));
}

TMPL static inline Series_vector
ser_ldiv (const Matrix_series& m, const Vector_series& f) {
  typedef ldiv_vec_series_rep<M,V> Div_rep;
  return recursive (Series_vector ((series_rep<Vector,V>*) 
                                   new Div_rep (m, f)));
}

TMPL static inline Series_vector
ser_ldiv (const Matrix_series& m, const Vector_series& f, const Vector& init) {
  typedef ldiv_vec_series_rep<M,V> Div_rep;
  return recursive (Series_vector ((series_rep<Vector,V>*) 
                                   new Div_rep (m, f, init)));
}


}; // implementation<series_vector_divide,V,series_naive>

TMPL static inline Series_vector
lmul_quo (const Matrix& m, const Series_vector& sv) {
  typedef implementation<series_vector_divide,V> Ser;
  return Ser::template ser_lmul_quo_sc<M,V> (m, sv);
}

TMPL inline Series_vector
lmul_rem (const Matrix& m, const Series_vector& sv) {
  typedef implementation<series_vector_divide,V> Ser;
  return Ser::template ser_lmul_rem_sc<M,V> (m, sv);
}

TMPL static inline Series_vector
ldiv (const Matrix& m, const Vector_series& vs) {
  typedef implementation<series_vector_divide,V> Ser;
  return Ser::template ser_ldiv_sc<M,V> (m, vs);
}

TMPL static inline Series_vector
ldiv (const Matrix& m, const Vector_series& vs, const Matrix& inv) {
  typedef implementation<series_vector_divide,V> Ser;
  return Ser::template ser_ldiv_sc<M,V> (m, vs, inv);
}

TMPL static inline Series_vector
ldiv (const Matrix& m, const Vector_series& vs, const Matrix& inv,
      const Vector& init) {
  typedef implementation<series_vector_divide,V> Ser;
  return Ser::template ser_ldiv_sc<M,V> (m, vs, inv, init);
}

TMPL static inline Series_vector
ldiv (const Matrix_series& m, const Vector_series& vs) {
  typedef implementation<series_vector_divide,V> Ser;
  return Ser::template ser_ldiv<M,V> (m, vs);
}

TMPL static inline Series_vector
ldiv (const Matrix_series& m, const Vector_series& vs, const Vector& init) {
  typedef implementation<series_vector_divide,V> Ser;
  return Ser::template ser_ldiv<M,V> (m, vs, init);
}

/******************************************************************************
* Multiplication quotient and remainder between matrices
******************************************************************************/

struct series_matrix_divide {};

template<typename U>
struct implementation<series_matrix_divide,U,series_naive> {

template<typename M, typename V, typename X>
class matrix_carry_mul_quo_series_rep: 
    public series_rep<Matrix,V> {
  typedef Lift_type(M) L;
protected:
  Series_matrix f;
  matrix<L> x, fn, c;
  Matrix q;
  L p;
  nat cx, rx, cf, rf;
public:
  inline matrix_carry_mul_quo_series_rep (const matrix<X>& x2, 
                                          const Series_matrix& f2):
    series_rep<Matrix,V> (get_format (vec<Matrix> 
                                      (get_sample (CF(f2))))),
    f(f2), cx(cols(x2)), rx (rows (x2)), cf (cols (f2[0])), rf (rows (f2[0]))
  {
    ASSERT ((cx == rx) && (cx == rf), "Matrices sizes do not match");
    p = L (* M::get_modulus());
    c = matrix<L> (L (0), rf, cf);
    fn= matrix<L> (L (0), rf, cf);
    x = matrix<L> (L (0), rx, cx);
    for (nat i=0; i < rx; i++)
      for (nat j=0; j < cx; j++)
        x(i,j)= L (* x2(i,j));
  }
  syntactic expression (const syntactic& z) const {
    return apply (syntactic ("matrix_carry_mul_quo"), flatten (f, z),
                  flatten (x)); }
  virtual void Increase_order (nat l) {
    series_rep<Matrix,V>::Increase_order (l);
    increase_order (f, l); }
  virtual Matrix next () {
    for (nat i=0; i < rf; i++)
      for (nat j=0; j < cf; j++)
        fn(i,j)= L (* f[this->n](i,j));
    fn= quo (x * fn, p) + c;
    q = as<Matrix> (rem (fn, p));
    c = quo (fn, p);
    return q; }
};

TMPL
class matrix_carry_mul_rem_series_rep: public series_rep<Matrix, V> {
protected:
  Matrix m;
  Series_matrix sm;
public:
  inline matrix_carry_mul_rem_series_rep (const Matrix& m2,
                                          const Series_matrix& sm2):
    series_rep<Matrix, V> (CF(m2)), m (m2), sm(sm2) {}
  syntactic expression (const syntactic& z) const {
    return apply (apply (syntactic ("matrix_carry_mul_rem"), flatten (m),
                         flatten (sm,z))); }
  virtual void Increase_order (nat l) {
    series_rep<Matrix, V>::Increase_order (l);
    increase_order (sm, l); }
  virtual Matrix next () {
    return m*sm[this->n]; }
};
TMPL
class ldiv_sc_vec_series_rep: public recursive_series_rep<Vector,V> {
protected:
  Vector_series f;
  Matrix c;
  nat nf, rc, cc;
  Matrix inv;
public:
  ldiv_sc_vec_series_rep (const Matrix& c2, const Vector_series& f2):
    recursive_series_rep<Vector,V> (get_vector_format (f2)),
    f(f2), c(c2), nf (N(f2)), rc (rows (c2)), cc (cols (c2)), inv (invert (c))
  {}
  ldiv_sc_vec_series_rep (const Matrix& c2, const Vector_series& f2,
                          const Matrix& inv2):
    recursive_series_rep<Vector,V> (get_vector_format (f2)),
    f(f2), c(c2), nf (N(f2)), rc (rows (c2)), cc (cols (c2)), inv (inv2) {}
  syntactic expression (const syntactic& z) const {
    return invert (flatten (c)) * flatten (f,z); }
  virtual void Increase_order (nat l) {
    recursive_series_rep<Vector,V>::Increase_order (l);
    for (nat i = 0; i < nf; i++)
      increase_order (f[i],l);
  }
  Series_vector initialize () {
    ASSERT ((cc == rc) && (cc == nf), "Matrix and vector sizes do not match");
    this->initial (0)= inv * access (f, 0);
    Series_vector ls_me = lshiftz_series_vector (this -> me (), nf);
    Vector_series tmp= f - as_vector (lmul_quo (c, ls_me));
    return lmul_rem (inv, asseries (tmp));
  }
};

TMPL static inline Series_vector
ser_ldiv_sc_vec (const Matrix& c, const Vector_series& f) {
  typedef ldiv_sc_vec_series_rep<M,V> Div_sc_vec_rep;
  /*if (is_exact_zero (asseries (f)))
    return Vector_series (Series (CF(c)), cols(c));*/
  return recursive (Series_vector ((series_rep<Vector,V>*)
                                   new Div_sc_vec_rep (c,f)));
}

TMPL static inline Series_vector
ser_ldiv_sc_vec (const Matrix& c, const Vector_series& f, const Matrix& inv) {
  typedef ldiv_sc_vec_series_rep<M,V> Div_sc_vec_rep;
  /*if (is_exact_zero (asseries (f)))
    return Vector_series (Series (CF(c)), cols(c));*/
  return recursive (Series_vector ((series_rep<Vector,V>*)
                                   new Div_sc_vec_rep (c,f, inv)));
}

TMPL
class ldiv_vec_series_rep: public recursive_series_rep<Vector,V> {
protected:
  Vector_series f;
  Matrix_series m;
  nat nf, rm, cm;
public:
  ldiv_vec_series_rep (const Matrix_series& m2, const Vector_series& f2):
    recursive_series_rep<Vector,V> (get_vector_format (f2)), f(f2), m(m2),
    nf(N(f2)), rm(rows (m2)), cm(cols (m2)) {}
  syntactic expression (const syntactic& z) const {
    return invert (flatten (m)) * flatten (f,z); }
  virtual void Increase_order (nat l) {
    recursive_series_rep<Vector,V>::Increase_order (l);
    for (nat i = 0; i < nf; i++) {
      increase_order (f[i],l);
      for (nat j = 0; j < cm; j++)
        increase_order (m(i,j),l);
    }
  }
  Series_vector initialize () {
    ASSERT ((cm == rm) && (cm == nf), "Matrix and vector sizes do not match");
    Matrix inv= invert (access (m, 0));
    this->initial (0)= inv * access (f, 0);
    Vector_series va = f - lshiftz (rshiftz (m) * as_vector (this-> me ()));
    //Peut-on faire un produit semi-détendu ici ^ ?
    return ser_ldiv_sc_vec (access (m, 0), va, inv);
  }
};

TMPL static inline Series_vector
ser_ldiv_vec (const Matrix_series& m, const Vector_series& f) {
  typedef ldiv_vec_series_rep<M,V> Div_vec_rep;
  /*
  if (is_exact_zero (asseries (f))) {
    Series zero (get_format1 (CF(f)));
    return Vector_series (zero, cols(m), cols(f));
  }
  */
  return recursive (Series_vector ((series_rep<Vector,V>*) 
                                   new Div_vec_rep (m,f)));
}

/******************************************************************************
* Matrix equations resolution
******************************************************************************/

TMPL
class ldiv_sc_mat_series_rep: public recursive_series_rep<Matrix,V> {
protected:
  Matrix_series f;
  Matrix c;
  nat rf, cf, rc, cc;
  Matrix inv;
  bool init;
  Matrix minit;
public:
  ldiv_sc_mat_series_rep (const Matrix& c2, const Matrix_series& f2):
    recursive_series_rep<Matrix,V> (get_matrix_format (f2)), f(f2), c(c2),
    rf(rows(f2)), cf(cols(f2)), rc(rows(c2)), cc(cols(c2)), inv(invert(c)),
    init(false) {}
  ldiv_sc_mat_series_rep (const Matrix& c2, const Matrix_series& f2,
                          const Matrix& inv2):
    recursive_series_rep<Matrix,V> (get_matrix_format (f2)), f (f2), c (c2),
    rf(rows(f2)), cf(cols(f2)), rc(rows(c2)), cc(cols(c2)), inv(inv2),
    init(false) {}
  ldiv_sc_mat_series_rep (const Matrix& c2, const Matrix_series& f2,
                          const Matrix& inv2, const Matrix& minit2):
    recursive_series_rep<Matrix,V> (get_matrix_format (f2)), f(f2), c(c2),
    rf(rows(f2)), cf(cols(f2)), rc(rows(c2)), cc(cols(c2)), inv(inv2),
    init(true), minit(minit2) {}
  syntactic expression (const syntactic& z) const {
    return invert (flatten (c)) * flatten (f,z); }
  virtual void Increase_order (nat l) {
    recursive_series_rep<Matrix,V>::Increase_order (l);
    for (nat i = 0; i < rf; i++)
      for (nat j = 0; j < cf; j++)
        increase_order (f(i,j), l);
  }
  Series_matrix initialize () {
    ASSERT ((cc == rc) && (cc == rf), "Matrices sizes do not match");
    if (init) this->initial (0)= minit;
    else this->initial (0)= inv * access (f, 0);
    Series_matrix ls_me = lshiftz_series_matrix (this -> me (), rc, cf);
    Matrix_series tmp= f - asmatrix (lmul_quo (c, ls_me));
    return lmul_rem (inv, asseries (tmp));
  }
};
TMPL
class ldiv_mat_series_rep: public recursive_series_rep<Matrix,V> {
protected:
  Matrix_series f;
  Matrix_series m;
  nat rf, cf, rm, cm;
  bool init;
  Matrix minit;
public:
  ldiv_mat_series_rep (const Matrix_series& m2, const Matrix_series& f2):
    recursive_series_rep<Matrix,V> (get_matrix_format (m2)), f(f2), m(m2),
    rf(rows(f2)), cf(cols(f2)), rm(rows(m2)), cm(cols(m2)), init(false) {}
  ldiv_mat_series_rep (const Matrix_series& m2, const Matrix_series& f2,
                       const Matrix& minit2):
    recursive_series_rep<Matrix,V> (get_matrix_format (m2)), f(f2), m(m2),
    rf(rows(f2)), cf(cols(f2)), rm(rows(m2)), cm(cols(m2)), init(true),
    minit (minit2) {}
  syntactic expression (const syntactic& z) const {
    return invert (flatten (m)) * flatten (f,z); }
  virtual void Increase_order (nat l) {
    recursive_series_rep<Matrix,V>::Increase_order (l);
    for (nat i = 0; i < rf; i++) {
      for (nat j = 0; j < cf; j++)
        increase_order (f(i,j),l);
      for (nat j = 0; j < cm; j++)
        increase_order (m(i,j),l);
    }
  }
  Series_matrix initialize () {
    ASSERT ((cm == rm) && (cm == rf), "Matrices sizes do not match");
    Matrix inv= invert (access (m, 0));
    if (init) this->initial (0)= minit;
    else this->initial (0)= inv * access (f, 0);
    Matrix_series ma = f - lshiftz (rshiftz (m) * asmatrix (this-> me ()));
    return ldiv (access (m, 0), ma, inv);
  }
};

TMPL static inline Series_matrix
ser_lmul_quo_sc (const Matrix& c, const Series_matrix& f) {
  return Series_matrix ((series_rep<Matrix,V> *)
                        new matrix_carry_mul_quo_series_rep<M,V,M> (c, f));
}

TMPL static inline Series_matrix
ser_lmul_rem_sc (const Matrix& m, const Series_matrix& sm) {
  return Series_matrix ((series_rep<Matrix, V> *) 
                        new matrix_carry_mul_rem_series_rep<M,V> (m, sm));
}

TMPL static inline Series_matrix
ser_ldiv_sc (const Matrix& c, const Matrix_series& f) {
  typedef ldiv_sc_mat_series_rep<M,V> Div_sc_rep;
  return recursive (Series_matrix ((series_rep<Matrix,V>*)
                                   new Div_sc_rep (c, f)));
}

TMPL static inline Series_matrix
ser_ldiv_sc (const Matrix& c, const Matrix_series& f, const Matrix& inv) {
  typedef ldiv_sc_mat_series_rep<M,V> Div_sc_rep;
  return recursive (Series_matrix ((series_rep<Matrix,V>*)
                                   new Div_sc_rep (c, f, inv)));
}

TMPL static inline Series_matrix
ser_ldiv_sc (const Matrix& c, const Matrix_series& f, const Matrix& inv,
             const Matrix& init) {
  typedef ldiv_sc_mat_series_rep<M,V> Div_sc_rep;
  return recursive (Series_matrix ((series_rep<Matrix,V>*)
                                   new Div_sc_rep (c, f, inv, init)));
}

TMPL static inline Series_matrix
ser_ldiv (const Matrix_series& m, const Matrix_series& f) {
  typedef ldiv_mat_series_rep<M,V> Div_rep;
  return recursive (Series_matrix ((series_rep<Matrix,V>*) 
                                   new Div_rep (m, f)));
}

TMPL static inline Series_matrix
ser_ldiv (const Matrix_series& m, const Matrix_series& f, const Matrix& init) {
  typedef ldiv_mat_series_rep<M,V> Div_rep;
  return recursive (Series_matrix ((series_rep<Matrix,V>*) 
                                   new Div_rep (m, f, init)));
}

}; // implementation<series_matrix_carry_divide,series_naive>

TMPL static inline Series_matrix
lmul_quo (const Matrix& m, const Series_matrix& sm) {
  typedef implementation<series_matrix_divide,V> Ser;
  return Ser::template ser_lmul_quo_sc<M,V> (m, sm);
}

TMPL inline Series_matrix
lmul_rem (const Matrix& m, const Series_matrix& sm) {
  typedef implementation<series_matrix_divide,V> Ser;
  return Ser::template ser_lmul_rem_sc<M,V> (m, sm);
}

TMPL static inline Series_matrix
ldiv (const Matrix& m, const Matrix_series& ms) {
  typedef implementation<series_matrix_divide,V> Ser;
  return Ser::template ser_ldiv_sc<M,V> (m, ms);
}

TMPL static inline Series_matrix
ldiv (const Matrix& m, const Matrix_series& ms, const Matrix& inv) {
  typedef implementation<series_matrix_divide,V> Ser;
  return Ser::template ser_ldiv_sc<M,V> (m, ms, inv);
}

TMPL static inline Series_matrix
ldiv (const Matrix& m, const Matrix_series& ms, const Matrix& inv,
      const Matrix& init) {
  typedef implementation<series_matrix_divide,V> Ser;
  return Ser::template ser_ldiv_sc<M,V> (m, ms, inv, init);
}

TMPL static inline Series_matrix
ldiv (const Matrix_series& m, const Matrix_series& ms) {
  typedef implementation<series_matrix_divide,V> Ser;
  return Ser::template ser_ldiv<M,V> (m, ms);
}

TMPL static inline Series_matrix
ldiv (const Matrix_series& m, const Matrix_series& ms, const Matrix& init) {
  typedef implementation<series_matrix_divide,V> Ser;
  return Ser::template ser_ldiv<M,V> (m, ms, init);
}

/******************************************************************************
* Monoblock linear algebra
******************************************************************************/

#define TMPLB template<typename M, typename V, nat s, typename BV, nat t>
#define Monoblock_series series <Monoblock_modular,BV>
#define Monoblock_transformer series_carry_monoblock_transformer<M,V,s,BV>
#define Monoblock_modular modular<modulus<Lift_type(M)>,                \
                                  modular_global_series_carry_monoblock \
                                  <M,s,BV> >

template<typename U,typename W,nat s,typename BV,nat t>
struct implementation<series_vector_divide,U,
                      series_carry_monoblock<W,s,BV,t> > {

template<typename M, typename V>
class ldiv_vec_monoblock_series_rep: public Series_vector_rep {
protected:
  Vector_series f;
  Matrix_series m;
  nat nf, rm, cm;

  Monoblock_transformer blocker;
  bool h1_init;
  Series_vector h0, h1;  

public:
  ldiv_vec_monoblock_series_rep (const Matrix_series& m2,
                                 const Vector_series& f2):
    Series_vector_rep (get_vector_format (f2)), f(f2), m(m2), nf(N(f2)),
    rm(rows(m2)), cm(cols(m2)), h1_init(false) {
    h0 = ldiv (m, f);
  }
  syntactic expression (const syntactic& z) const {
    return invert (flatten (m)) * flatten (f,z); }
  virtual void Increase_order (nat l) {
    Series_vector_rep::Increase_order (l);
    if (l <= t)
      increase_order (h0, l);
    else {
      if (! h1_init) {
        typedef Monoblock_modular MM;
        typedef Lift_type (M) L;
        MM::set_modulus (binpow (L(* M::get_modulus ()), s));
        matrix<Monoblock_series> newm (Monoblock_series (), rm, rm);
        vector<Monoblock_series> newf (Monoblock_series (), nf);
        vector<MM> init (MM (0), rm);
        for (nat i=0; i<rm; i++) {
          newf[i]= to_monoblock (f[i], blocker);
          init[i]= as<L> (truncate (access (h0, i), s));
          for (nat j=0; j<rm; j++)
            newm(i,j)= to_monoblock (m(i,j), blocker);
        }
        vector<Monoblock_series> newvh1= as_vector (ldiv (newm, newf, init));
        Vector_series vh1 (Series(), nf);
        for (nat i=0; i<rm; i++)
          vh1[i]= from_monoblock (newvh1[i], blocker);
        h1= asseries (vh1);
        h1_init= true;
      }
      increase_order (h1, l); } }
  Vector next () {
    return this->n < t ? h0[this->n] : h1[this->n]; }
};

//The next functions should be inherited
TMPL static inline Series_vector
ser_lmul_quo_sc (const Matrix& c, const Series_vector& f) {
  typedef implementation<series_vector_divide,W> Ser;
  return Ser::ser_lmul_quo_sc (c, f); }

TMPL static inline Series_vector
ser_lmul_rem_sc (const Matrix& c, const Series_vector& f) {
  typedef implementation<series_vector_divide,W> Ser;
  return Ser::ser_lmul_rem_sc (c, f); }

TMPL static inline Series_vector
ser_ldiv_sc (const Matrix& c, const Vector_series& f) {
  typedef implementation<series_vector_divide,W> Ser;
  return Ser::ser_ldiv_sc (c, f); }

TMPL static inline Series_vector
ser_ldiv_sc (const Matrix& c, const Vector_series& f, const Matrix& inv) {
  typedef implementation<series_vector_divide,W> Ser;
  return Ser::ser_ldiv_sc (c, f, inv); }
// --inherited

TMPL static inline Series_vector
ser_ldiv (const Matrix_series& c, const Vector_series& f) {
  typedef ldiv_vec_monoblock_series_rep<M,W> Div_rep;
  series<vector<M>,W> s =
    series<vector<M>,W> ((series_rep<vector<M>,W>*) new Div_rep
                         (as<matrix<series<M,W> > > (c),
                          as<vector<series<M,W> > > (f)));
  return as<Series_vector> (s);
}

}; // implementation<series_vector_divide,U,series_carry_monoblock<W> >

template<typename U,typename W,nat s,typename BV,nat t>
struct implementation<series_matrix_divide,U,
                      series_carry_monoblock<W,s,BV,t> > {

template<typename M, typename V>
class ldiv_mat_monoblock_series_rep: public Series_matrix_rep {
protected:
  Matrix_series f;
  Matrix_series m;
  nat rf, cf, rm, cm;

  Monoblock_transformer blocker;
  bool h1_init;
  Series_matrix h0, h1;  
public:
  ldiv_mat_monoblock_series_rep (const Matrix_series& m2,
                                 const Matrix_series& f2):
    Series_matrix_rep (get_matrix_format (f2)), f(f2), m(m2), rf(rows(f2)),
    cf(cols(f2)), rm(rows(m2)), cm(cols(m2)), h1_init(false) {
    h0 = ldiv (m, f);
  }
  syntactic expression (const syntactic& z) const {
    return invert (flatten (m)) * flatten (f,z); }
  virtual void Increase_order (nat l) {
    Series_matrix_rep::Increase_order (l);
    if (l <= t)
      increase_order (h0, l);
    else {
      if (! h1_init) {
        typedef Monoblock_modular MM;
        typedef Lift_type (M) L;
        MM::set_modulus (binpow (L(* M::get_modulus ()), s));
        matrix<Monoblock_series> newm (Monoblock_series (), rm, cm);
        matrix<Monoblock_series> newf (Monoblock_series (), rf, cf);
        matrix<MM> init (MM (0), rm, cf);
        for (nat i=0; i<rm; i++) {
          for (nat j=0; j<cm; j++)
            newm(i,j)= to_monoblock (m(i,j), blocker);
          for (nat j=0; j<cf; j++) {
            newf(i,j)= to_monoblock (f(i,j), blocker);
            init(i,j)= as<L> (truncate (access (h0, i, j), s));
          }
        }
        matrix<Monoblock_series> newmh1=  asmatrix (ldiv (newm, newf, init));
        Matrix_series mh1 (Series(), rm, cf);
        for (nat i=0; i<rm; i++)
          for (nat j=0; j<cf; j++) {
            mh1(i,j)= from_monoblock (newmh1(i,j), blocker);
          }
        h1= asseries (mh1);
        h1_init= true;
      }
      increase_order (h1, l); } }
  Matrix next () {
    return this->n < t ? h0[this->n] : h1[this->n]; }
};

TMPL static inline Series_matrix
ser_lmul_quo_sc (const Matrix& c, const Series_matrix& f) {
  typedef implementation<series_matrix_divide,W> Ser;
  return Ser::ser_lmul_quo_sc (c, f); }

TMPL static inline Series_matrix
ser_lmul_rem_sc (const Matrix& c, const Series_matrix& f) {
  typedef implementation<series_matrix_divide,W> Ser;
  return Ser::ser_lmul_rem_sc (c, f); }

TMPL static inline Series_matrix
ser_ldiv_sc (const Matrix& c, const Matrix_series& f) {
  typedef implementation<series_matrix_divide,W> Ser;
  return Ser::ser_ldiv_sc (c, f);
}
  
TMPL static inline Series_matrix
ser_ldiv_sc (const Matrix& c, const Matrix_series& f, const Matrix& inv) {
  typedef implementation<series_matrix_divide,W> Ser;
  return Ser::ser_ldiv_sc (c, f, inv);
}

TMPL static inline Series_matrix
ser_ldiv (const Matrix_series& c, const Matrix_series& f) {
  typedef ldiv_mat_monoblock_series_rep<M,W> Div_rep;
  series<matrix<M>,W> s =
    series<matrix<M>,W> ((series_rep<matrix<M>,W>*) new Div_rep
                         (as<matrix<series<M,W> > > (c),
                          as<matrix<series<M,W> > > (f)));
  return as<Series_matrix> (s);
}
  
}; // implementation<series_matrix_divide,U,series_carry_monoblock<W> >

#undef Monoblock_modular
#undef Monoblock_transformer
#undef Monoblock_series
#undef TMPLB

/******************************************************************************
* System of polynomial SLP regular root lifting
******************************************************************************/

#if 0
TMPL
class slp_polynomial_system_regular_root_series_rep:
    public recursive_series_rep<Vector, V> {
public:
  //protected:
  vector<generic> out, vars;
  Vector y0;

  inline vector<Series>
  rec_add (const vector<Series>& p, const vector<Series>& q) {
    return p + q;
  }
  inline vector<Series>
  rec_add (const vector<Series>& ep, vector<vector<L> >& tp,
           const vector<Series>& eq, const vector<vector<L> >& tq) {
    tp= tp + tq;
    return ep + eq;
  }

  inline vector<Series>
  rec_minus (const vector<Series>& p) {
    return -p;
  }
  inline vector<Series>
  rec_minus (const vector<Series>& ep, vector<vector<L> >& tp) {
    tp= -tp;
    return -ep;
  }

  inline vector<Series>
  rec_minus (const vector<Series>& p, const vector<Series>& q) {
    return p - q;
  }
  inline vector<Series>
  rec_minus (const vector<Series>& ep, vector<vector<L> >& tp,
             const vector<Series>& eq, const vector<vector<L> >& tq) {
    tp= tp - tq;
    return ep - eq;
  }

  inline vector<Series>
  rec_prod (const vector<Series>& p, const vector<Series>& q,
            const vector<Series>& y) {
    nat l= N(p);
    ASSERT (l == N(q) && (l-2) == N(y), "Wrong size of arguments");
    const Series& ep=p[0], ap=p[1], eq=q[0], aq=q[1];
    vector<Series> bp = range (p, 2, l);
    vector<Series> bq = range (q, 2, l);
    //Beware of as_vector (rshiftz (y)) which call
    // as_vector (y, N(f[1]))
    vector<Series> ry  = as_vector (rshiftz (as_series (y)), N(y));
    vector<Series> lry = as_vector (lshiftz_series_vector
                                    (rshiftz (as_series (y)), l-2));
    Series tp = ap + dot (bp, lry);
    Series tq = aq + dot (bq, lry);
    Series er = (ep * eq + lshiftz (tp * rshiftz (eq, 2), 2) + 
                 lshiftz (rshiftz (ep, 2) * tq, 2) +
                 lshiftz (dot (ry, bp) * dot (bq, ry), 2));
    return (vec<Series> (er, ap * aq) << (ap * bq + bp * aq));
  }
  inline vector<Series>
  rec_prod (const vector<Series>& ep, vector<vector<L> >& tp,
            const vector<Series>& eq, const vector<vector<L> >& tq,
            const vector<Series>& rme, const vector<Series>& pr_rme) {
    nat r= N(rme);
    ASSERT ((N(ep) == r+1) && (N(eq) == r+1), "Wrong size of tau");
    const vector<L>& ap= tp[0], aq= tq[0];
    vector<Series> aps (Series (), r);
    vector<Series> bps (Series (), r);
    for (nat i = 0; i < r; i++) {
      aps[i]= integer_to_series_carry<M,V> (ap[i]);  
      bps[i]= integer_to_series_carry<M,V> (bp[i]);
    }
    vector<vector<L> > bp (vector<L> (), r);
    vector<vector<L> > bq (vector<L> (), r);
    vector<vector<Series> > bps (vector<Series> (), r);
    vector<vector<Series> > bqs (vector<Series> (), r);
    for (nat i = 0; i < r; i++) {
      bp[i]= range (tp, 1, r+1);
      bq[i]= range (tq, 1, r+1);
      for (nat j= 0; j < r; j++) {
        bps[i][j]= integer_to_series_carry<M,V> (bp[i][j]);
        bqs[i][j]= integer_to_series_carry<M,V> (bq[i][j]);
      }
    }
    vector<Series> c1 (Series (), r);
    vector<Series> c2 (Series (), r);
    for (nat i = 0; i < r; i++) {
      c1[i]= (dot (bp[i], lshiftz (rme)) * rshiftz (eq[i], 2)
              + dot (bq[i], lshiftz (rme)) * rshiftz (ep[i], 2));
      c2[i]= aps[i] * rshiftz (eq[i], 2) + bps[i] * rshiftz (ep[i], 2);
    }
    vector<Series> res (Series (), r);
    for (nat i = 0; i < r; i++) {
      tp[i]= vec<L> (ap[i] * aq[i]) << (ap[i] * bq[i] + bp[i] * aq[i]);
      res[i]= c1 + c2;
      for (nat j = 0; j < r; j++)
        for (nat k = j; k < r; k++)
          res[i] += ((bp[i][j] * bq[i][k] + bp[i][k] * bq[i][j])
                     * pr_rme [i*r+j-i*(i+1)/2]);
      res[i] = ep[i] * eq[i] + lshiftz (res[i], 2);
    }
    return res;
  }

  inline vector<Series>
  rec_square (const vector<Series>& p, const vector<Series>& y) {
    nat l= N(p);
    ASSERT ((l-2) == N(y), "Wrong size of arguments");
    Series ser_2 = integer_to_series_carry<M,V> (2);
    const Series& ep=p[0], ap=p[1];
    vector<Series> bp = range (p, 2, l);
    vector<Series> ry = as_vector (rshiftz (as_series (y)), N(y));
    vector<Series> lry = as_vector (lshiftz_series_vector
                                    (rshiftz (as_series (y)), l-2));
    Series tp = ap + dot (bp, lry);
    Series er = (square (ep) + ser_2 * lshiftz (tp * rshiftz (ep, 2), 2) +
                 lshiftz (square (dot (bp, ry)), 2));
    return (vec<Series> (er, square (ap)) << (ser_2 * ap * bp));
  }

  inline vector<Series>
  rec_lin (const Series& a, const vector<Series>& b) {
    return vec<Series> (Series (CF(a)), a) << b;
  }

  inline vector<Series>
  rec_cst (const Series& p, nat n) {
    return rec_lin (p, vector<Series> (Series (CF(p)), n));
  }


  static Series
  _ev_der (const generic& f, const vector<Series>& a,
           const vector<generic>& x, vector<Series>& grad) {
    grad = vector<Series> (Series (get_format1 (CF(a))), N(a));
    if (is<literal> (f))
      for (nat i = 0; i < N(a); i++)
        if (f == x[i]) {
          grad[i]= Series (M(1));
          return a[i];
        }
    if (is<Series> (f))
      return as<Series> (f);
    if (is<compound> (f)) {
      Series o;
      if (N(f) == 2) {
        if (f[0] == GEN_MINUS) {
          o = -_ev_der (f[1], a, x, grad);
          grad = -grad;
          return o;
        }
        if (f[0] == GEN_SQUARE) {
          o = _ev_der (f[1], a, x, grad);
          Series ser_2=  integer_to_series_carry<M,V> (2);
          grad = ser_2 * o * grad;
          return square (o);
        }
      }
      if (N(f) == 3) {
        vector<Series> grad2;
        if (f[0] == GEN_MINUS) {
          o = _ev_der (f[1], a, x, grad) - _ev_der (f[2], a, x, grad2);
          grad = grad - grad2;
          return o;
        }
        if (f[0] == GEN_PLUS) {
          o = _ev_der (f[1], a, x, grad) + _ev_der (f[2], a, x, grad2);
          grad = grad + grad2;
          return o;
        }
        if (f[0] == GEN_TIMES) {
          o =  _ev_der (f[1], a, x, grad);
          Series o2 = _ev_der (f[2], a, x, grad2);
          grad = o*grad2 + grad*o2;
          return o*o2;
        } } }
    ERROR ("bug, unexpected type with generic"); }

  static Series_vector
  _ev_der (const vector<generic>& f, const vector<Series>& a,
           const vector<generic>& x, matrix<Series, matrix_naive>& m) {
    m = matrix<Series, matrix_naive> (Series (get_format1 (CF(a))),
                                      N(f), N(a));
    vector<Series> v (Series (get_format1 (CF(a))), N(f));
    vector<Series> grad;
    for (nat i = 0; i < N(f); i++) {
      v[i] = _ev_der (f[i], a, x, grad);
      for (nat j = 0; j < N(a); j++)
        m (i,j) = grad [j];
    }
    return as_series (v); }
  
  vector<Series> 
  _eps (const generic& f, const vector<generic>& x, const Series_vector& y) {
    nat n = N(x);
    if (is<Series> (f))
      return rec_cst (as<Series> (f), n);
    if (is<literal> (f)) {
      for (nat i = 0; i < n; i++)
        if (f == x[i]) {
          vector<Series> tmp (Series (M(0) /*FIXME*/), n);
          tmp [i] = Series (M(1));
          return rec_lin (Series (y0[i]), tmp);
        }
    }
    if (is<compound> (f)) {
      if (N(f) == 2) {
        if (f[0] == GEN_MINUS)
          return rec_minus (_eps (f[1], x, y));
        if (f[0] == GEN_SQUARE)
          return rec_square (_eps (f[1], x, y), as_vector (y));
      }
      if (N(f) == 3) {
        if (f[0] == GEN_MINUS)
          return rec_minus (_eps (f[1], x, y), _eps (f[2], x, y));
        if (f[0] == GEN_PLUS)
          return rec_add (_eps (f[1], x, y), _eps (f[2], x, y));
        if (f[0] == GEN_TIMES)
          return rec_prod (_eps (f[1], x, y), _eps (f[2], x, y),
                           as_vector (y));
      }
    }
    ERROR ("bug, unexpected type with generic"); }

  vector<Series>
  _eps (const vector<generic>& e, const vector<generic>& x, 
        const Series_vector& y) {
    vector<Series> veps (Series (M(0) /*FIXME*/), N(e));
    for (nat i = 0; i < N(e); i++) veps [i] = _eps (e[i], x, y)[0];
    return veps;
  }

  static void
  increase_order_generic (generic f, nat l) {
    if (is<Series> (f))
      increase_order (as<Series> (f), l);
    if (is<compound> (f)) {
      if (N(f) == 2)
        increase_order_generic (f[1], l);
      if (N(f) == 3) {
        increase_order_generic (f[1], l);
        increase_order_generic (f[2], l);
      }
    }
  }

  static inline void
  increase_order_generic (vector<generic> f, nat l) {
    for (nat i = 0; i < N(f); i++)
      increase_order_generic (f[i], l);
  }
  
  
  public:
  slp_polynomial_system_regular_root_series_rep (const vector<generic>& out2,
                                                 const vector<generic>& vars2,
                                                 const Vector& y02):
    recursive_series_rep<Vector,V> (get_format (y02)),
    out (out2),  vars (vars2), y0 (y02) {}
  syntactic expression (const syntactic& z) const {
    return z; }
    //return flatten (f, z) / flatten (c); }
  virtual void Increase_order (nat l) {
    recursive_series_rep<Vector,V>::Increase_order (l); 
    increase_order_generic (out, l); }
  Series_vector initialize () {
    //double start = mmx_time ();
    //double inter;
    
    nat n= N(out);
    ASSERT (N(vars) == n, "numbers of equations and variables don't match");
    this->initial (0) = y0;
    matrix<Series, matrix_naive> jac;
    vector<Series> num = as_vector (_ev_der (out, as<vector<Series> > (y0),
                                             vars, jac));
    //mmout << "init: " << (int) (mmx_time() - start) << " + ";
    //inter = mmx_time ();
    vector<Series> eps = _eps (out, vars, this->me());
    //mmout << (int) (mmx_time () - inter) << " + ";
    //inter = mmx_time ();
    Series_vector res = ldiv (jac, (jac * as<vector<Series> > (y0)
                                    - num - eps));
    //mmout << (int) (mmx_time () - inter)  << " = ";
    //mmout << (int) (mmx_time () - start)  << " ms" << lf;
    return res;
  }
};

TMPL static inline vector<Series>
slp_polynomial_system_regular_root_series (const vector<generic>& f,
                                           const vector<generic>& x,
                                           const Vector& y0) {
  typedef slp_polynomial_system_regular_root_series_rep<M,V> Sys_root_rep;
  Series_vector s= (series_rep<Vector,V>*) new Sys_root_rep (f, x, y0);
  return as_vector (recursive(s));
}
#endif


#undef Series_vector_rep
#undef Series_matrix_rep
#undef Vector
#undef Series_vector
#undef Vector_series
#undef Series_matrix
#undef Matrix_series
#undef Matrix
#undef Series
#undef Series_rep
#undef C
#undef TMPL
} // namespace mmx
#endif // __MMX_SERIES_CARRY_LINEAR_ALGEBRA_HHP

---