include/algebramix/series.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : series.hpp
* DESCRIPTION: Univariate power series
* COPYRIGHT  : (C) 2004  Joris van der Hoeven
*******************************************************************************
* 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_HPP
#define __MMX_SERIES_HPP
#include <algebramix/series_naive.hpp>

namespace mmx {
#define Series_variant(C) series_variant_helper<C>::SV
#define TMPL template<typename C, typename V>
#define TMPLK template<typename C, typename V, typename K>
#define Format format<C>
#define Series series<C,V>
#define Series_rep series_rep<C,V>
#define Recursive_series_rep recursive_series_rep<C,V>
#define Polynomial polynomial<C, typename series_polynomial_helper<C,V>::PV>
#define PolynomialV polynomial<C,V>
#define Truncated_series Polynomial
#define Completed_series Series
#define Truncated_polynomial Polynomial
#define Completed_polynomial Series
TMPL syntactic flatten (const Series& f, const syntactic& z);
TMPL Series lshiftz (const Series& f, const int& shift= 1);
TMPL Series rshiftz (const Series& f, const int& shift= 1);

extern xnat mmx_bit_precision;
  
/******************************************************************************
* Dense series
******************************************************************************/

TMPL
class series_rep REP_STRUCT_1(C) {
public:   // should be protected
  C*  a;  // coefficients
  nat n;  // number of computed coefficients
  nat l;  // number of allocated coefficients
          // a[n],...,a[l-1] may be used by relaxed computations

  inline series_rep (const Format& fm):
    Format (fm), a (mmx_new<C> (0)), n (0), l (0) {}
  inline virtual ~series_rep () { mmx_delete<C> (a, l); }
  inline C zero () { return promote (0, this->tfm ()); }
  inline C one () { return promote (1, this->tfm ()); }
  inline Series me () const;
  virtual C next () = 0;
  virtual syntactic expression (const syntactic& z) const = 0;

public:
  virtual void Set_order (nat l2);
  virtual void Increase_order (nat l2=0);
  virtual inline bool test_exact_zero () const { return false; }
  friend class Series;
};

TMPL
class series {
INDIRECT_PROTO_2 (series, series_rep, C, V)
  typedef implementation<series_defaults,V> Ser;
  typedef typename Ser::template global_variables<Series> S;

public:
  static inline generic get_variable_name () {
    return S::get_variable_name (); }
  static inline void set_variable_name (const generic& x) {
    S::set_variable_name (x); }

  static inline nat get_output_order () {
    return S::get_output_order (); }
  static inline void set_output_order (const nat& x) {
    S::set_output_order (x); }

  static inline nat get_cancel_order () {
    return S::get_cancel_order (); }
  static inline void set_cancel_order (const nat& x) {
    S::set_cancel_order (x); }

  static inline bool get_formula_output () {
    return S::get_formula_output (); }
  static inline void set_formula_output (const bool& x) {
    S::set_formula_output (x); }

public:
  series ();
  series (const Format& fm);
  series (const C& c);
  template<typename T> series (const T& c);
  template<typename T> series (const T& c, const Format& fm);
  template<typename T> series (const T& c, nat deg);
  template<typename T> series (const polynomial<T>& P);
  template<typename T> series (const polynomial<T>& P, const Format& fm);
  template<typename T, typename W>
  series (const series<T,W>& f);
  template<typename T, typename W>
  series (const series<T,W>& f, const Format& fm);
  series (const vector<C>& coeffs);
  series (const iterator<C>& it, const string& name= "explicit");
  series (C (*coeffs) (nat), const string& name= "explicit");
  const C& operator [] (nat n) const;
  const C* operator () (nat start, nat end) const;
};
INDIRECT_IMPL_2 (series, series_rep, typename C, C, typename V, V)
DEFINE_UNARY_FORMAT_2 (series)

STYPE_TO_TYPE(TMPL,scalar_type,Series,C);

TMPL inline Format CF (const Series& f) { return f->tfm (); }

TMPL inline void set_variable_name (const Series&, const generic& x) {
  return Series::set_variable_name (x); }
TMPL inline void set_output_order (const Series&, const nat& n) {
  return Series::set_output_order (n); }
TMPL inline void set_cancel_order (const Series&, const nat& n) {
  return Series::set_cancel_order (n); }
TMPL inline void set_formula_output (const Series&, const bool& b) {
  return Series::set_formula_output (b); }

/******************************************************************************
* Recursive series
******************************************************************************/

TMPL
class recursive_series_rep: public Series_rep {
public:
  Series eq;
public:
  inline recursive_series_rep (const Format& fm):
    Series_rep (fm) {}
  virtual Series initialize () = 0;
  C& initial (nat n2) {
    if (n2>=this->n) { this->n = n2+1; this -> Set_order (this->n); }
    return this->a[n2]; }
  virtual void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (eq, l); }
  inline C next () { return eq[this->n]; }
};

