/Users/mourrain/Devel/mmx/algebramix/include/algebramix/series_naive.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_NAIVE__HPP
#define __MMX_SERIES_NAIVE__HPP
#include <algebramix/polynomial.hpp>

namespace mmx {
#define TMPL template<typename C, typename V>
#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>

template<typename C>
struct series_variant_helper;

template <typename C, typename V=typename series_variant_helper<C>::SV>
class series_rep;

template <typename C, typename V=typename series_variant_helper<C>::SV>
class series;

template <typename C, typename V=typename series_variant_helper<C>::SV>
class recursive_series_rep;

/******************************************************************************
* Naive variant for series
******************************************************************************/

struct series_naive {};

template<typename C>
struct series_variant_helper {
  typedef series_naive SV;
};

template<typename C, typename V>
struct series_polynomial_helper {
  typedef typename Polynomial_variant(C) PV;
};

/******************************************************************************
* Global variables, default is R[z]
******************************************************************************/

struct series_defaults {};

template<typename V>
struct implementation<series_defaults,V,series_naive> {

  template<typename S>
  class global_variables {

    static inline generic& dyn_name () {
      static generic name = "z";
      return name; }

    static inline nat& dyn_output_order () {
      static nat n = 10;
      return n; }

    static inline nat& dyn_cancel_order () {
      static nat n = 5;
      return n; }

    static inline bool& dyn_formula_output () {
      static bool b = false;
      return b; }

  public:
    static inline void set_variable_name (const generic& x) {
      dyn_name () = x; }
    static inline generic get_variable_name () {
      return dyn_name (); }
    static inline void set_output_order (const nat& x) {
      dyn_output_order () = x; }
    static inline nat get_output_order () {
      return dyn_output_order (); }
    static inline void set_cancel_order (const nat& x) {
      dyn_cancel_order () = x; }
    static inline nat get_cancel_order () {
      return dyn_cancel_order (); }
    static inline void set_formula_output (const bool& x) {
      dyn_formula_output () = x; }
    static inline bool get_formula_output () {
      return dyn_formula_output (); }
  };
}; // implementation<series_defaults,V,series_naive>

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

struct series_scalar_abstractions {};

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

template<typename Op,typename C,typename V,typename X>
class binary_scalar_series_rep: public Series_rep {
protected:
  Series f;
  X x;
public:
  inline binary_scalar_series_rep (const Series& f2, const X& x2):
    Series_rep (CF(f2)), f(f2), x (x2) {}
  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 C next () { return Op::op (f[this->n], x); }
};

template<typename Op,typename C,typename V,typename X,typename Y>
class ternary_scalar_series_rep: public Series_rep {
protected:
  Series f;
  X x;
  Y y;
public:
  inline ternary_scalar_series_rep (const Series& f2, const X& a, const Y& b):
    Series_rep (CF(f2)), f (f2), x (a), y (b) {}
  syntactic expression (const syntactic& z) const {
    return Op::op (flatten (f, z), flatten (x), flatten (y)); }
  virtual void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (f, l); }
  virtual C next () { return Op::op (f[this->n], x, y); }
};
}; // implementation<series_abstractions,V,series_naive>

/******************************************************************************
* Abstract maps on series
******************************************************************************/

struct series_map_as_abstractions {};

