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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : series_matrix.hpp
* DESCRIPTION: Vectors and matrices of 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_MATRIX_HPP
#define __MMX_SERIES_MATRIX_HPP
#include <basix/vector.hpp>
#include <algebramix/matrix.hpp>
#include <algebramix/series.hpp>
#include <algebramix/series_vector.hpp>

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

/******************************************************************************
* Variant
******************************************************************************/

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

/******************************************************************************
* Matrices of series
******************************************************************************/

STYPE_TO_TYPE(TMPL,as_matrix_type,Series_matrix,Matrix_series);

TMPL
class matrix_access_series_rep: public Series_rep {
  Series_matrix f;
  nat i, j;
public:
  matrix_access_series_rep (const Series_matrix& f2, nat i2, nat j2):
    Series_rep (get_format1 (CF(f2))), f (f2), i (i2), j (j2) {}
  syntactic expression (const syntactic& z) const {
    return access (flatten (f, z), flatten (i), flatten (j)); }
  void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (f, l); }
  C next () { return f[this->n](i,j); }
};

TMPL Series
access (const Series_matrix& f, nat i, nat j) {
  return (Series_rep*) new matrix_access_series_rep<C,V,U> (f, i, j);
}

TMPL Matrix_series
as_matrix (const Series_matrix& f) {
  nat nr= rows (f[0]), nc= cols (f[0]);
  Matrix_series m (Series (0), nr, nc);
  for (nat i=0; i<nr; i++)
    for (nat j=0; j<nc; j++)
      m (i, j)= access (f, i, j);
  return m;
}

TMPL syntactic
flatten (const Matrix_series& m, const syntactic& z) {
  syntactic g;
  nat i, j, nr= rows (m), nc= cols (m);
  matrix<syntactic,U> r (g, nr, nc);
  for (i=0; i<nr; i++)
    for (j=0; j<nc; j++)
      r (i, j)= flatten (m (i, j), z);
  return flatten (r);
}

TMPL format<Matrix >
get_matrix_format (const Matrix_series& m) {
  format<Series > fm1= CF(m);
  format<C> fm2= get_format1 (fm1);
  C zero= promote (0, fm2);
  Matrix r (zero, rows (m), cols (m));
  return get_format (r);
}

TMPL
class matrix_series_rep: public series_rep<Matrix,V> {
  const Matrix_series m;
public:
  matrix_series_rep (const Matrix_series& m2):
    series_rep<Matrix,V> (get_matrix_format (m2)), m (m2) {}
  syntactic expression (const syntactic& z) const {
    return flatten (m, z); }
  void Increase_order (nat l) {
    series_rep<Matrix,V>::Increase_order (l);
    for (nat i=0; i<rows(m); i++)
      for (nat j=0; j<cols(m); j++)
        increase_order (m (i, j), l); }
  Matrix next () {
    nat i, j, nr= rows (m), nc= cols (m);
    Matrix coeff (C (0), nr, nc);
    for (i=0; i<nr; i++)
      for (j=0; j<nc; j++)
        coeff (i, j)= m (i, j) [this->n];
    return coeff; }
};

TMPL Series_matrix
as_series (const Matrix_series& m) {
  return (series_rep<Matrix,V>*) new matrix_series_rep<C,V,U> (m);
}

TMPL Series_matrix
from_matrix (const Matrix_series& m) {
  return (series_rep<Matrix,V>*) new matrix_series_rep<C,V,U> (m);
}

TMPL
class lshiftz_series_matrix_rep: public series_rep<Matrix,V> {
  const Series_matrix f;
  const nat r,c;
  const int shift;
public:
  lshiftz_series_matrix_rep (const Series_matrix& f2, nat r2, nat c2,
                             int shift2):
    series_rep<Matrix,V> (CF(f2)),
    f (f2),  r (r2), c (c2), shift (shift2) {}
  syntactic expression (const syntactic& z) const {
    return lshiftz (flatten (f, z), flatten (shift)); }
  void Increase_order (nat l) {
    series_rep<Matrix,V>::Increase_order (l);
    increase_order (f, max (0, ((int) l) - shift)); }
  Matrix next () { return shift > ((int) this->n)? 
      Matrix (C (), r, c): f[this->n - shift]; }
};

TMPL inline Series_matrix
lshiftz_series_matrix (const Series_matrix& f, const nat& r, const nat& c,
                       const int& shift=1) {
  if (is_exact_zero (f)) {
    Series zero (get_format1 (CF(f)));
    return as_series (Matrix_series (zero, r, c));
  }
  return (series_rep<Matrix,V>*) 
    new lshiftz_series_matrix_rep<C,V,U> (f, r, c, shift);
}

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

TMPL Series_matrix
solve_matrix_lde_init (const Series_matrix& f, const Matrix& c) {
  return unary_recursive_series<solve_matrix_lde_op> (f, c);
}

TMPL Matrix_series
solve_lde (const Matrix_series& f) {
  ASSERT (is_square_matrix (f), "square matrix expected");
  Matrix c (C(1), rows (f), cols (f));
  return as_matrix (solve_matrix_lde_init (as_series (f), c));
}

TMPL Matrix_series
solve_lde_init (const Matrix_series& f, const Matrix& c) {
  ASSERT (is_square_matrix (f), "square matrix expected");
  ASSERT (rows (f) == rows (c), "unequal matrix dimensions");
  return as_matrix (solve_matrix_lde_init (as_series (f), c));
}



/******************************************************************************
* Matrix equations resolution
******************************************************************************/

#define D       typename C::C