TMPL
class recursive_container_series_rep: public Series_rep {
  Series f;
public:
  recursive_container_series_rep (const Series& f2):
    Series_rep (CF(f2)), f (f2) {
    Recursive_series_rep* rep= (Recursive_series_rep*) f.operator -> ();
    rep->eq= rep->initialize (); }
  ~recursive_container_series_rep () {
    Recursive_series_rep* rep= (Recursive_series_rep*) f.operator -> ();
    rep->eq= Series (); }
  syntactic expression (const syntactic& z) const { return flatten (f, z); }
  virtual void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (f, l); }
  C next () { return f[this->n]; }
};

TMPL Series
recursive (const Series& f) {
  return (Series_rep*) new recursive_container_series_rep<C,V> (f);
}

/******************************************************************************
* Routines on dense series
******************************************************************************/

TMPL inline generic var (const Series& f) {
  (void) f; return Series::get_variable_name ();
}

TMPL inline Series
Series_rep::me () const {
  return Series (this, true);
}

TMPL void
Series_rep::Set_order (nat l2) {
  typedef typename series_polynomial_helper<C,V>::PV PV;
  typedef implementation<polynomial_linear,PV> Pol;
  if (l2 <= l) return;
  C* b= mmx_new<C> (l2);
  Pol::copy (b, a, l);
  Pol::clear (b+l, l2-l);
  mmx_delete<C> (a, l);
  a= b;
  l= l2;
}

TMPL void
Series_rep::Increase_order (nat l2) {
  l2= max (l2, l << 1);
  Set_order (max ((nat) 1, l2));
}

TMPL void
increase_order (const Series& f, nat l) {
  if (f->l < l) inside (f) -> Increase_order (l);
}

TMPL const C&
Series::operator [] (nat n) const {
  if (n <  rep->n) return rep->a[n];
  if (n >= rep->l) rep->Increase_order (n+1);
  while (rep->n <= n) {
    rep->a[rep->n]= rep->next ();
    //mmout << "{" << rep->n << ";" << this << "}";
    rep->n++;
  }
  // mmout << (*this) << " [" << n << "] -> " << rep->a[n] << "\n";
  return rep->a[n];
}

TMPL const C*
Series::operator () (nat start, nat end) const {
  if (end <= rep->n) return rep->a + start;
  if (end >= rep->l) rep->Increase_order (end);
  while (rep->n < end) {
    rep->a[rep->n]= rep->next ();
    //mmout << "{" << rep->n << ";" << this << "}";
    rep->n++;
  }
  return rep->a + start;
}

TMPL Polynomial
truncate (const Series& f, nat n) {
  typedef typename series_polynomial_helper<C,V>::PV PV;
  nat l= aligned_size<C,PV> (n);
  C* coeffs= mmx_formatted_new<C> (l, CF(f));
  if (n>0) (void) f[n-1];
  for (nat i=0; i<n; i++) coeffs[i]= f[i];
  return Polynomial (coeffs, n, l, CF(f));
}

TMPL Polynomial
range (const Series& f, nat start, nat end) {
  typedef typename series_polynomial_helper<C,V>::PV PV;
  nat n= (end >= start? end - start: 0);
  nat l= aligned_size<C,PV> (n);
  C* coeffs= mmx_formatted_new<C> (l, CF(f));
  if (end>start) (void) f[end-1];
  for (nat i=0; i<n; i++) coeffs[i]= f[i + start];
  return Polynomial (coeffs, n, l, CF(f));
}

/******************************************************************************
* Series iterator
******************************************************************************/

TMPL
class series_iterator_rep: public iterator_rep<C> {
  Series f;  // the series to traverse
  nat i;     // current position
protected:
  bool is_busy () { return true; }
  void advance () { i++; }
  C current () { return f[i]; }
public:
  series_iterator_rep (const Series& f2):
    iterator_rep<C> (CF(f2)), f(f2), i(0) {}
};

TMPL iterator<C>
iterate (const Series& f) {
  return iterator<C> (new series_iterator_rep<C,V> (f));
}

/******************************************************************************
* Hashing and equality testing
******************************************************************************/

HARD_TO_EXACT_IDENTITY_SUGAR(TMPL,Series)

template<typename Op, typename C, typename V> nat
unary_hash (const Series& s) {
  register nat i, h= 7531;
  for (i=0; i< Series::get_cancel_order (); i++)
    h= (h<<1) ^ (h<<5) ^ (h>>27) ^ Op::op (s[i]);
  return h;
}

TMPL inline nat
hash (const Series& s) {
  return unary_hash<hash_op> (s);
}

template<typename Op,typename C,typename V> inline bool
binary_test (const Series& f1, const Series& f2) {
  for (nat n=0; n< Series::get_cancel_order (); n++)
    if (Op::not_op (f1[n], f2[n])) return false;
  return true;
}

TMPL bool operator == (const Series& f1, const Series& f2) {
  return binary_test<equal_op> (f1, f2); }
TMPL bool operator != (const Series& f1, const Series& f2) {
  return !binary_test<equal_op> (f1, f2); }

TMPL bool is_exact_zero (const Series& f) {
  return f->test_exact_zero (); }

/******************************************************************************
* Printing series
******************************************************************************/

TMPL syntactic
flatten (const Series& f, const syntactic& z) {
  if (Series::get_formula_output ()) return f->expression (z);
  syntactic s= 0;
  nat order= Series::get_output_order ();
  for (nat i=0; i<order; i++) {
    (void) f[order-1];
    C coeff= f[i];
    s= s + flatten (coeff) * pow (z, i);
  }
  s= s + apply ("O", pow (z, syntactic (order)));
  return s;
}

