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  J. Berthomieu, R. Lebreton and G. 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_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/table.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>

/******************************************************************************
* 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;
public:
  inline carry_mul_sc_series_rep (const Series& f2, const X& x2):
    Series_rep (CF(f2)), f(f2), x (as<M> (x2)) {}
  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);
    C unused= 0;
    C ret= 0;
    mul_mod (unused, *f[this->n-1], *x, M::get_modulus (), ret);
    return M(ret); }
};

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_rquo_sc (const Series& f, const M& c) {
  return ser_rdiv_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>

/******************************************************************************
* sign aware unsigned_p_expansion
******************************************************************************/

TMPL inline Series 
integer_to_series_carry (integer a) {
  if (a == integer (0)) return Series ();
  if (a > integer (0))
    return Series (coefficients (as<default_p_expansion (M)> (a)));
  return -Series (coefficients (as<default_p_expansion (M)> (-a)));
}

/******************************************************************************
* Polynomial SLP regular root lifting
******************************************************************************/
struct series_slp_polynomial_regular_root {};

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

template<typename M, typename V, typename L>
class slp_polynomial_regular_root_series_rep: public Recursive_series_rep {
public:
  //protected:
  generic f, x;
  M y0;
  Series rme, sq_rme;           // square of rshiftz(this->me ())
  table<Series,generic> ht_ser;
  table<vector<L>,generic> ht_tau;
  table<L,generic> ht_eval;
  
  static inline Series
  rec_add (const Series& ep, vector<L>& tp,
           const Series& eq, const vector<L>& tq) {
    tp= tp + tq;
    return ep + eq;
  }
  
  static inline Series
  rec_minus (const Series& ep, vector<L>& tp) {
    tp= -tp;
    return -ep;
  }
 
  static inline Series
  rec_minus (const Series& ep, vector<L>& tp,
             const Series& eq, const vector<L>& tq) {
    tp= tp - tq;
    return ep - eq;
  }

  static Series
  rec_prod (const Series& ep, vector<L>& tp,
            const Series& eq, const vector<L>& tq,
            const Series& rme, const Series& sq_rme) {
    ASSERT (N(tp) == 2 && N (tq) == 2, "Wrong size of tau");
    const L& ap= tp[0], bp= tp[1];
    const L& aq= tq[0], bq= tq[1];
    Series c1= (integer_to_series_carry<M,V> (bp) * rshiftz (eq, 2) 
                + integer_to_series_carry<M,V> (bq) * rshiftz (ep, 2));
    Series c2= (integer_to_series_carry<M,V> (ap) * rshiftz (eq, 2) 
                + integer_to_series_carry<M,V> (aq) * rshiftz (ep, 2));
    tp= vec<L> (ap * aq, ap * bq + bp * aq);
    return (ep * eq + lshiftz (c1 * lshiftz (rme) + c2
                               + integer_to_series_carry<M,V> (bp * bq)
                               * sq_rme, 2));
  }
  
  static inline Series
  rec_square (const Series& ep, vector<L>& tp, 
              const Series& rme, const Series sq_rme) {
    ASSERT (N(tp) == 2, "Wrong size of tau");
    const L& ap= tp[0], bp= tp[1];
    Series taup = (integer_to_series_carry<M,V> (2 * ap) + 
                   integer_to_series_carry<M,V> (2 * bp) * lshiftz (rme));
    tp= vec<L> (square (ap), 2 * ap * bp);
    return (square (ep) + lshiftz (taup * rshiftz (ep, 2), 2) +
            lshiftz (integer_to_series_carry<M,V> (square (bp)) * sq_rme, 2));
  } 
  
  //Old function for roots of polynomial with series coefficients
  static Series
  _eval2 (const generic& e, const Series& a) {
    if (is<literal> (e))
      return a;
    if (is<L> (e))
      return integer_to_series_carry<M,V> (as<L> (e));
    if (is<compound> (e)) {
      if (N(e) == 2) {
        if (e[0] == GEN_MINUS)
          return - _eval2 (e[1], a);
        if (e[0] == GEN_SQUARE)
          return square (_eval2 (e[1], a));
      }
      if (N(e) == 3) {
        if (e[0] == GEN_MINUS)
          return _eval2 (e[1], a) - _eval2 (e[2], a);
        if (e[0] == GEN_PLUS)
          return _eval2 (e[1], a) + _eval2 (e[2], a);
        if (e[0] == GEN_TIMES)
          return _eval2 (e[1], a) * _eval2 (e[2], a);
      }
    }
    ERROR ("bug, unexpected type with generic"); }
  
