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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : series_vector.hpp
* DESCRIPTION: Vectors and matrices of univariate power series
* COPYRIGHT  : (C) 2011  J. Berthomieu, J. van der Hoeven and R. 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_VECTOR_HPP
#define __MMX_SERIES_VECTOR_HPP
#include <basix/vector.hpp>
#include <algebramix/series.hpp>
#include <algebramix/matrix.hpp>
#include <algebramix/series_matrix.hpp>

namespace mmx {
#define TMPL template<typename C, typename V, typename W>
#define TMPLU template<typename C, typename V, typename U>
#define TMPLWU template<typename C, typename V, typename W, typename U>
#define Series series<C,V>
#define Vector vector<C,W>
#define Matrix matrix<C,U>
#define Series_rep series_rep<C,V>
#define Series_vector series<Vector,V>
#define Vector_series vector<Series,W>
#define Series_matrix series<Matrix,V>
#define Matrix_series matrix<Series,U>

/******************************************************************************
* Variants
******************************************************************************/

template<typename C, typename V>
struct series_variant_helper<vector<C,V> > {
  typedef typename series_variant_helper<C>::SV SV;
};

/******************************************************************************
* Vectors of series
******************************************************************************/

STYPE_TO_TYPE(TMPL,as_vector_type,Series_vector,Vector_series);

TMPL
class vector_access_series_rep: public Series_rep {
  const Series_vector f;
  const nat i;
public:
  vector_access_series_rep (const Series_vector& f2, nat i2):
    Series_rep (get_format1 (CF(f2))), f (f2), i (i2) {}
  syntactic expression (const syntactic& z) const {
    return access (flatten (f, z), flatten (i)); }
  void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (f, l); }
  C next () {
    if (is_a_scalar (f[this->n])) return C();
    return i < N(f[this->n]) ? f[this->n][i] : C(); }
};

TMPL Series
access (const Series_vector& f, nat i) {
  return (Series_rep*) new vector_access_series_rep<C,V,W> (f, i);
}

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

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

TMPL syntactic
flatten (const Vector_series& v, const syntactic& z) {
  syntactic g;
  nat i, n= N(v);
  vector<syntactic,W> r (g, n);
  for (i=0; i<n; i++)
    r[i]= flatten (v[i], z);
  return flatten (r);
}

TMPL format<Vector >
get_vector_format (const Vector_series& v) {
  format<Series > fm1= CF(v);
  format<C> fm2= get_format1 (fm1);
  C zero= promote (0, fm2);
  Vector r (zero, N(v));
  return get_format (r);
}

TMPL
class vector_series_rep: public series_rep<Vector,V> {
  const Vector_series v;
public:
  vector_series_rep (const Vector_series& v2):
    series_rep<Vector,V> (get_vector_format (v2)), v (v2) {}
  syntactic expression (const syntactic& z) const {
    return flatten (v, z); }
  void Increase_order (nat l) {
    series_rep<Vector,V>::Increase_order (l);
    for (nat i=0; i<N(v); i++)
      increase_order (v[i], l); }
  Vector next () {
    nat i, n= N(v);
    Vector coeff (promote (0, get_format1 (CF(v))), n);
    for (i=0; i<n; i++)
      coeff[i]= v[i][this->n];
    return coeff; }
};

TMPL Series_vector
as_series (const Vector_series& v) {
  return (series_rep<Vector,V>*) new vector_series_rep<C,V,W> (v);
}

TMPL Series_vector
from_vector (const Vector_series& v) {
  return (series_rep<Vector,V>*) new vector_series_rep<C,V,W> (v);
}

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

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

/******************************************************************************
* Linear differential equations
******************************************************************************/

TMPL Series_vector
solve_vector_lde_init (const Series_vector& f, const Vector& c) {
  return unary_recursive_series<solve_vector_lde_op> (f, c);
}

TMPL Vector_series
solve_lde_init (const Vector_series& f, const Vector& c) {
  return as_vector (solve_vector_lde_init (as_series (f), c));
}

TMPL Series // FIXME: obsolete
solve_lde (const Vector_series& f, const Vector& c) {
  return read (solve_lde_init (f, c), 0);
}

/******************************************************************************
* Linear equations resolution
******************************************************************************/
#define D       typename C::C