TMPL syntactic
flatten (const Series& f) {
  return flatten (f, as_syntactic (var (f)));
}

/******************************************************************************
* Comparisons
******************************************************************************/

TMPL int
val (const Series& f) {
  for (nat n=0; n< Series::get_cancel_order (); n++)
    if (f[n] != 0) return n;
  return (int) (((nat) (-1)) >> 1);
}

TMPL int
sign (const Series& f) {
  for (nat n=0; n< Series::get_cancel_order (); n++) {
    int sgn= sign (f[n]);
    if (sgn != 0) return sgn;
  }
  return 0;
}

TMPL int
compare (const Series& f1, const Series& f2) {
  for (nat n=0; n< Series::get_cancel_order (); n++) {
    int sgn= compare (f1[n], f2[n]);
    if (sgn != 0) return sgn;
  }
  return 0;
}

EQUAL_INT_SUGAR(TMPL,Series)
COMPARE_SUGAR(TMPL,Series)
COMPARE_INT_SUGAR(TMPL,Series)
EQUAL_SCALAR_SUGAR(TMPL,Series,C)
COMPARE_SCALAR_SUGAR(TMPL,Series,C)
EQUAL_SCALAR_SUGAR_BIS(TMPL,Series,C)
COMPARE_SCALAR_SUGAR_BIS(TMPL,Series,C)

/******************************************************************************
* Explicit series
******************************************************************************/

TMPL
class zero_series_rep: public Series_rep {
public:
  zero_series_rep (const Format& fm):
    Series_rep (fm) {}
  syntactic expression (const syntactic&) const {
    return flatten (0); }
  bool test_exact_zero () const {
    return true; }
  C next () { return this->zero (); }
};

TMPL Series::series () {
  rep= new zero_series_rep<C,V> (format<C> (no_format ())); }
TMPL Series::series (const Format& fm) {
  rep= new zero_series_rep<C,V> (fm); }

TMPL
class scalar_series_rep: public Series_rep {
  C c;
public:
  scalar_series_rep (const C& c2):
    Series_rep (get_format (c2)), c(c2) {}
  syntactic expression (const syntactic& z) const {
    (void) z; return flatten (c); }
  bool test_exact_zero () const {
    return is_exact_zero (c); }
  C next () { return this->n == 0 ? c : this->zero (); }
};

TMPL Series::series (const C& c) {
  rep= new scalar_series_rep<C,V> (c); }
TMPL template<typename T> Series::series (const T& c) {
  rep= new scalar_series_rep<C,V> (as<C> (c)); }
TMPL template<typename T> Series::series (const T& c, const Format& fm) {
  rep= new scalar_series_rep<C,V> (promote (c, fm)); }

TMPL
class polynomial_series_rep: public Series_rep {
  Polynomial P;
public:
  polynomial_series_rep (const polynomial<C>& P2):
    Series_rep (CF(P2)), P (as<Polynomial> (P2)) {}
  polynomial_series_rep (const polynomial<C>& P2, const Format& fm):
    Series_rep (fm), P (as<Polynomial> (P2)) {}
  syntactic expression (const syntactic& z) const {
    return flatten (P, z); }
  bool test_exact_zero () const {
    return is_exact_zero (P); }
  C next () {
    return this->n < N(P)? P[this->n]: this->zero (); }
};

TMPL template<typename T>
Series::series (const T& c, nat deg) {
  rep= new polynomial_series_rep<C,V> (polynomial<C> (as<C> (c), deg));
}

TMPL template<typename T>
Series::series (const polynomial<T>& P) {
  rep= new polynomial_series_rep<C,V> (polynomial<C> (P));
}

TMPL template<typename T>
Series::series (const polynomial<T>& P, const Format& fm) {
  rep= new polynomial_series_rep<C,V> (polynomial<C> (P), fm);
}

TMPL
Series::series (const vector<C>& coeffs) {
  rep= new polynomial_series_rep<C,V> (polynomial<C> (coeffs));
}

TMPL
class iterator_series_rep: public Series_rep {
  iterator<C> it;
  string name;
public:
  iterator_series_rep (const iterator<C>& it2, const string& name2):
    Series_rep (CF(it)), it (it2), name (name2) {}
  syntactic expression (const syntactic& z) const {
    return apply (syntactic (name), z); }
  C next () {
    if (busy (it)) {
      C c= *it; ++it; return c; }
    else return this->zero (); }
};

TMPL
Series::series (const iterator<C>& it, const string& name) {
  rep= new iterator_series_rep<C,V> (it, name);
}

/******************************************************************************
* Casting
******************************************************************************/

template<typename C, typename V>
class fast_series_rep: public series_rep<Fast_type(C) > {
  Series f;
public:
  fast_series_rep (const Series& f2):
    series_rep<Fast_type(C) > (fast (CF(f2))), f (f2) {}
  syntactic expression (const syntactic& z) const {
    return apply (GEN_FAST, flatten (f, z)); }
  void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (f, l); }
  Fast_type(C) next () { return fast (f[this->n]); }
};

