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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : series_carry_naive.hpp
* DESCRIPTION: lazy p-adic numbers
* COPYRIGHT  : (C) 2011  Jeremy Berthomieu, Gregoire Lecerf 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_NAIVE_HPP
#define __MMX_SERIES_CARRY_NAIVE_HPP
#include <algebramix/p_expansion.hpp>
#include <algebramix/series_naive.hpp>
#include <algebramix/root_modular.hpp>
#include <algebramix/series_vector.hpp>
#include <basix/literal.hpp>
#include <basix/compound.hpp>
#include <basix/timer.hpp>

namespace mmx {
#define Series     series<M,V>
#define Series_rep series_rep<M,V>
#define Recursive_series_rep recursive_series_rep<M,V>
#define C       typename M::C
#define Modulus typename M::modulus
#define TMPL    template<typename M, typename V>
#define TMPLW   template<typename M, typename V, typename W>
#define Vector        vector<M,W>
#define Series_vector series<Vector, V>  

/******************************************************************************
* Variant helpers
******************************************************************************/

struct series_carry_naive {};

template<typename V,typename W>
struct implementation<V,W,series_carry_naive>:
  public implementation<V,W,series_naive> {};

template<typename M>
struct series_polynomial_helper<M,series_carry_naive> {
  typedef typename P_expansion_variant(M) PV;
};

template<typename M>
struct series_carry_variant_helper {
  typedef series_carry_naive SV;
};

#define Series_carry_variant(C) series_carry_variant_helper<C>::SV

/******************************************************************************
* Output
******************************************************************************/

template<typename U>
struct implementation<series_defaults,U,series_carry_naive> {
  typedef implementation<series_defaults,U,series_naive> Ser;

  template<typename S>
  class global_variables :
    public Ser::template global_variables<S> {

    static inline generic& dyn_name () {
      static generic name = "p";
      return name; }
    
  public:
    static inline void set_variable_name (const generic& x) {
      dyn_name () = x; }
    static inline generic get_variable_name () {
      return dyn_name (); }
  };
};

/******************************************************************************
* Abstract scalar operations on series
******************************************************************************/

template<typename U>
struct implementation<series_scalar_abstractions,U,series_carry_naive>:
  public implementation<series_defaults,U>
{

template<typename Op,typename M,typename V,typename X>
class binary_scalar_series_rep: public Series_rep {
protected:
  Series f;
  M x;
  C c;
public:
  inline binary_scalar_series_rep (const Series& f2, const X& x2):
    Series_rep (CF(f2)), f(f2), x (as<M> (x2)), c (0) {}
  syntactic expression (const syntactic& z) const {
    return Op::op (flatten (f, z), flatten (x)); }
  virtual void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (f, l); }
  virtual M next () {
    return Op::op_mod (f[this->n].rep, x.rep, M::get_modulus (), c); }
};

}; // implementation<series_scalar_abstractions,V,series_carry_naive>

/******************************************************************************
* Abstract operations on p_adic
******************************************************************************/

template<typename U>
struct implementation<series_abstractions,U,series_carry_naive>:
  public implementation<series_recursive_abstractions,U>
{

template<typename Op,typename M,typename V>
class unary_series_rep: public Series_rep {
protected:
  const Series f;
  C carry;
public:
  inline unary_series_rep (const Series& f2):
    Series_rep (CF(f2)), f(f2), carry (C(0)) {}
  syntactic expression (const syntactic& z) const {
    return Op::op (flatten (f, z)); }
  virtual void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (f, l); }
  virtual M next () {
    return Op::op_mod (f[this->n].rep, M::get_modulus (), carry); }
};

