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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : series_sugar.hpp
* DESCRIPTION: Recursive and implicit power series via routines
* COPYRIGHT  : (C) 2009  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_SUGAR_HPP
#define __MMX_SERIES_SUGAR_HPP
#include <basix/routine.hpp>
#include <algebramix/series_implicit.hpp>
namespace mmx {
#define TMPL_DEF template<typename C, typename V=typename Series_variant(C)>
#define TMPL template<typename C, typename V>
#define Series series<C,V>
#define Series_rep series_rep<C,V>
#define VC vector<C>
#define VSeries series<vector<C> >
#define VSeries_rep series_rep<vector<C> >
#define UC unknown<C,V>
#define UC_rep unknown_rep<C,V>
#define USeries series<UC >
#define USeries_rep series_rep<UC >
#define Solver_rep solver_series_rep<C,V>

/******************************************************************************
* Fixed point series
******************************************************************************/

template<typename C>
class fixed_point_series_rep: public recursive_series_rep<C> {
  routine   fun;
  vector<C> c;
public:
  fixed_point_series_rep (const routine& fun2, const vector<C>& c2):
    recursive_series_rep<C> (CF(c2)), fun (fun2), c (c2) {}
  syntactic expression (const syntactic& z) const {
    (void) z;
    return apply ("fixed_point", flatten (fun), flatten (c)); }
  series<C> initialize () {
    // mmout << "Initialize " << c << "\n";
    for (nat i=0; i<N(c); i++) this->initial(i)= c[i];
    return as<series<C> > (fun->apply (as<generic> (this->me ()))); }
};

template<typename C> inline series<C>
fixed_point_series (const routine& fun, const vector<C>& c) {
  series_rep<C>* rep= new fixed_point_series_rep<C> (fun, c);
  return recursive (series<C> (rep));
}

template<typename C> inline series<C>
fixed_point_series (const routine& fun, const C& c) {
  return fixed_point_series (fun, vec<C> (c));
}

template<typename C> inline series<C>
integrate_series (const routine& fun, const C& c) {
  return fixed_point_series (integrate (fun), vec<C> (c));
}

template<typename C>
class fixed_point_vector_series_rep: public recursive_series_rep<vector<C> > {
  routine            fun;
  vector<vector<C> > c;
public:
  fixed_point_vector_series_rep (const routine& fun2,
                                 const vector<vector<C> >& c2):
    recursive_series_rep<vector<C> > (CF(c2)), fun (fun2), c (c2) {}
  syntactic expression (const syntactic& z) const {
    (void) z;
    return apply ("fixed_point", flatten (fun), flatten (c)); }
  series<vector<C> > initialize () {
    for (nat i=0; i<N(c); i++) this->initial(i)= c[i];
    vector<series<C> > in= as_vector (this->me ());
    generic genout= fun->apply (as<generic> (as<vector<generic> > (in)));
    vector<series<C> > out=
      as<vector<series<C> > > (as<vector<generic> > (genout));
    return from_vector (out); }
};

template<typename C> inline series<vector<C> >
fixed_point_series_vector (const routine& fun, const vector<vector<C> >& c) {
  series_rep<vector<C> >* rep= new fixed_point_vector_series_rep<C> (fun, c);
  return recursive (series<vector<C> > (rep));
}

template<typename C> inline vector<series<C> >
fixed_point_vector_series (const routine& fun, const vector<vector<C> >& c) {
  return as_vector (fixed_point_series_vector (fun, c));
}

template<typename C> inline series<vector<C> >
fixed_point_series_vector (const routine& fun, const vector<C>& c) {
  return fixed_point_series_vector (fun, vec<vector<C> > (c));
}
template<typename C> inline vector<series<C> >
fixed_point_vector_series (const routine& fun, const vector<C>& c) {
  return fixed_point_vector_series (fun, vec<vector<C> > (c));
}

template<typename C> inline vector<generic>
gen_fixed_point_vector_series (const routine& fun, const vector<C>& c) {
  return as<vector<generic> > (fixed_point_vector_series (fun, c));
}

template<typename C> inline vector<generic>
gen_integrate_vector_series (const routine& fun, const vector<C>& c) {
  return gen_fixed_point_vector_series (integrate (fun), c);
}

/******************************************************************************
* Implicit series
******************************************************************************/

template<typename C> vector<vector<C> >
singleton_vector (const vector<C>& v) {
  nat n= N(v);
  vector<vector<C> > r= fill<vector<C> > (n);
  for (nat i=0; i<n; i++) r[i]= vec<C> (v[i]);
  return r;
}

TMPL_DEF
class implicit_series_rep: public Solver_rep {
  routine   fun;
  vector<C> c;
public:
  implicit_series_rep (const routine& fun2, const vector<C>& c2):
    Solver_rep (1, singleton_vector (c2)), fun (fun2), c (c2) {}
  syntactic expression (const syntactic& z) const {
    (void) z;
    return apply ("implicit", flatten (fun), flatten (c)); }
  vector<USeries > initialize () {
    USeries in = this->me () [0];
    USeries out= as<USeries > (fun->apply (as<generic> (in)));
    return vec<USeries > (out); }
};

template<typename C> inline series<C>
implicit_series (const routine& fun, const vector<C>& c) {
  series_rep<VC >* rep= (series_rep<VC >*) new implicit_series_rep<C> (fun, c);
  return access (solver (series<VC > (rep)), 0);
}

template<typename C> inline series<C>
implicit_series (const routine& fun, const C& c) {
  return implicit_series (fun, vec<C> (c));
}

TMPL_DEF
class implicit_vector_series_rep: public Solver_rep {
  routine fun;
  vector<vector<C> > c;
public:
  implicit_vector_series_rep (const routine& fun2,
                              const vector<vector<C> >& c2):
    Solver_rep (N(c2[0]), c2), fun (fun2), c (c2) {}
  syntactic expression (const syntactic& z) const {
    (void) z;
    return apply ("implicit", flatten (fun), flatten (c)); }
  vector<USeries > initialize () {
    /*
    vector<USeries> mmm= this->me ();
    for (nat i=0; i<N(mmm); i++)
      for (nat j=0; j<3; j++)
        mmerr << i << ", " << j << " -> " << mmm[i][j] << "\n";
    */
    vector<generic> in =
      as<vector<generic> > (this->me ());
    vector<generic> out=
      as<vector<generic> > (fun->apply (as<generic> (in)));
    /*
    vector<USeries> rrr= as<vector<USeries > > (out);
    for (nat i=0; i<N(rrr); i++)
      for (nat j=0; j<3; j++)
        mmerr << i << ", " << j << " -> " << rrr[i][j] << "\n";    
    */
    return as<vector<USeries > > (out); }
};

template<typename C> inline vector<series<C> >
implicit_vector_series (const routine& fun, const vector<vector<C> >& c) {
  ASSERT (N(c) > 0, "at least one initial condition required");
  series_rep<VC >* rep=
    (series_rep<VC >*) new implicit_vector_series_rep<C> (fun, c);
  return as_vector (solver (series<VC > (rep)));
}

template<typename C> inline vector<series<C> >
implicit_vector_series (const routine& fun, const vector<C>& c) {
  return implicit_vector_series (fun, vec<vector<C> > (c));
}

template<typename C> inline vector<generic>
gen_implicit_vector_series (const routine& fun, const vector<C>& c) {
  return as<vector<generic> > (implicit_vector_series (fun, c));
}

#undef TMPL_DEF
#undef TMPL
#undef Series
#undef Series_rep
#undef VC
#undef VSeries
#undef VSeries_rep
#undef UC
#undef UC_rep
#undef Iseries
#undef Iseries_rep
#undef Solver_rep
} // namespace mmx
#endif // __MMX_SERIES_IMPLICIT_HPP