template<typename C, typename V, typename X>
class carry_mul_quorem_series_rep: public series_rep<vector<C>,V> {
protected:
  Series f;
  C x;
  D q;
public:
  inline carry_mul_quorem_series_rep (const X& x2, const Series& f2):
    series_rep<vector<C>,V> (get_format (vec<C> (get_sample (CF(f2))))),
    f(f2), x (as<C> (x2)), q (0) {}
  syntactic expression (const syntactic& z) const {
    return apply (syntactic ("carry_mul_quorem"), flatten (f, z),
                  flatten (x)); }
  virtual void Increase_order (nat l) {
    series_rep<vector<C>,V>::Increase_order (l);
    increase_order (f, l); }
  virtual vector<C> next () {
    C r= C (rem (x.rep * f[this->n].rep, * C::get_modulus (), q));
    return vec<C> (C (q), r); }
};

template<typename C, typename V> static inline Series
carry_mul_rem_series (const C& c, const Series& f, Series& q) {
  series<vector<C>,V> s= (series_rep<vector<C>,V> *)
    new carry_mul_quorem_series_rep<C,V,C> (c, f);
  vector<Series> v= as_vector (s);
  q= v[0];
  return v[1];
}

TMPLWU
static inline Vector_series
matrix_mul_nocarry (const Matrix& mat, const Vector_series& v) {
  typedef typename Series_variant(D) SV;
  return as<vector<Series, W> > (as<matrix<series<C,SV>, U> > (mat) * 
                                 as<vector<series<C,SV>, W> > (v));
}

template<typename C, typename V>
class carry_add_quorem_series_rep: public series_rep<vector<C>,V> {
protected:
  Series a, b;
  D r;
public:
  inline carry_add_quorem_series_rep (const Series& a2, const Series& b2):
    series_rep<vector<C>,V> (get_format (vec<C> (get_sample (CF(a2))))),
    a(a2), b (b2), r (0) {}
  syntactic expression (const syntactic& z) const {
    return apply (syntactic ("carry_add_quorem"), flatten (a, z),
                  flatten (b,z)); }
  virtual void Increase_order (nat l) {
    series_rep<vector<C>,V>::Increase_order (l);
    increase_order (a, l);
    increase_order (b, l); }
  virtual vector<C> next () {
    D r = a[this->n].rep + b[this->n].rep;
    if (r < * C::get_modulus ())
      return vec<C> (C(0), C(r));
    return vec<C> (C(1), C(r- (* C::get_modulus ())) ); }
};


template<typename C, typename V> static inline Series
carry_add_rem (const Series& a, const Series& b, Series& q) {
  series<vector<C>,V> s= (series_rep<vector<C>,V> *)
    new carry_add_quorem_series_rep<C,V> (a, b);
  vector<Series> v= as_vector (s);
  q= v[0];
  return v[1];
}

template<typename C, typename V>
class carry_special_add_series_rep: public Series_rep {
protected:
  Series a, b, c;
  D r;
public:
  inline carry_special_add_series_rep (const Series& a2, const Series& b2,
                                       const Series& c2):
    Series_rep (CF(a2)), a(a2), b (b2), c (c2), r (0) {}
  syntactic expression (const syntactic& z) const {
    return apply (apply (syntactic ("carry_special_add"), flatten (a, z),
                         flatten (b,z)), flatten (c,z)); }
  virtual void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (a, l);
    increase_order (b, l);
    increase_order (c, l); }
  virtual C next () {
    D s = a[this->n].rep + b[this->n].rep + c[this->n].rep + r;
    if (s < * C::get_modulus ()) {
      r = 0;
      return C (s);
    }
    r = 1;
    return C(s - (* C::get_modulus ())); }
};

template<typename C, typename V> static inline Series
carry_special_add (const Series& a, const Series& b, const Series& c) {
  return (Series_rep *) new carry_special_add_series_rep<C,V> (a, b, c);
}


TMPLWU
static Vector_series
matrix_mul_quo (const Matrix& mat, const Vector_series& v) {
  nat m= rows(mat), n= cols(mat);
  ASSERT (N(v) == n, "sizes don't match");
  if (is_a_scalar (mat)) return as<Series> (mat.scalar()) * v;
  ASSERT (is_non_scalar (v), "non-scalar vector expected");
  Series zero (get_format1 (CF(v)));
  Vector_series w (zero, m);
  Series mr, mq, c, ar, aq;
  for (nat k=0; k < m; k++) {
    mr = carry_mul_rem_series (mat(k,0), v[0], mq);
    c  = mq;
    ar = mr;
    for (nat l=1; l < n; l++) {
      if (mat (k,l) == C(0)) {
        mr = zero;
        mq = zero;
        aq = zero;
        continue;
      }
      mr = carry_mul_rem_series (mat(k,l), v[l], mq);
      ar = carry_add_rem (ar, mr, aq);
      c  = carry_special_add (c, mq, aq);  
    }
    w[k] = c;
  }
  return w;
}