template<typename Op,typename M,typename V>
class binary_series_rep: public Series_rep {
protected:
  Series f, g;
  C carry;
public:
  inline binary_series_rep (const Series& f2, const Series& g2):
    Series_rep (CF(f2)), f(f2), g (g2), carry (0) {}
  syntactic expression (const syntactic& z) const {
    return Op::op (flatten (f, z), flatten (g, z)); }
  virtual void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (f, l);
    increase_order (g, l); }
  virtual M next () {
    return Op::op_mod (f[this->n].rep, g[this->n].rep,
                       M::get_modulus (), carry); }
};

}; // implementation<series_abstractions,U,series_carry_naive>:

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

template<typename U>
struct implementation<series_multiply,U,series_carry_naive>:
  public implementation<series_abstractions,U>
{
  typedef implementation<series_abstractions,U> Ser;

TMPL
class mul_series_rep: public Ser::template binary_series_rep<mul_op,M,V> {
  typedef Lift_type (M) L;
  L carry;
public:
  inline mul_series_rep (const Series& f, const Series& g):
    Ser::template binary_series_rep<mul_op,M,V> (f, g), carry (0) {}
  M next () {
    Modulus p= M::get_modulus ();
    L acc= carry;
    for (nat i= 0; i <= this->n; i++)
      acc += L ((this->f[i]).rep) * L((this->g[this->n-i]).rep);
    return rem (acc, * p, carry); }
};

TMPL
class truncate_mul_series_rep:
    public Ser::template binary_series_rep<mul_op,M,V> {
  nat nf, ng;
  typedef Lift_type (M) L;
  L carry;
public:
  inline truncate_mul_series_rep (const Series& f, const Series& g,
                                  nat nf2, nat ng2):
    Ser::template binary_series_rep<mul_op,M,V> (f, g), nf (nf2), ng (ng2),
    carry (0) {}
  M next () {
    nat i, n= this->n;
    Modulus p= M::get_modulus ();
    L acc= carry;
    for (i= n >= ng ? n - ng + 1 : 0; i < min (1 + n, nf); i++)
      acc += L ((this->f[i]).rep) * L((this->g[n-i]).rep);
    return rem (acc, * p, carry); }
};

TMPL static inline Series
ser_mul (const Series& f, const Series& g) {
  typedef mul_series_rep<M,V> Mul_rep;
  if (is_exact_zero (f) || is_exact_zero (g))
    return Series (CF(f));
  return (Series_rep*) new Mul_rep (f, g); }

TMPL static inline Series
ser_truncate_mul (const Series& f, const Series& g, nat nf, nat ng) {
  typedef truncate_mul_series_rep<M,V> Mul_rep;
  if (is_exact_zero (f) || is_exact_zero (g) || nf == 0 || ng == 0)
    return Series (CF(f));
  return (Series_rep*) new Mul_rep (f, g, nf, ng); }

}; // implementation<series_multiply,U,series_carry_naive>:

/******************************************************************************
* Multiplication with lift & reduction to series
******************************************************************************/

template<typename V>
struct series_carry_lift {};

template<typename U,typename V,typename W>
struct implementation<V,U,series_carry_lift<W> >:
  public implementation<V,U,W> {};