template<typename C, typename V>
class slow_series_rep: public Series_rep {
  series<Fast_type(C) > f;
public:
  slow_series_rep (const series<Fast_type(C) >& f2):
    Series_rep (slow<C> (CF(f2))), f (f2) {}
  syntactic expression (const syntactic& z) const {
    return apply (GEN_SLOW, flatten (f, z)); }
  void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (f, l); }
  C next () { return slow<C> (f[this->n]); }
};

TMPL
struct fast_helper<Series > {
  typedef series<Fast_type(C) > fast_type;
  static inline fast_type dd (const Series& f) {
    return (series_rep<Fast_type(C) >*) new fast_series_rep<C,V> (f); }
  static inline Series uu (const fast_type& f) {
    return (Series_rep*) new slow_series_rep<C,V> (f); }
};

template<typename C, typename V, typename K, typename W>
class cast_series_rep: public Series_rep {
  series<K,W> f;
public:
  cast_series_rep (const series<K,W>& f2):
    Series_rep (as<Format> (CF(f2))), f (f2) {}
  cast_series_rep (const series<K,W>& f2, const Format& fm):
    Series_rep (fm), f (f2) {}
  syntactic expression (const syntactic& z) const { return flatten (f, z); }
  void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (f, l); }
  C next () { return as<C> (f[this->n]); }
};

TMPL template<typename T, typename W>
Series::series (const series<T,W>& f) {
  rep= new cast_series_rep<C,V,T,W> (f);
}

TMPL template<typename T, typename W>
Series::series (const series<T,W>& f, const Format& fm) {
  rep= new cast_series_rep<C,V,T,W> (f, fm);
}

template<typename T, typename F, typename TV, typename FV>
struct as_helper<series<T,TV>,series<F,FV> > {
  static inline series<T,TV>
  cv (const series<F,FV>& f) {
    return (series_rep<T,TV>*) new cast_series_rep<T,TV,F,FV> (f); }
};

template<typename T, typename F, typename TV, typename FV> inline void
set_as (series<T,TV>& r, const series<F,FV>& f) {
  r= series<T,TV> (f, CF(r));
}

template<typename C, typename V, typename T> inline void
set_as (Series& r, const T& x) {
  r= Series (x, CF(r));
}

/******************************************************************************
* Scalar abstractions
******************************************************************************/

template<typename Op,typename C,typename V,typename X> inline Series
binary_scalar_series (const Series& f, const X& x) {
  typedef implementation<series_scalar_abstractions,V> Ser;
  typedef typename Ser::template binary_scalar_series_rep<Op,C,V,X>
    Binary_rep;
  return (Series_rep*) new Binary_rep (f, x);
}

template<typename Op,typename C,typename V,typename X, typename Y>
inline Series
ternary_scalar_series (const Series& f, const X& x, const Y& y) {
  typedef implementation<series_scalar_abstractions,V> Ser;
  typedef typename Ser::template ternary_scalar_series_rep<Op,C,V,X,Y>
    Ternary_rep;
  return (Series_rep*) new Ternary_rep (f, x, y);
}

/******************************************************************************
* Map abstractions
******************************************************************************/

template<typename Op, typename C, typename V, typename S, typename SV> Series
unary_map_as (const series<S,SV>& f) {
  typedef implementation<series_map_as_abstractions,V> Ser;
  typedef typename Ser::
    template unary_map_as_series_rep<Op,C,V,S,SV> Map_as_rep;
  return (Series_rep*) new Map_as_rep (f);
}

/******************************************************************************
* Dynamic mappers
******************************************************************************/

template<typename C, typename V, typename S, typename SV, typename Fun>
class map_series_rep: public Series_rep {
protected:
  Fun fun;
  series<S,SV> f;
public:
  inline map_series_rep (const Fun& fun2,
                         const series<S,SV>& f2, const Format& fm):
    Series_rep (fm), fun (fun2), f (f2) {}
  syntactic expression (const syntactic& z) const {
    return syn (flatten (fun), flatten (f, z)); }
  virtual void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (f, l); }
  virtual C next () {
    return fun (f[this->n]); }
};

template<typename S1, typename D, typename Fun> series<D>
map (const Fun& fun, const series<S1>& f, const format<D>& fm) {
  typedef map_series_rep<D ,typename Series_variant(D ),
                         S1,typename Series_variant(S1),Fun> Mapper;
  return (series_rep<D>*) new Mapper (fun, f, fm);
}

/******************************************************************************
* Recursive abstractions
******************************************************************************/

template<typename Op,typename C> inline series<C>
nullary_recursive_series (const C& c= Op::template op<C> ()) {
  typedef typename Series_variant (C) V;
  typedef implementation<series_recursive_abstractions,V> Ser;
  typedef typename Ser::template nullary_recursive_series_rep<Op,C,V> Nullary;
  series_rep<C>* rep= new Nullary (c);
  return recursive (series<C> (rep));
}

template<typename Op,typename C,typename V> inline Series
unary_recursive_series (const Series& f) {
  typedef implementation<series_recursive_abstractions,V> Ser;
  typedef typename Ser::template unary_recursive_series_rep<Op,C,V> Unary;
  Series_rep* rep= new Unary (f);
  return recursive (Series (rep));
}