  static L
  _eval (const generic& e, const L& a, table<L,generic>& ht) {
    if (is<literal> (e))
      return a;
    if (is<L> (e))
      return as<L> (e);
    if (contains (ht, e))
      return ht[e];
    L l;
    if (is<compound> (e)) {
      if (N(e) == 2) {
        if (e[0] == GEN_MINUS)
          l = - _eval (e[1], a, ht);
        if (e[0] == GEN_SQUARE)
          l = square (_eval (e[1], a, ht));
      }
      if (N(e) == 3) {
        if (e[0] == GEN_MINUS)
          l = _eval (e[1], a, ht) - _eval (e[2], a, ht);
        if (e[0] == GEN_PLUS)
          l = _eval (e[1], a, ht) + _eval (e[2], a, ht);
        if (e[0] == GEN_TIMES)
          l = _eval (e[1], a, ht) * _eval (e[2], a, ht);
      }
      ht[e] = l;
      return l; }
    ERROR ("bug, unexpected type with generic"); }

  static generic
  _derive (const generic& e, const generic& x, 
           table<generic, generic>& ht_derive) {
    if (e == x) return as<generic> (L (1));
    if (is<L> (e))
      return as<generic> (as<L> (0));
    if (contains (ht_derive, e)) {
      return ht_derive [e];
    }
    generic der;
    if (is<compound> (e)) {
      if (N(e) == 2) {
        if (e[0] == GEN_MINUS)
          der = gen (GEN_MINUS, _derive (e[1], x, ht_derive));
        else if (e[0] == GEN_SQUARE) {
          der = gen (GEN_TIMES, as<generic> (L (2)),
                      gen (GEN_TIMES, _derive (e[1], x, ht_derive), e[1]));
        }
      }
      if (N(e) == 3) {
        if (e[0] == GEN_MINUS)
          der = gen (GEN_MINUS, _derive (e[1], x, ht_derive), 
                     _derive (e[2], x, ht_derive));
        else if (e[0] == GEN_PLUS)
          der = gen (GEN_PLUS, _derive (e[1], x, ht_derive), 
                     _derive (e[2], x, ht_derive));
        else if (e[0] == GEN_TIMES)
          der = gen (GEN_PLUS, 
                     gen (GEN_TIMES, _derive (e[1], x, ht_derive), e[2]),
                     gen (GEN_TIMES, e[1], _derive (e[2], x, ht_derive)));
      }
      ht_derive [e] = der;
      return der;
    }
    ERROR ("bug, unexpected type with generic"); }

  Series
  _eps (const generic& e, const generic& x, vector<L>& tau) {
    if (contains (ht_ser, e)) {
      tau = ht_tau[e];
      return ht_ser[e];
    }
    Series ser;
    if (is<L> (e)) {
      tau= vec<L> (as<L> (e), L (0));
      ser = Series ();
      ht_tau[e] = tau;
      ht_ser[e] = ser;
      return ser;
    }
    if (e == x) {
      tau= vec<L> (L (* y0), L(1));
      ser = Series ();
      ht_tau[e] = tau;
      ht_ser[e] = ser;
      return ser;
    }
    if (is<compound> (e)) {
      if (N(e) == 2) {
        Series s= _eps (e[1], x, tau);
        if (e[0] == GEN_MINUS)
          ser = rec_minus (s, tau);
        if (e[0] == GEN_SQUARE)
          ser = rec_square (s, tau, rme, sq_rme);
      }
      if (N(e) == 3) {
        vector<L> tau2 (L (0), 2);
        Series s1= _eps (e[1], x, tau), s2= _eps (e[2], x, tau2);
        if (e[0] == GEN_MINUS)
          ser = rec_minus (s1, tau, s2, tau2);
        if (e[0] == GEN_PLUS)
          ser = rec_add (s1, tau, s2, tau2);
        if (e[0] == GEN_TIMES)
          ser = rec_prod (s1, tau, s2, tau2, rme, sq_rme);
      }
      ht_tau[e] = tau;
      ht_ser[e] = ser;
      return ser; }
    ERROR ("bug, unexpected type with generic"); }