template<typename U,typename W>
struct implementation<series_multiply,U,series_carry_lift<W> >:
  public implementation<series_abstractions,U>
{
  typedef implementation<series_abstractions,U> Ser;

template<typename M,typename V>
class mul_series_rep: public Ser::template binary_series_rep<mul_op,M,V> {
  typedef C L;
  series<L> H;
  L c;
public:
  inline mul_series_rep (const Series& f, const Series& g):
    Ser::template binary_series_rep<mul_op,M,V> (f, g), c(0) {
    H= as<series<L> > (f) * as<series<L> > (g); }
  M next () {
    c += H[this->n];
    return M (as<C> (rem (c, * M::get_modulus (), c)), true); }
};

template<typename M,typename V>
class truncate_mul_series_rep:
    public Ser::template binary_series_rep<mul_op,M,V> {
  typedef C L;
  series<L> H;
  L c;
public:
  inline truncate_mul_series_rep (const Series& f, const Series& g,
                                   nat nf, nat ng):
    Ser::template binary_series_rep<mul_op,M,V> (f, g), c(0) {
    H= truncate_mul (as<series<L> > (f), as<series<L> > (g), nf, ng); }
  M next () {
    c += H[this->n];
    return M (as<C> (rem (c, * M::get_modulus (), c)), true); }
};

template<typename M,typename V> static inline Series
ser_mul (const Series& f, const Series& g) {
  typedef mul_series_rep<M,V> Mul_rep;
  if (is_exact_zero (f) || is_exact_zero (g))
    return Series (CF(f));
  return (Series_rep*) new Mul_rep (f, g); }

TMPL static inline Series
ser_truncate_mul (const Series& f, const Series& g, nat nf, nat ng) {
  typedef truncate_mul_series_rep<M,V> Mul_rep;
  if (is_exact_zero (f) || is_exact_zero (g) || nf == 0 || ng == 0)
    return Series (CF(f));
  return (Series_rep*) new Mul_rep (f, g, nf, ng); }

}; // implementation<series_multiply,V,series_carry_lift<U> >

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

template<typename U>
struct implementation<series_divide,U,series_carry_naive>:
  public implementation<series_recursive_abstractions,U>
{

template<typename M,typename V,typename X>
class carry_mul_sc_series_rep: public Series_rep {
protected:
  Series f;
  M x;
  C c;
public:
  inline carry_mul_sc_series_rep (const Series& f2, const X& x2):
    Series_rep (CF(f2)), f(f2), x (as<M> (x2)), c (0) {}
  syntactic expression (const syntactic& z) const {
    return apply (syntactic ("carry_mul"), flatten (f, z), flatten (x)); }
  virtual void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (f, l); }
  virtual M next () {
    if (this->n == 0) return M (0);
    return (quo (f[this->n-1].rep * x.rep, * M::get_modulus ())); }
};

template<typename M,typename V> static inline Series
carry_mul_series (const Series& f, const M& c) {
  return (Series_rep*) new carry_mul_sc_series_rep<M,V,M> (f, c);
}

template<typename M,typename V,typename X>
class rdiv_sc_series_rep: public Recursive_series_rep {
protected:
  Series f;
  M c;
public:
  rdiv_sc_series_rep (const Series& f2, const X& c2):
    Recursive_series_rep (CF(f2)), f(f2), c(as<M> (c2)) {}
  syntactic expression (const syntactic& z) const {
    return flatten (f, z) / flatten (c); }
  virtual void Increase_order (nat l) {
    Recursive_series_rep::Increase_order (l);
    increase_order (f, l); }
  Series initialize () {
    typedef typename Series_variant(C) SV;
    ASSERT (c != 0, "division by zero");
    M u= invert (c);
    this->initial (0)= u * f[0];
    Series tmp= f - carry_mul_series (this -> me (), c);
    return as<Series> (u * as<series<M,SV> > (tmp));
  }
};

TMPL static inline Series
ser_rdiv_sc (const Series& f, const M& c) {
  typedef rdiv_sc_series_rep<M,V,M> Div_sc_rep;
  if (is_exact_zero (f)) return Series (CF(f));
  return recursive (Series ((Series_rep*) new Div_sc_rep (f, c))); }

template<typename M,typename V>
class div_series_rep: public Recursive_series_rep {
protected:
  Series f, g;

public:
  div_series_rep (const Series& f2, const Series& g2):
    Recursive_series_rep (CF(f2)), f(f2), g(g2) {}
  syntactic expression (const syntactic& z) const {
    return flatten (f, z) / flatten (g, z); }
  virtual void Increase_order (nat l) {
    Recursive_series_rep::Increase_order (l);
    increase_order (f, l);
    increase_order (g, l); }
  Series initialize () {
    typedef typename Series_variant(C) SV;
    ASSERT (g[0] != 0, "division by zero");
    M u= invert (g[0]);
    this->initial (0)= u * f[0];
    Series tmp= f - lshiftz (this -> me () * rshiftz (g))
      - carry_mul_series (this -> me (), g[0]);
    return as<Series> (u * as<series<M,SV> > (tmp)); }
};