template<typename Op,typename C,typename V> inline Series
unary_recursive_series (const Series& f, const C& c) {
  typedef implementation<series_recursive_abstractions,V> Ser;
  typedef typename Ser::template unary_recursive_series_rep<Op,C,V> Unary;
  Series_rep* rep= new Unary (f, c);
  return recursive (Series (rep));
}

template<typename Op,typename C,typename V> inline Series
binary_recursive_series (const Series& f, const Series& g) {
  typedef implementation<series_recursive_abstractions,V> Ser;
  typedef typename Ser::template binary_recursive_series_rep<Op,C,V> Binary;
  Series_rep* rep= new Binary (f, g);
  return recursive (Series (rep));
}

template<typename Op,typename C,typename V,typename X> inline Series
binary_scalar_recursive_series (const Series& f, const X& x) {
  typedef implementation<series_recursive_abstractions,V> Ser;
  typedef typename Ser::
    template binary_scalar_recursive_series_rep<Op,C,V,X> Binary;
  Series_rep* rep= new Binary (f, x);
  return recursive (Series (rep));
}

template<typename Op,typename C,typename V,typename X> inline Series
binary_scalar_recursive_series (const Series& f, const X& x, const C& init) {
  typedef implementation<series_recursive_abstractions,V> Ser;
  typedef typename Ser::
    template binary_scalar_recursive_series_rep<Op,C,V,X> Binary;
  Series_rep* rep= new Binary (f, x, init);
  return recursive (Series (rep));
}

template<typename Op,typename C,typename V,typename L> inline Series
unary_polynomial_recursive_series (const polynomial<L>& P) {
  typedef implementation<series_recursive_abstractions,V> Ser;
  typedef typename Ser::
    template unary_polynomial_recursive_series_rep<Op,C,V,L> Unary;
  Series_rep* rep= new Unary (P);
  return recursive (Series (rep));
}

template<typename Op,typename C,typename V,typename L> inline Series
unary_polynomial_recursive_series (const polynomial<L>& P, const C& init) {
  typedef implementation<series_recursive_abstractions,V> Ser;
  typedef typename Ser::
    template unary_polynomial_recursive_series_rep<Op,C,V,L> Unary;
  Series_rep* rep= new Unary (P, init);
  return recursive (Series (rep));
}

/******************************************************************************
* Function abstractions
******************************************************************************/

template<typename Op,typename C,typename V> inline Series
unary_series (const Series& f) {
  typedef implementation<series_abstractions,V> Ser;
  typedef typename Ser::template unary_series_rep<Op,C,V> Unary;
  return (Series_rep*) new Unary (f);
}

template<typename Op,typename C,typename V> inline Series
binary_series (const Series& f, const Series& g) {
  typedef implementation<series_abstractions,V> Ser;
  typedef typename Ser::template binary_series_rep<Op,C,V> Binary;
  return (Series_rep*) new Binary (f, g);
}

/******************************************************************************
* Coefficient-wise operations
******************************************************************************/

TMPL inline Series
operator + (const Series& f, const Series& g) {
  if (is_exact_zero (f)) return g;
  if (is_exact_zero (g)) return f;
  return binary_series<add_op> (f, g);
}

TMPL inline Series
operator + (const C& c, const Series& f) {
  if (is_exact_zero (c)) return f;
  if (is_exact_zero (f)) return Series (c);
  return binary_series<add_op> (Series (c), f);
}

TMPL inline Series
operator + (const Series& f, const C& c) {
  if (is_exact_zero (c)) return f;
  if (is_exact_zero (f)) return Series (c);
  return binary_series<add_op> (f, Series (c));
}

TMPL inline Series
operator - (const Series& f) {
  if (is_exact_zero (f)) return f;
  return unary_series<neg_op> (f);
}

TMPL inline Series
operator - (const Series& f, const Series& g) {
  if (is_exact_zero (f)) return -g;
  if (is_exact_zero (g)) return f;
  return binary_series<sub_op> (f, g);
}

TMPL inline Series
operator - (const C& c, const Series& f) {
  if (is_exact_zero (c)) return -f;
  if (is_exact_zero (f)) return Series (c);
  return binary_series<sub_op> (Series (c), f);
}

TMPL inline Series
operator - (const Series& f, const C& c) {
  if (is_exact_zero (f)) return Series (-c);
  if (is_exact_zero (c)) return f;
  return binary_series<sub_op> (f, Series (c));
}

TMPL inline Series
operator * (const Series& f, const C& c) {
  if (is_exact_zero (f) || is_exact_zero (c))
    return Series (get_format (c));
  return binary_scalar_series<rmul_op> (f, c);
}

TMPL inline Series
operator * (const C& c, const Series& f) {
  if (is_exact_zero (f) || is_exact_zero (c))
    return Series (get_format (c));
  return binary_scalar_series<lmul_op> (f, c);
}

/******************************************************************************
* Gcd -- over fields only
******************************************************************************/