TMPLWU
class ldiv_sc_mat_series_rep: public recursive_series_rep<Vector,V> {
protected:
  Vector_series f;
  Matrix c;
public:
  ldiv_sc_mat_series_rep (const Matrix& c2, const Vector_series& f2):
    recursive_series_rep<Vector,V> (get_vector_format (f2)),
    f(f2), c(c2) {}
  syntactic expression (const syntactic& z) const {
    return invert (flatten (c)) * flatten (f,z); }
  virtual void Increase_order (nat l) {
    recursive_series_rep<Vector,V>::Increase_order (l);
    for (nat i = 0; i < N(f); i++)
      increase_order (f[i],l);
  }
  Series_vector initialize () {
    Matrix u= invert (c);
    this->initial (0)= u * as_series (f)[0];
    Series_vector ls_me = lshiftz_series_vector (this -> me (), N(f));
    Vector_series tmp= f - matrix_mul_quo (c, as_vector (ls_me));
    return as_series (matrix_mul_nocarry (u, tmp));
  }
};


TMPLWU static inline Vector_series
ser_ldiv_sc_mat (const Matrix& c, const Vector_series& f) {
  typedef ldiv_sc_mat_series_rep<C,V,W,U> Div_sc_mat_rep;
  if (is_exact_zero (as_series (f)))
    return Vector_series (Series (CF(c)), cols(c));
  Series_vector s = (series_rep<Vector,V>*) new Div_sc_mat_rep (c,f);
  return as_vector (recursive (s)); 
}



TMPLWU
class ldiv_mat_series_rep: public recursive_series_rep<Vector,V> {
protected:
  Vector_series f;
  Matrix_series mat;
public:
  ldiv_mat_series_rep (const Matrix_series& mat2, const Vector_series& f2):
    recursive_series_rep<Vector,V> (get_vector_format (f2)),
    f(f2), mat(mat2) {}
  syntactic expression (const syntactic& z) const {
    return invert (flatten (mat)) * flatten (f,z); }
  virtual void Increase_order (nat l) {
    recursive_series_rep<Vector,V>::Increase_order (l);
    for (nat i = 0; i < N(f); i++) {
      increase_order (f[i],l);
      for (nat j = 0; j < cols(mat); j++)
        increase_order (mat(i,j),l);
    }
  }
  Series_vector initialize () {
    Series_matrix smat = as_series (mat);
    //TODO create a matrix variant helper for Matrix_series
    //matrix<C> smat0 = smat[0]; 
    this->initial (0)= invert (smat[0]) * as_series (f)[0];
    
    Vector_series v = f - 
      as_vector (lshiftz_series_vector (as_series (as_matrix (rshiftz (smat))*
                                                   as_vector (this-> me ())),
                                        N(f)));
    //Conversions series <-> vector, matrix dues à
    //-produit matrice vecteur prend en entrée matrix, vector
    //-rshiftz sur une matrice prend en entrée series
    //-appel à ser_ldiv_sc_smat prend en entrée matrix, vector
    return as_series (ser_ldiv_sc_mat (smat[0], v));
    /* Pour essayer de simplifier l'appel à ser_ldiv_sc_mat
    typedef ldiv_sc_mat_series_rep<C,V,W,matrix_naive> Div_sc_mat_rep;
    Series_vector s = (series_rep<Vector,V>*) Div_sc_mat_rep (smat0, v);
    return s; */
  }
};

TMPLWU static inline Vector_series
ser_ldiv_mat (const Matrix_series& mat, const Vector_series& f) {
  typedef ldiv_mat_series_rep<C,V,W,U> Div_mat_rep;
  if (is_exact_zero (as_series (f))) {
    Series zero (get_format1 (CF(mat)));
    return Vector_series (zero, cols(mat));
  }
  Series_vector s = (series_rep<Vector,V>*) new Div_mat_rep (mat,f);
  return as_vector (recursive (s)); 
}

#undef D

#undef Series_matrix
#undef Matrix_series
#undef Vector_series
#undef Series_vector
#undef Vector
#undef Matrix
#undef Series
#undef Series_rep
#undef TMPLWU
#undef TMPLU
#undef TMPL
} // namespace mmx
#endif // __MMX_SERIES_VECTOR_HPP