TMPL static inline Series
ser_div (const Series& f, const Series& g) {
  typedef div_series_rep<M,V> Div_rep;
  if (is_exact_zero (f)) return Series (CF(f));
  return recursive (Series ((Series_rep*) new Div_rep (f, g))); }

TMPL static inline Series
ser_quo (const Series& f, const M& c) {
  return ser_div_sc (f, c); }
  
TMPL static inline Series
ser_quo (const Series& f, const Series& g) {
  return ser_div (f, g); }

TMPL static inline Series
ser_rrem_sc (const Series& f, const M& c) {
  return Series (CF(f)); }

}; // implementation<series_divide,U,series_carry_naive>:

/******************************************************************************
* separable_root
******************************************************************************/

struct ser_carry_separable_root_op {
  // regular case when r is invertible in the residual field

  template<typename S, typename I> static inline S
  binpow_no_tangent (const S& me, const I& r) {
    typedef Scalar_type(S) M;
    // me^r - r * (me - me[0]) * me[0]^(r-1) - me[0]^r
    //M me0= me[0]; S me1= lshiftz (me,1);
    if (r <= 1) return S (0);
    if (r == 2) return lshiftz (square (rshiftz (me, 1)), 2);

    I h= r >> 1;
    S s_h (coefficients (as<default_p_expansion(M)> (integer (h))));
    M x= me[0];
    S x_h_1= binpow (S(x), h-1);
    S y= me - me[0];
    S z= rshiftz (binpow_no_tangent (me, h), 1);
    S t0= s_h * y;
    S t1= square (rshiftz (t0 * x_h_1, 1));
    S t2= x_h_1 * z * (t0 + x);
    S w= lshiftz (square (z) + t1, 2) + lshiftz (t2 + t2, 1);
    if ((r & 1) == 1) {
      I h2= h << 1;
      S s_h2 (coefficients (as<default_p_expansion(M)> (integer (h2))));
      w= lshiftz (me * rshiftz (w, 1), 1)
        + lshiftz (s_h2 * square (rshiftz (y * x_h_1, 1)) * x, 2);
    }
    return w; }

  static generic name () { return "separable_root"; }

  template<typename S> static inline S
  op (const S& x, nat r) { return separable_root (x, r); }

  static inline syntactic
  op (const syntactic& x, const syntactic& r) {
    return apply ("separable_root", x, r); }

  template<typename R, typename S> static inline void
  set_op (R& x, const S& y, nat r) { x= separable_root (y, r); }

  template<typename S, typename I> static inline S
  op_init (const S& x, nat r, const I& i) {
    return separable_root_init (x, r, i); }

  static inline nat nr_init () { return 1; }

  template<typename S> static inline S
  def (const S& me, const S& a, nat r) {
    typedef Scalar_type(S) M;
    if (r == 2) {
      return ((a - S (me[0]) * me[0] - binpow_no_tangent (me, r))
              / M (r)) / me[0];
    }
    S ser_r (coefficients (as<default_p_expansion(M)> (integer (r))));
    S me_0_r_1= binpow (S (me[0]), r - 1);
    S me_0_r= me[0] * me_0_r_1;
    return (a - me_0_r - binpow_no_tangent (me, r))
      / (ser_r * me_0_r_1); }

}; // struct ser_carry_separable_root_op

template<typename U>
struct implementation<series_separable_root,U,series_carry_naive> {

TMPL static inline Series
sep_root (const Series& a, nat r) {
  ASSERT (M(r) != 0, "wrong argument");
  if (is_exact_zero (a)) return Series (CF(a));
  return binary_scalar_recursive_series
    <ser_carry_separable_root_op> (a, r); }

TMPL static inline Series
sep_root (const Series& a, nat r, const M& init) {
  ASSERT (M(r) != 0, "wrong argument");
  if (is_exact_zero (a)) return Series (CF(a));
  return binary_scalar_recursive_series
    <ser_carry_separable_root_op> (a, r, init); }

}; // implementation<series_separable_root,V,series_carry_naive>