  /*
  //Old function for benchmark purposes
  static Series
  _eps_bench (const generic& e, const generic& x, vector<L>& tau, 
              const Series& y) {
    if (is<L> (e)) {
      tau= vec<L> (as<L> (e), L (0));
      return Series ();
    }
    if (e == x) {
      tau= vec<L> (L (* y[0]), L(1));
      return Series ();
    }
    if (is<compound> (e)) {
      if (N(e) == 2) {
        Series s= _eps_bench (e[1], x, tau, y);
        if (e[0] == GEN_MINUS)
          return rec_minus (s, tau);
        if (e[0] == GEN_SQUARE)
          return rec_square (s, tau, rshiftz (y), square (rshiftz (y)));
      }
      if (N(e) == 3) {
        vector<L> tau2 (L (0), 2);
        Series s1= _eps_bench (e[1], x, tau , y);
        Series s2= _eps_bench (e[2], x, tau2, y);
        if (e[0] == GEN_MINUS)
          return rec_minus (s1, tau, s2, tau2);
        if (e[0] == GEN_PLUS)
          return rec_add (s1, tau, s2, tau2);
        if (e[0] == GEN_TIMES)
          return rec_prod (s1, tau, s2, tau2, rshiftz (y), 
                           square (rshiftz (y)));
      }  }
    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:
  slp_polynomial_regular_root_series_rep (const generic& f2, const generic& x2,
                                          const M& y02):
    Recursive_series_rep (get_format (y02)), f (f2), x (x2), y0 (y02) {
    rme= rshiftz (this->me ());
    sq_rme= square (rme);
  }
  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 () {
    ht_ser = table<Series,generic> ();
    ht_tau = table<vector<L>,generic> ();
    table<L,generic> ht_eval = table<L,generic> ();
    L Ly0 (* y0);
    this->initial (0) = y0;
    table<generic, generic> ht_derive = 
      table<generic, generic> ();
    generic der = _derive (f, x, ht_derive);
    Series evder= integer_to_series_carry<M,V> ( _eval (der, Ly0, ht_eval));
    ht_eval = table<L,generic> ();
    L deriv = _eval (der, Ly0, ht_eval);
    ht_eval = table<L,generic> ();
    Series cst= integer_to_series_carry<M,V> (Ly0 * deriv 
                                                  - _eval (f, Ly0, ht_eval));
    //Beware :
    // The integer behind deriv and cst can grow exponentially in the 
    // multiplicative complexity of f. By computing all their p_adic 
    // decomposition at once, we have an exponential initialization time.
    // If we had Series instead of integer, this problem would not exist.
    vector<L> tau (L (0), 2);
    return (cst - _eps (f, x, tau)) / evder;
  }
};

template<typename M, typename V, typename L> static inline Series
slp_pol_root (const generic& f, const generic& x, const M& y0) {
  return recursive (Series ((Series_rep*)
                            new slp_polynomial_regular_root_series_rep<M,V,L>
                            (f, x, y0)));
}

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

template<typename M, typename V, typename L> static inline Series
slp_polynomial_regular_root (const generic& f, const generic& x, const M& y0) {
  typedef implementation<series_slp_polynomial_regular_root,V> Ser;
  return Ser::template slp_pol_root<M,V,L> (f, x, y0);
}


/******************************************************************************
* Polynomial regular root
******************************************************************************/

struct ser_carry_polynomial_regular_root_op {
  // regular case when the derivative does not vanish
  static generic name () { return "polynomial_regular_root"; }

  template<typename M> static inline M
  op (const polynomial<Lift_type (M)>& P) {
    typedef Lift_type (M) L;
    L p (* M::get_modulus());
    vector<M> cP (M(0), N(P));
    for (nat i= 0; i < N(P); i++)
      cP[i]= M (rem (P[i], p));
    return polynomial_modular_root (polynomial<M> (cP)); }
  
  template<typename M> static inline M
  op (const polynomial<M>& P, const M& init) {
    return init; }

  template<typename M, typename V, typename L> static inline Series
  op (const polynomial<L>& P) {
    return polynomial_regular_root<M,V,L> (P); }
  
  template<typename M, typename V, typename L> static inline Series
  op (const polynomial<L>& P, const M& init) {
    return polynomial_regular_root<M,V,L> (P,init); }
  
  static inline syntactic
  op (const syntactic& P) {
    return apply ("polynomial_regular_root", P); }
  
  template<typename R, typename L> static inline void
  set_op (R& x, const polynomial<L>& P) {
    x= polynomial_regular_root (P); }
  
  template<typename M, typename V, typename L> static inline Series
  op_init (const polynomial<L>& P, const M& i) {
    return polynomial_regular_root_init<M,V,L> (P, i); }

  static inline nat nr_init () { return 1; }

  template<typename M, typename V, typename L> static inline Series
  def (const Series& me, const polynomial<L>& P) {
    L y0 (* me[0]);
    polynomial<L> Q;
    polynomial<L> sqr (vec<L> (square (y0), - 2 * y0, 1)); // (X-y_O)²
    polynomial<L> R= rem (P, sqr, Q);
    Series Py0= integer_to_series_carry<M, V> (R[0]);
    Series dPy0= integer_to_series_carry<M, V> (R[1]);
    nat n = N(Q);
    vector<Series> vs (Series(), n);
    for (nat i= 0; i<n; i++)
      vs[i] = integer_to_series_carry<M, V> (Q[i]);
    polynomial<Series> sQ (vs);
    return ((Py0 - me[0] * dPy0 + lshiftz (square (rshiftz (me))
                                           * eval (sQ, me), 2)) / (-dPy0));
  }
}; // struct ser_carry_polynomial_regular_root_op

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

template<typename M, typename V, typename L> static inline Series
pol_root (const polynomial<L>& P) {
  ASSERT (N(P) > 1, "Polynomial must be non constant");
  return unary_polynomial_recursive_series
    <ser_carry_polynomial_regular_root_op,M,V,L> (P); }

template<typename M, typename V, typename L> static inline Series
pol_root (const polynomial<L>& P, const M& init) {
  ASSERT (N(P) > 1, "Polynomial must be non constant");
  return unary_polynomial_recursive_series
    <ser_carry_polynomial_regular_root_op,M,V,L> (P, init); }
}; // implementation<series_polynomial_regular_root,V,series_carry_naive>


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

} // namespace mmx
#endif // __MMX_SERIES_CARRY_NAIVE_HPP

---