TMPL
class gcd_series_rep: public Series_rep {
protected:
  Series f, g;
  int v; // gcd is z^v
public:
  inline gcd_series_rep (const Series& f2, const Series& g2):
    Series_rep (CF(f2)), f (f2), g (g2), v (-1) {}
  syntactic expression (const syntactic& z) const {
    return gcd (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 C next () {
    if (v >= 0) return this->zero ();
    if (f[this->n] != 0 || g[this->n] != 0) {
      v= this->n;
      return this->one ();
    }
    return this->zero (); }
};

TMPL inline Series
gcd (const Series& f, const Series& g) {
  return (Series_rep*) new gcd_series_rep<C,V> (f, g);
}

TMPL
class lcm_series_rep: public Series_rep {
protected:
  Series f, g;
  int v_f, v_g; // respective valuations
public:
  inline lcm_series_rep (const Series& f2, const Series& g2):
    Series_rep (CF(f2)), f(f2), g (g2), v_f (-1), v_g (-1) {}
  syntactic expression (const syntactic& z) const {
    return lcm (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 C next () {
    if (v_f >= 0 && v_g >= 0) return this->zero ();
    if (v_f < 0 && f[this->n] != 0) v_f= this->n;
    if (v_g < 0 && g[this->n] != 0) v_g= this->n;
    if (v_f >= 0 && v_g >= 0) return this->one ();
    return this->zero ();
  }
};

TMPL inline Series
lcm (const Series& f, const Series& g) {
  return (Series_rep*) new lcm_series_rep<C,V> (f, g);
}

/******************************************************************************
* Other linear time operations
******************************************************************************/

TMPL
class lshiftz_series_rep: public Series_rep {
  Series f;
  int shift;
public:
  lshiftz_series_rep (const Series& f2, int shift2):
    Series_rep (CF(f2)), f (f2), shift (shift2) {}
  syntactic expression (const syntactic& z) const {
    return lshiftz (flatten (f, z), flatten (shift)); }
  void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (f, max (0, ((int) l) - shift)); }
  C next () {
    return shift > ((int) this->n)? this->zero (): f[this->n - shift]; }
};

TMPL inline Series
lshiftz (const Series& f, const int& shift) {
  if (is_exact_zero (f)) return Series (CF(f));
  return (Series_rep*) new lshiftz_series_rep<C,V> (f, shift);
}

TMPL inline Series
rshiftz (const Series& f, const int& shift) {
  if (is_exact_zero (f)) return Series (CF(f));
  return (Series_rep*) new lshiftz_series_rep<C,V> (f, -shift);
}

TMPL
class restrict_series_rep: public Series_rep {
  Series f;
  nat start, end;
public:
  restrict_series_rep (const Series& f2, nat start2, nat end2):
    Series_rep (CF(f2)), f (f2), start (start2), end (end2) {}
  syntactic expression (const syntactic& z) const {
    return apply ("restrict", flatten (f, z),
                  flatten (start), flatten (end)); }
  void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (f, min (l, end)); }
  C next () {
    return ((this->n>=start && this->n<end)? f[this->n]: this->zero ()); }
};

TMPL inline Series
restrict (const Series& f, nat start, nat end) {
  if (is_exact_zero (f)) return Series (CF(f));
  return (Series_rep*) new restrict_series_rep<C,V> (f, start, end);
}

TMPL
class piecewise_series_rep: public Series_rep {
  Series f, g;
  nat pos;
public:
  piecewise_series_rep (const Series& f2, const Series &g2, nat pos2):
    Series_rep (CF(f2)), f (f2), g (g2), pos (pos2) {}
  syntactic expression (const syntactic& z) const {
    return apply ("piecewise", flatten (f, z), flatten (g, z),
                  flatten (pos)); }
  void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (f, min (l, pos));
    increase_order (g, l <= pos ? 0 : l); }
  C next () { return this->n < pos ? f[this->n] : g[this->n]; }
};

TMPL inline Series
piecewise (const Series& f, const Series& g, nat pos) {
  return (Series_rep*) new piecewise_series_rep<C,V> (f, g, pos);
}

TMPL
class derive_series_rep: public Series_rep {
  Series f;
public:
  derive_series_rep (const Series& f2):
    Series_rep (CF(f2)), f (f2) {}
  syntactic expression (const syntactic& z) const {
    return derive (flatten (f, z)); }
  void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (f, l+1); }
  C next () {
    return promote ((int) (this->n + 1), CF(f)) * f[this->n + 1]; }
};

TMPL inline Series
derive (const Series& f) {
  if (is_exact_zero (f)) return Series (CF(f));
  return (Series_rep*) new derive_series_rep<C,V> (f);
}

template<typename C, typename V, typename T> inline Series
derive_coeffs (const Series& f, const T& v) {
  return binary_scalar_series<derive_op> (f, v);
}

TMPL
class xderive_series_rep: public Series_rep {
  Series f;
public:
  xderive_series_rep (const Series& f2):
    Series_rep (CF(f2)), f (f2) {}
  syntactic expression (const syntactic& z) const {
    return xderive (flatten (f, z)); }
  void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (f, l); }
  C next () { return promote ((int) this->n, CF(f)) * f[this->n]; }
};

TMPL inline Series
xderive (const Series& f) {
  if (is_exact_zero (f)) return Series (CF(f));
  return (Series_rep*) new xderive_series_rep<C,V> (f);
}