/******************************************************************************
* simple_root
******************************************************************************/
TMPL inline vector<Series>
rec_add (const vector<Series>& p, const vector<Series>& q) {
  return p + q;
}

TMPL inline vector<Series>
rec_minus (const vector<Series>& p) {
  return -p;
}

TMPL inline vector<Series>
rec_minus (const vector<Series>& p, const vector<Series>& q) {
  return p - q;
}

TMPL inline vector<Series>
rec_prod (const vector<Series>& p, const vector<Series>& q,
          const Series& y) {
  ASSERT (N(p) == N(q) && N(p) == 3, "Wrong size of arguments");
  const Series& ep=p[0], ap=p[1], eq=q[0], aq=q[1];
  const Series& bp=p[2], bq=q[2];
  Series ry = rshiftz (y); 
  Series tp = ap + bp * lshiftz (ry);
  Series tq = aq + bq * lshiftz (ry);
  //TODO stocker récursivement rshiftz (er, 2) plutôt 
  Series er = (ep * eq + lshiftz (tp * rshiftz (eq, 2), 2) + 
               lshiftz (rshiftz (ep, 2) * tq, 2) + 
               bp * bq * lshiftz (square (ry), 2));
  return vec<Series> (er, ap * aq, ap * bq + bp * aq);   
}

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

TMPL inline vector<Series>
rec_square (const vector<Series>& p, const Series& y) {
  ASSERT (N(p) == 3, "Wrong size of arguments");
  Series ser_2 (coefficients (as<default_p_expansion(M)> (2)));
  const Series& ep=p[0], ap=p[1], bp=p[2];
  Series ry = rshiftz (y); 
  Series tp = ap + bp * lshiftz (ry);
  Series er = (square (ep) + ser_2 * lshiftz (tp * rshiftz (ep, 2), 2) +
               lshiftz (square (bp * ry), 2));
  return vec<Series> (er, square (ap), ser_2 * ap * bp);   
}


TMPL 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 (coefficients (as<default_p_expansion(M)> (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));
}

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

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

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

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


TMPL
class root_series_rep: public Recursive_series_rep {
public:
  //protected:
  generic f, x;
  M y0;
  static Series _eval (const generic& e, const Series& a) {
    if (is<literal> (e))
      return a;
    if (is<Series> (e))
      return as<Series> (e);
    if (is<compound> (e)) {
      if (N(e) == 2) {
        if (e[0] == GEN_MINUS)
          return - _eval (e[1], a);
        if (e[0] == GEN_SQUARE)
          return square (_eval (e[1], a));
      }
      if (N(e) == 3) {
        if (e[0] == GEN_MINUS)
          return _eval (e[1], a) - _eval (e[2], a);
        if (e[0] == GEN_PLUS)
          return _eval (e[1], a) + _eval (e[2], a);
        if (e[0] == GEN_TIMES)
          return _eval (e[1], a) * _eval (e[2], a);
      }
    }
    ERROR ("bug, unexpected type with generic"); }