template<typename C,typename V,typename U,typename UU>
static inline Matrix_series
matrix_mul_nocarry (const matrix<C,UU>& mat, const Matrix_series& m) {
  typedef typename Series_variant(D) SV;
  return as<Matrix_series> (as<matrix<series<C,SV>, U> > (mat) * 
                            as<matrix<series<C,SV>, U> > (m));
}

template<typename C,typename V,typename U,typename UU>
static Matrix_series
matrix_mul_quo (const matrix<C,UU>& Mat, const Matrix_series& ma) {
  Matrix mat= as<Matrix> (Mat);
  nat m= rows(mat), n= cols(mat);
  ASSERT (rows(ma) == n, "sizes don't match");
  if (is_a_scalar (mat)) return as<Series> (mat.scalar()) * ma;
  ASSERT (is_non_scalar (ma), "non-scalar matrix expected");
  Series zero (CF(Mat));
  Matrix_series w (zero, m, cols(ma));
  Series mr, mq, c, ar, aq;
  for (nat k=0; k < m; k++) 
    for (nat j=0; j < cols(ma); j++) {
      mr = carry_mul_rem_series (mat(k,0), ma(0,j), mq);
      c  = mq;
      ar = mr;
      for (nat l=1; l < n; l++) {
        if (is_exact_zero (mat (k,l))) {
          mr = zero;
          mq = zero;
          aq = zero;
          continue;
        }
        mr = carry_mul_rem_series (mat(k,l), ma(l,j), mq);
        ar = carry_add_rem (ar, mr, aq);
        c  = carry_special_add (c, mq, aq);  
      }
      w(k,j) = c;
    }
  return w;
}


template<typename C,typename V,typename U,typename UU>
class ldiv_sc_mat_mat_series_rep: public recursive_series_rep<Matrix,V> {
protected:
  Matrix_series f;
  Matrix c;
public:
  ldiv_sc_mat_mat_series_rep (const matrix<C,UU>& c2, const Matrix_series& f2):
    recursive_series_rep<Matrix,V> (get_matrix_format (f2)),
    f (f2), c (as<Matrix> (c2)) {}
  syntactic expression (const syntactic& z) const {
    return invert (flatten (c)) * flatten (f,z); }
  virtual void Increase_order (nat l) {
    recursive_series_rep<Matrix,V>::Increase_order (l);
    for (nat i = 0; i < rows(f); i++)
      for (nat j = 0; j < cols(f); j++)
        increase_order (f(i,j),l);
  }
  Series_matrix initialize () {
    Matrix u= invert (c);
    this->initial (0)= u * as_series (f)[0];
    Series_matrix ls_me = lshiftz_series_matrix (this -> me (),
                                                 cols(c), cols(f));
    Matrix_series tmp= f - matrix_mul_quo (c, as_matrix (ls_me));
    return as_series (matrix_mul_nocarry (u, tmp));
  }
};

template<typename C,typename V,typename U,typename UU> 
static inline Matrix_series
ser_ldiv_sc_mat_mat (const matrix<C,UU>& c, const Matrix_series& f) {
  typedef ldiv_sc_mat_mat_series_rep<C,V,U,UU> Div_sc_mat_mat_rep;
  if (is_exact_zero (as_series (f))) {
    Series zero (get_format1 (CF(f)));
    return Matrix_series (zero, cols(c), cols(f));
  }
  Series_matrix s = (series_rep<Matrix,V>*) new Div_sc_mat_mat_rep (c,f);
  return as_matrix (recursive (s));
}

TMPL
class ldiv_mat_mat_series_rep: public recursive_series_rep<Matrix,V> {
protected:
  Matrix_series f;
  Matrix_series mat;
public:
  ldiv_mat_mat_series_rep (const Matrix_series& mat2, const Matrix_series& f2):
    recursive_series_rep<Matrix,V> (get_matrix_format (mat2)),
    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<Matrix,V>::Increase_order (l);
    for (nat i = 0; i < rows(f); i++) {
      for (nat j = 0; j < cols(f); j++)
        increase_order (f(i,j),l);
      for (nat j = 0; j < cols(mat); j++)
        increase_order (mat(i,j),l);
    }
  }
  Series_matrix initialize () {
    Series_matrix smat = as_series (mat);
    //matrix<C> smat0 = smat[0]; 
    //this->initial (0)= invert (smat0) * ((matrix<C>) as_series (f)[0]);
    this->initial (0)= invert (smat[0]) * as_series (f)[0];
    Matrix_series ma = f - 
      as_matrix (lshiftz_series_matrix (as_series (as_matrix (rshiftz (smat))*
                                                   as_matrix (this-> me ())),
                                        cols(mat), cols(f)));
    return as_series (ser_ldiv_sc_mat_mat (smat[0], ma));
  }
};

TMPL static inline Matrix_series
ser_ldiv_mat_mat (const Matrix_series& mat, const Matrix_series& f) {
  typedef ldiv_mat_mat_series_rep<C,V,U> Div_mat_rep;
  if (is_exact_zero (as_series (f))) {
    Series zero (get_format1 (CF(f)));
    return Matrix_series (zero, cols(mat), cols(f));
  }
  Series_matrix s = (series_rep<Matrix,V>*) new Div_mat_rep (mat,f);
  return as_matrix (recursive (s)); 
}


#undef D

#undef Series_matrix
#undef Matrix_series
#undef Matrix
#undef Series
#undef Series_rep
#undef TMPLWU
#undef TMPL
} // namespace mmx
#endif // __MMX_SERIES_MATRIX_HPP