TMPL
class integrate_series_rep: public Series_rep {
  Series f;
  C c;
public:
  integrate_series_rep (const Series& f2):
    Series_rep (CF(f2)), f (f2), c (promote (0, CF(f2))) {}
  integrate_series_rep (const Series& f2, const C& c2):
    Series_rep (CF(f2)), f (f2), c (c2) {}
  syntactic expression (const syntactic& z) const {
    if (c == promote (0, c)) return integrate (flatten (f, z));
    else return integrate_init (flatten (f, z), flatten (c)); }
  void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (f, max (0, ((int) l) - 1)); }
  C next () {
    if (this->n == 0) return c;
    return invert (promote ((int) this->n, c)) * f[this->n - 1]; }
};

TMPL inline Series
integrate (const Series& f) {
  if (is_exact_zero (f)) return Series (CF(f));
  return (Series_rep*) new integrate_series_rep<C,V> (f);
}

TMPL inline Series
integrate_init (const Series& f, const C& c) {
  return (Series_rep*) new integrate_series_rep<C,V> (f, c);
}

template<typename C, typename V, typename T> inline Series
integrate_coeffs (const Series& f, const T& v) {
  return binary_scalar_series<integrate_op> (f, v);
}

TMPL
class shift_series_rep: public Series_rep {
  // heuristic shifting
  Series f;
  C c;
  nat order;
public:
  shift_series_rep (const Series& f2, const C& c2, nat order2):
    Series_rep (CF(f2)), f (f2), c (c2), order (order2) {}
  syntactic expression (const syntactic& z) const {
    return apply ("shift", flatten (f, z), flatten (c)); }
  void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (f, l + order); }
  C next () {
    C s= this->zero (), p= this->one ();
    for (nat i=0; i<order; i++) {
      s += p * f[this->n + i];
      p *= c * (promote ((int) (this->n + i + 1), c) /
                promote ((int) (i + 1), c));
    }
    return s; }
};

TMPL inline Series
shift (const Series& f, const C& q, nat order) {
  if (is_exact_zero (f)) return Series (CF(f));
  return (Series_rep*) new shift_series_rep<C,V> (f, q, order);
}

TMPL inline Series
series_shift_default (const Series& s, const C& sh) {
  return shift (s, sh, mmx_bit_precision);
}

TMPL
class q_difference_series_rep: public Series_rep {
  Series f;
  C q;
  C p;
public:
  q_difference_series_rep (const Series& f2, const C& q2):
    Series_rep (CF(f2)), f (f2), q (q2), p (promote (1, q2)) {}
  syntactic expression (const syntactic& z) const {
    return apply ("q_difference", flatten (f, z), flatten (q)); }
  void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (f, l); }
  C next () {
    C c= p * f[this->n];
    p *= q;
    return c; }
};

TMPL inline Series
q_difference (const Series& f, const C& q) {
  if (is_exact_zero (f)) return Series (CF(f));
  return (Series_rep*) new q_difference_series_rep<C,V> (f, q);
}

TMPL
class dilate_series_rep: public Series_rep {
  Series f;
  nat p;
public:
  dilate_series_rep (const Series& f2, int p2):
    Series_rep (CF(f)), f (f2), p (p2) {}
  syntactic expression (const syntactic& z) const {
    return apply ("dilate", flatten (f, z), flatten (p)); }
  void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (f, l/p); }
  C next () {
    return this->n % p == 0? f[this->n / p]: this->zero (); }
};

TMPL inline Series
dilate (const Series& f, nat p) {
  if (is_exact_zero (f)) return Series (CF(f));
  return (Series_rep*) new dilate_series_rep<C,V> (f, p);
}

TMPL
class deflate_series_rep: public Series_rep {
  Series f;
  nat p;
public:
  deflate_series_rep (const Series& f2, int p2):
    Series_rep (CF(f)), f (f2), p (p2) {}
  syntactic expression (const syntactic& z) const {
    return apply ("deflate", flatten (f, z), flatten (p)); }
  void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (f, l*p); }
  C next () { return f[(this->n)*p]; }
};

TMPL inline Series
deflate (const Series& f, nat p) {
  if (is_exact_zero (f)) return Series (CF(f));
  return (Series_rep*) new deflate_series_rep<C,V> (f, p);
}

/******************************************************************************
* Multiplication and division
******************************************************************************/

TMPL inline Series
operator * (const Series& f, const Series& g) {
  typedef implementation<series_multiply,V> Ser;
  return Ser::ser_mul (f, g);
}

TMPL inline Series
truncate_mul (const Series& f, const Series& g, nat nf, nat ng) {
  // Product of the only nf first terms of f by the ng first ones of g
  typedef implementation<series_multiply,V> Ser;
  return Ser::ser_truncate_mul (f, g, nf, ng);
}

TMPL inline Series
operator / (const Series& f, const C& c) {
  typedef implementation<series_divide,V> Ser;
  return Ser::ser_rdiv_sc (f, c);
}

TMPL inline Series
operator / (const Series& f, const Series& g) {
  typedef implementation<series_divide,V> Ser;
  return Ser::ser_div (f, g);
}