  static generic _derive (const generic& e, const generic& x) {
    if (e == x) return as<generic> (Series (1));
    if (is<Series> (e))
      return as<generic> (as<Series> (0));
    if (is<compound> (e)) {
      if (N(e) == 2) {
        if (e[0] == GEN_MINUS)
          return gen (GEN_MINUS, _derive (e[1], x));
        if (e[0] == GEN_SQUARE) {
          Series ser_2 (coefficients (as<default_p_expansion(M)> (2)));
          return gen (GEN_TIMES, as<generic> (ser_2),
                      gen (GEN_TIMES, _derive (e[1], x), e[1]));
        }
      }
      if (N(e) == 3) {
        if (e[0] == GEN_MINUS)
          return gen (GEN_MINUS, _derive (e[1], x), _derive (e[2], x));
        if (e[0] == GEN_PLUS)
          return gen (GEN_PLUS, _derive (e[1], x), _derive (e[2], x));
        if (e[0] == GEN_TIMES)
          return gen (GEN_PLUS, gen (GEN_TIMES, _derive (e[1], x), e[2]),
                      gen (GEN_TIMES, e[1], _derive (e[2], x)));
      }
    }
    ERROR ("bug, unexpected type with generic"); }

  vector<Series> _eps (const generic& e, const generic& x) {
    if (is<Series> (e))
      return rec_cst (as<Series> (e));
    if (e == x)
      return rec_lin (Series (y0), Series (M(1)));
    if (is<compound> (e)) {
      if (N(e) == 2) {
        if (e[0] == GEN_MINUS)
          return rec_minus (_eps (e[1], x));
        if (e[0] == GEN_SQUARE)
          return rec_square (_eps (e[1], x), this->me ());
      }
      if (N(e) == 3) {
        if (e[0] == GEN_MINUS)
          return rec_minus (_eps (e[1], x), _eps (e[2], x));
        if (e[0] == GEN_PLUS)
          return rec_add (_eps (e[1], x), _eps (e[2], x));
        if (e[0] == GEN_TIMES)
          return rec_prod (_eps (e[1], x), _eps (e[2], x), this->me ());
      }
    }
    ERROR ("bug, unexpected type with generic"); }

  static void
  increase_order_generic (generic ff, nat l) {
    if (is<Series> (ff))
      increase_order (as<Series> (ff), l);
    if (is<compound> (ff)) {
      if (N(ff) == 2)
        increase_order_generic (ff[1], l);
      if (N(ff) == 3) {
        increase_order_generic (ff[1], l);
        increase_order_generic (ff[2], l);
      }
    }
  }
  
  public:
  root_series_rep (const generic& f2, const generic& x2, const M& y02):
    Recursive_series_rep (get_format (y02)), f (f2), x (x2), 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::Increase_order (l);
    increase_order_generic (f, l); }
  Series initialize () {
    this->initial (0) = y0;
    generic der = _derive (f, x);
    Series evder = _eval (der, y0);
    Series num = (Series (y0) * evder - _eval(f, Series (y0)) 
                  - _eps (f, x)[0]) / evder;
    return num;
  }
};

TMPL static inline Series
root_series (const generic& f, const generic& x, const M& y0) {
  return recursive(Series ((Series_rep*) new root_series_rep<M,V> (f, x, y0)));
}



TMPLW
class system_root_series_rep: public recursive_series_rep<Vector, V> {
public:
  //protected:
  vector<generic> out, vars;
  Vector y0;

  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 (coefficients (as<default_p_expansion(M)> (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:
  system_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 = as_series (ser_ldiv_mat (jac, (jac * \
                                   as<vector<Series> > (y0) - num - eps)));
    //mmout << (int) (mmx_time () - inter)  << " = ";
    //mmout << (int) (mmx_time () - start)  << " ms" << lf;
    return res;
  }
};

TMPLW static inline vector<Series,W>
system_root_series (const vector<generic>& f, const vector<generic>& x,
                    const Vector& y0) {
  typedef system_root_series_rep<M,V,W> Sys_root_rep;
  Series_vector s= (series_rep<Vector,V>*) new Sys_root_rep (f, x, y0);
  return as_vector (recursive(s));
}

#undef Series
#undef Series_rep
#undef Recursive_series_rep
#undef C
#undef Modulus
#undef TMPL
#undef TMPLW
#undef Vector
#undef Series_vector

} // namespace mmx
#endif // __MMX_SERIES_CARRY_NAIVE_HPP