template<typename U>
struct implementation<series_map_as_abstractions,U,series_naive>:
  public implementation<series_scalar_abstractions,U>
{

template<typename Op, typename C, typename V, typename S, typename SV>
class unary_map_as_series_rep: public Series_rep {
protected:
  series<S,SV> f;
public:
  inline unary_map_as_series_rep (const series<S,SV>& f2, const Format& fm):
    Series_rep (fm), f (f2) {}
  inline unary_map_as_series_rep (const series<S,SV>& f2):
    Series_rep (unary_map<Op> (CF(f2))), f (f2) {}
  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 C next () {
    C r= get_sample (this->tfm ());
    Op::set_op (r, f[this->n]);
    return r; }
};

template<typename Op, typename C, typename V, typename S, typename SV> 
static inline Series
unary_map_as (const series<S,SV>& f, const Format& fm) {
  typedef unary_map_as_series_rep<Op,C,V,S,SV> Unary;
  return (Series_rep*) new Unary (f, fm);
}

template<typename Op, typename C, typename V, typename S, typename SV> 
static inline Series
unary_map_as (const series<S,SV>& f) {
  return unary_map_as<Op,C,V,S,SV> (f, unary_map<Op> (CF(f)));
}

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

/******************************************************************************
* Abstract recursive operations on series
******************************************************************************/

struct series_recursive_abstractions {};

template<typename U>
struct implementation<series_recursive_abstractions,U,series_naive>:
  public implementation<series_map_as_abstractions,U>
{

template<typename Op,typename C,
         typename V=typename series_variant_helper<C>::SV>
class nullary_recursive_series_rep: public Recursive_series_rep {
protected:
  C c;
public:
  nullary_recursive_series_rep (const C& c2):
    Recursive_series_rep (get_format (c)), c (c2) {}
  syntactic expression (const syntactic& z) const {
    (void) z;
    return Op::op_init (flatten (c)); }
  Series initialize () {
    this->initial (0)= c;
    return Op::def (Series (this->me ())); }
};

template<typename Op,typename C,typename V>
class unary_recursive_series_rep: public Recursive_series_rep {
protected:
  Series f;
  bool with_init;
  C c;
public:
  unary_recursive_series_rep (const Series& f2):
    Recursive_series_rep (CF(f2)), f (f2), with_init (false) {}
  unary_recursive_series_rep (const Series& f2, const C& c2):
    Recursive_series_rep (CF(f2)), f (f2), with_init (true), c (c2) {}
  syntactic expression (const syntactic& z) const {
    return Op::op (flatten (f, z)); }
  virtual void Increase_order (nat l) {
    Recursive_series_rep::Increase_order (l);
    increase_order (f, l); }
  Series initialize () {
    if (with_init) this->initial (0)= c;
    else {
      ASSERT (Op::nr_init () <= 1, "wrong number of initial conditions");
      if (Op::nr_init () == 1)
        this->initial (0)= Op::op (f[0]);
    }
    return Op::def (Series (this->me ()), f); }
};

template<typename Op,typename C,typename V>
class binary_recursive_series_rep: public Recursive_series_rep {
protected:
  Series f, g;
public:
  binary_recursive_series_rep (const Series& f2, const Series& g2):
    Recursive_series_rep (CF(f2)), f(f2), g(g2) {}
  syntactic expression (const syntactic& z) const {
    return Op::op (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 () {
    ASSERT (Op::nr_init () <= 1, "wrong number of initial conditions");
    if (Op::nr_init () == 1) this->initial (0)= Op::op (f[0], g[0]);
    return Op::def (Series (this->me ()), f, g);
  }
};

template<typename Op,typename C,typename V,typename X>
class binary_scalar_recursive_series_rep: public Recursive_series_rep {
protected:
  Series f;
  bool with_init;
  X x;
  C c;

public:
  binary_scalar_recursive_series_rep (const Series& f2, const X& x2):
    Recursive_series_rep (CF(f2)), f (f2), with_init (false), x (x2) {}
  binary_scalar_recursive_series_rep (const Series& f2, const X& x2,
                                      const C& c2):
    Recursive_series_rep (CF(f2)), f (f2), with_init (true), x (x2), c (c2) {}
  syntactic expression (const syntactic& z) const {
    return Op::op (flatten (f, z), flatten (x)); }
  virtual void Increase_order (nat l) {
    Recursive_series_rep::Increase_order (l);
    increase_order (f, l); }
  Series initialize () {
    if (with_init) this->initial (0)= c;
    else {
      ASSERT (Op::nr_init () <= 1, "wrong number of initial conditions");
      if (Op::nr_init () == 1)
        this->initial (0)= Op::op (f[0], x);
    }
    return Op::def (Series (this->me ()), f, x); }
};

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

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

struct series_abstractions {};

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

template<typename Op,typename C,typename V>
class unary_series_rep: public Series_rep {
protected:
  Series f;
public:
  inline unary_series_rep (const Series& f2):
    Series_rep (CF(f2)), f(f2) {}
  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 C next () { return Op::op (f[this->n]); }
};

template<typename Op,typename C,typename V>
class binary_series_rep: public Series_rep {
protected:
  Series f, g;
public:
  inline binary_series_rep (const Series& f2, const Series& g2):
    Series_rep (CF(f2)), f(f2), g (g2) {}
  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 C next () { return Op::op (f[this->n], g[this->n]); }
};

}; // implementation<series_abstractions,V,series_naive>
  
/******************************************************************************
* Multiplication
******************************************************************************/

struct series_multiply {};

template<typename U>
struct implementation<series_multiply,U,series_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,C,V> {
public:
  inline mul_series_rep (const Series& f, const Series& g):
    Ser::template binary_series_rep<mul_op,C,V> (f, g) {}
  C next () {
    nat i;
    C acc= this->f[0] * this->g[this->n];
    for (i=1; i<=this->n; i++)
      acc += this->f[i] * this->g[this->n-i];
    return acc; }
};

TMPL
class truncate_mul_series_rep:
    public Ser::template binary_series_rep<mul_op,C,V> {
  nat nf, ng;
public:
  inline truncate_mul_series_rep (const Series& f, const Series& g,
                                  nat nf2, nat ng2):
    Ser::template binary_series_rep<mul_op,C,V> (f, g), nf (nf2), ng (ng2) {}
  C next () {
    nat i, n= this->n;
    C acc= this->zero ();
    for (i= n >= ng ? n - ng + 1 : 0; i < min (1 + n, nf); i++)
      acc += this->f[i] * this->g[n-i];
    return acc; }
};

TMPL static inline Series
ser_mul (const Series& f, const Series& g) {
  typedef mul_series_rep<C,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<C,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_naive>

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

struct series_divide {};

template<typename U>
struct implementation<series_divide,U,series_naive>:
  public implementation<series_multiply,U>
{

TMPL static inline Series
ser_rdiv_sc (const Series& f, const C& c) {
  typedef implementation<series_scalar_abstractions,U> Ser;
  typedef typename Ser::
    template binary_scalar_series_rep<rdiv_op,C,V,C> Div_sc_rep;
  if (is_exact_zero (f)) return Series (CF(f));
  return (Series_rep*) new Div_sc_rep (f, c); }

TMPL static inline Series
ser_div (const Series& f, const Series& g) {
  typedef implementation<series_recursive_abstractions,U> Ser;
  typedef typename Ser::
    template binary_recursive_series_rep<div_op,C,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 C& c) {
  typedef implementation<series_scalar_abstractions,U> Ser;
  typedef typename Ser
    ::template binary_scalar_series_rep<rquo_op,C,V,C> Rquo_rep;
  if (is_exact_zero (f)) return Series (CF(f));
  return (Series_rep*) new Rquo_rep (f, c); }

TMPL static inline Series
ser_quo (const Series& f, const Series& g) {
  typedef implementation<series_recursive_abstractions,U> Ser;
  typedef typename Ser
    ::template binary_recursive_series_rep<quo_op,C,V> Quo_rep;
  if (is_exact_zero (f)) return Series (CF(f));
  return recursive (Series ((Series_rep*) new Quo_rep (f, g))); }

TMPL static inline Series
ser_rrem_sc (const Series& f, const C& c) {
  typedef implementation<series_scalar_abstractions,U> Ser;
  typedef typename Ser
     ::template binary_scalar_series_rep<rrem_op,C,V> Rrem_rep;
  if (is_exact_zero (f)) return Series (CF(f));
  return (Series_rep*) new Rrem_rep (f, c); }

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

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

struct series_compose {};

template<typename U>
struct implementation<series_compose,U,series_naive>:
    public implementation<series_multiply,U>
{

TMPL
class compose_series_rep: public Series_rep {
  Series f;
  Series g;
  Series cumul;
  Series power;
public:
  inline compose_series_rep (const Series& f2, const Series& g2):
    Series_rep (CF(f2)), f (f2), g (g2),
    cumul (CF(f2)), power (promote (1, CF(f2))) {}
  syntactic expression (const syntactic& z) const {
    return apply (GEN_COMPOSE, flatten (f, z), flatten (g, z)); }
  void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (f, l);
    increase_order (g, l);
    increase_order (cumul, l);
    increase_order (power, l); }
  C next () {
    cumul += f[this->n] * power;
    power *= g;
    return cumul[this->n];
  }
};

TMPL static inline Series
ser_compose (const Series& f, const Series& g) {
  typedef compose_series_rep<C,V> Compose_rep;
  if (is_exact_zero (f)) return Series (CF(f));
  ASSERT (g[0] == 0, "constant coefficient must be zero");  
  return (Series_rep*) new Compose_rep (f, g); }

TMPL
class reverse_series_rep: public Recursive_series_rep {
protected:
  Series f;
public:
  reverse_series_rep (const Series& f2):
    Recursive_series_rep (CF(f2)), f (f2) {}
  syntactic expression (const syntactic& z) const {
    return apply (GEN_REVERSE, flatten (f, z)); }
  virtual void Increase_order (nat l) {
    Recursive_series_rep::Increase_order (l);
    increase_order (f, l); }
  Series initialize () {
    ASSERT (f[0] == 0 && f[1] != 0, "invalid reversion");
    C      a= f[1];
    Series z (this->one (), 1);
    Series s= square (rshiftz (this->me (), 1));
    Series c= compose (rshiftz (this->f, 2), this->me ());
    return (z - lshiftz (c * s, 2)) / a; }
};

TMPL static inline Series
ser_reverse (const Series& f) {
  typedef reverse_series_rep<C,V> Reverse_rep;
  Series_rep* rep= new Reverse_rep (f);
  return recursive (Series (rep)); }

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

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

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

private:
  template<typename S> static S
  binpow_no_tangent (const S& me, nat r) {
    typedef Scalar_type(S) C;
    // me^r - r * (me - me[0]) * me[0]^(r-1) - me[0]^r
    if (r <= 1) return S (0);
    if (r == 2) return lshiftz (square (rshiftz (me, 1)), 2);

    nat h= r >> 1;
    C x= me[0];
    C x_h_1= binpow (x, h-1);
    S y= me - me[0];
    S z= rshiftz (binpow_no_tangent (me, h), 1);
    S t0= as<C> (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) {
      nat h2= h << 1;
      w= rshiftz (w, 1);
      w= lshiftz ((y + x) * w, 1)
        + lshiftz (as<C> (h2) * square (rshiftz (y * x_h_1, 1)) * x, 2);
    }
    return w; }

public:
  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) C;
    C me0_r_1= binpow (me[0], r - 1);
    C me0_r= me[0] * me0_r_1;
    return (a - me0_r - binpow_no_tangent (me, r)) / (as<C> (r) * me0_r_1); }

}; // struct ser_separable_root_op

struct series_separable_root {};

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

TMPL static inline Series
sep_root (const Series& a, nat r) {
  typedef implementation<series_recursive_abstractions,U> Ser;
  typedef typename Ser::template binary_scalar_recursive_series_rep
    <ser_separable_root_op,C,V,nat> Root;
  ASSERT (C(r) != 0, "wrong argument");
  if (is_exact_zero (a)) return Series (CF(a));
  Series_rep* rep= new Root (a, r);
  return recursive (Series (rep)); }

TMPL static inline Series
sep_root (const Series& a, nat r, const C& init) {
  typedef implementation<series_recursive_abstractions,U> Ser;
  typedef typename Ser::template binary_scalar_recursive_series_rep
    <ser_separable_root_op,C,V,nat> Root;
  ASSERT (C(r) != 0, "wrong argument");
  if (is_exact_zero (a)) return Series (CF(a));
  Series_rep* rep= new Root (a, r);
  return recursive (Series (rep)); }

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

/******************************************************************************
* pth root where p is the characteristic
******************************************************************************/

struct series_pth_root_reg {};
struct series_pth_root {};

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

TMPL static inline Series
unsep_root (const Series& f) { ERROR ("not implemented"); }

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

#undef TMPL
#undef Format
#undef Series
#undef Series_rep
#undef Recursive_series_rep
#undef Polynomial
} // namespace mmx
#endif // __MMX_SERIES_NAIVE__HPP