TMPL inline Series
operator / (const C& c, const Series& f) {
  typedef implementation<series_divide,V> Ser;
  if (is_exact_zero (c)) return Series (get_format (c));
  return Ser::ser_div (Series (c), f);
}

TMPL inline Series
invert (const Series& f) {
  return promote (1, CF(f)) / f;
}

TMPL inline Series
quo (const Series& f, const C& c) {
  typedef implementation<series_divide,V> Ser;
  return Ser::ser_rquo_sc (f, c);
}

TMPL inline Series
quo (const Series& f, const Series& g) {
  typedef implementation<series_divide,V> Ser;
  return Ser::ser_quo (f, g);
}

TMPL inline Series
rem (const Series& f, const C& c) {
  typedef implementation<series_divide,V> Ser;
  return Ser::ser_rrem_sc (f, c);
}

TMPL inline Series
rem (const Series& f, const Series& g) {
  if (is_exact_zero (f)) return Series (CF(f));
  return f - g * quo (f, g);
}

TMPL inline bool
divides (const Series& f, const Series& g) {
  return val (f) <= val (g);
}

ARITH_SCALAR_INT_SUGAR(TMPL,Series);

/******************************************************************************
* Separable nth root
******************************************************************************/

TMPL inline Series
separable_root (const Series& f, nat n) {
  typedef implementation<series_separable_root,V> Ser;
  return Ser::sep_root (f, n);
}

TMPL inline Series
separable_root_init (const Series& f, nat n, const C& init) {
  typedef implementation<series_separable_root,V> Ser;
  return Ser::sep_root (f, n, init);
}

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

template<typename C, typename V, typename L> inline Series
polynomial_regular_root (const polynomial<L>& P) {
  typedef implementation<series_polynomial_regular_root,V> Ser;
  if (P[0] == L(0))  return Series ();
  return Ser::template pol_root<C,V,L> (P);
}

template<typename C, typename V, typename L> inline Series
polynomial_regular_root_init (const polynomial<L>& P, const C& init) {
  typedef implementation<series_polynomial_regular_root,V> Ser;
  if (P[0] == L(0))  return Series ();
  return Ser::template pol_root<C,V,L> (P, init);
}
  
/******************************************************************************
* pth root where p is the characteristic of the residual ring
******************************************************************************/

TMPL inline Series
pth_root (const Series& f) {
  typedef implementation<series_pth_root,V> Ser;
  return Ser::unsep_root (f);
}

/******************************************************************************
* Composition and reversion
******************************************************************************/

TMPL inline Series
compose (const Series& f, const Series& g) {
  typedef implementation<series_compose,V> Ser;
  return Ser::ser_compose (f, g);
}

TMPL inline Series
reverse (const Series& f) {
  typedef implementation<series_compose,V> Ser;
  return Ser::ser_reverse (f);
}

/******************************************************************************
* Reliable and precision related routines
******************************************************************************/

TMPL
class change_precision_series_rep: public Series_rep {
  Series f;
  nat p;
public:
  inline change_precision_series_rep (const Series& f2, nat p2):
    Series_rep (CF(f2)), f (f2), p (p2) {}
  syntactic expression (const syntactic& z) const {
    return syn ("change_precision", flatten (f, z), flatten (p)); }
  void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (f, l); }
  C next () {
    return change_precision (f[this->n], p); }
};

TMPL inline Series
change_precision (const Series& f, nat p) {
  return (Series_rep*) new change_precision_series_rep<C,V> (f, p);
}

TMPL inline nat
precision (const Series& f) {
  return precision (f[0]);
}

TMPL inline Series sharpen (const Series& f) {
  return unary_map<sharpen_op> (f); }
TMPLK inline Series blur (const Series& f, const K& x) {
  return binary_map_scalar<blur_op> (f, x); }

/******************************************************************************
* Operations on series like types
******************************************************************************/

BINARY_RETURN_TYPE(TMPL,truncate_op,Series,nat,Polynomial);
UNARY_RETURN_TYPE(TMPL,complete_op,Series,Series);

TMPL inline Series
complete (const Series& f) { return f; }

BINARY_RETURN_TYPE(TMPL,truncate_op,PolynomialV,nat,PolynomialV);
UNARY_RETURN_TYPE(TMPL,complete_op,PolynomialV,series<C>);

TMPL inline PolynomialV
truncate (const PolynomialV& p, nat sigma) { return range (p, 0, sigma); }

TMPL inline series<C>
complete (const PolynomialV& p) {
  return (series_rep<C>*)
    new polynomial_series_rep<C,typename Series_variant (C)>
      (as<Polynomial> (p)); }

template<typename N,typename D> struct quotient;
#define Quotient quotient<PolynomialV,PolynomialV>
UNARY_RETURN_TYPE(TMPL,complete_op,Quotient,series<C>);
TMPL inline series<C>
complete (const Quotient& f) {
  return complete (numerator (f))
       / complete (denominator (f)); }
#undef Quotient

#undef TMPL
#undef TMPLK
#undef Format
#undef Series
#undef Series_rep
#undef Recursive_series_rep
#undef Polynomial
#undef PolynomialV
#undef Truncated_series
#undef Completed_series
#undef Truncated_polynomial
#undef Completed_polynomial
} // namespace mmx
#endif // __MMX_SERIES_HPP

---