algebramix: include/algebramix/series

00001 
00002 /******************************************************************************
00003 * MODULE     : series_sugar.hpp
00004 * DESCRIPTION: Recursive and implicit power series via routines
00005 * COPYRIGHT  : (C) 2009  Joris van der Hoeven
00006 *******************************************************************************
00007 * This software falls under the GNU general public license and comes WITHOUT
00008 * ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
00009 * If you don't have this file, write to the Free Software Foundation, Inc.,
00010 * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00011 ******************************************************************************/
00012 
00013 #ifndef __MMX_SERIES_SUGAR_HPP
00014 #define __MMX_SERIES_SUGAR_HPP
00015 #include <basix/routine.hpp>
00016 #include <algebramix/series_implicit.hpp>
00017 namespace mmx {
00018 #define TMPL_DEF template<typename C, typename V=typename Series_variant(C)>
00019 #define TMPL template<typename C, typename V>
00020 #define Series series<C,V>
00021 #define Series_rep series_rep<C,V>
00022 #define VC vector<C>
00023 #define VSeries series<vector<C> >
00024 #define VSeries_rep series_rep<vector<C> >
00025 #define UC unknown<C,V>
00026 #define UC_rep unknown_rep<C,V>
00027 #define USeries series<UC >
00028 #define USeries_rep series_rep<UC >
00029 #define Solver_rep solver_series_rep<C,V>
00030 
00031 /******************************************************************************
00032 * Fixed point series
00033 ******************************************************************************/
00034 
00035 template<typename C>
00036 class fixed_point_series_rep: public recursive_series_rep<C> {
00037   routine   fun;
00038   vector<C> c;
00039 public:
00040   fixed_point_series_rep (const routine& fun2, const vector<C>& c2):
00041     recursive_series_rep<C> (CF(c2)), fun (fun2), c (c2) {}
00042   syntactic expression (const syntactic& z) const {
00043     (void) z;
00044     return apply ("fixed_point", flatten (fun), flatten (c)); }
00045   series<C> initialize () {
00046     // mmout << "Initialize " << c << "\n";
00047     for (nat i=0; i<N(c); i++) this->initial(i)= c[i];
00048     return as<series<C> > (fun->apply (as<generic> (this->me ()))); }
00049 };
00050 
00051 template<typename C> inline series<C>
00052 fixed_point_series (const routine& fun, const vector<C>& c) {
00053   series_rep<C>* rep= new fixed_point_series_rep<C> (fun, c);
00054   return recursive (series<C> (rep));
00055 }
00056 
00057 template<typename C> inline series<C>
00058 fixed_point_series (const routine& fun, const C& c) {
00059   return fixed_point_series (fun, vec<C> (c));
00060 }
00061 
00062 template<typename C> inline series<C>
00063 integrate_series (const routine& fun, const C& c) {
00064   return fixed_point_series (integrate (fun), vec<C> (c));
00065 }
00066 
00067 template<typename C>
00068 class fixed_point_vector_series_rep: public recursive_series_rep<vector<C> > {
00069   routine            fun;
00070   vector<vector<C> > c;
00071 public:
00072   fixed_point_vector_series_rep (const routine& fun2,
00073                                  const vector<vector<C> >& c2):
00074     recursive_series_rep<vector<C> > (CF(c2)), fun (fun2), c (c2) {}
00075   syntactic expression (const syntactic& z) const {
00076     (void) z;
00077     return apply ("fixed_point", flatten (fun), flatten (c)); }
00078   series<vector<C> > initialize () {
00079     for (nat i=0; i<N(c); i++) this->initial(i)= c[i];
00080     vector<series<C> > in= as_vector (this->me ());
00081     generic genout= fun->apply (as<generic> (as<vector<generic> > (in)));
00082     vector<series<C> > out=
00083       as<vector<series<C> > > (as<vector<generic> > (genout));
00084     return from_vector (out); }
00085 };
00086 
00087 template<typename C> inline series<vector<C> >
00088 fixed_point_series_vector (const routine& fun, const vector<vector<C> >& c) {
00089   series_rep<vector<C> >* rep= new fixed_point_vector_series_rep<C> (fun, c);
00090   return recursive (series<vector<C> > (rep));
00091 }
00092 
00093 template<typename C> inline vector<series<C> >
00094 fixed_point_vector_series (const routine& fun, const vector<vector<C> >& c) {
00095   return as_vector (fixed_point_series_vector (fun, c));
00096 }
00097 
00098 template<typename C> inline series<vector<C> >
00099 fixed_point_series_vector (const routine& fun, const vector<C>& c) {
00100   return fixed_point_series_vector (fun, vec<vector<C> > (c));
00101 }
00102 template<typename C> inline vector<series<C> >
00103 fixed_point_vector_series (const routine& fun, const vector<C>& c) {
00104   return fixed_point_vector_series (fun, vec<vector<C> > (c));
00105 }
00106 
00107 template<typename C> inline vector<generic>
00108 gen_fixed_point_vector_series (const routine& fun, const vector<C>& c) {
00109   return as<vector<generic> > (fixed_point_vector_series (fun, c));
00110 }
00111 
00112 template<typename C> inline vector<generic>
00113 gen_integrate_vector_series (const routine& fun, const vector<C>& c) {
00114   return gen_fixed_point_vector_series (integrate (fun), c);
00115 }
00116 
00117 /******************************************************************************
00118 * Implicit series
00119 ******************************************************************************/
00120 
00121 template<typename C> vector<vector<C> >
00122 singleton_vector (const vector<C>& v) {
00123   nat n= N(v);
00124   vector<vector<C> > r= fill<vector<C> > (n);
00125   for (nat i=0; i<n; i++) r[i]= vec<C> (v[i]);
00126   return r;
00127 }
00128 
00129 TMPL_DEF
00130 class implicit_series_rep: public Solver_rep {
00131   routine   fun;
00132   vector<C> c;
00133 public:
00134   implicit_series_rep (const routine& fun2, const vector<C>& c2):
00135     Solver_rep (1, singleton_vector (c2)), fun (fun2), c (c2) {}
00136   syntactic expression (const syntactic& z) const {
00137     (void) z;
00138     return apply ("implicit", flatten (fun), flatten (c)); }
00139   vector<USeries > initialize () {
00140     USeries in = this->me () [0];
00141     USeries out= as<USeries > (fun->apply (as<generic> (in)));
00142     return vec<USeries > (out); }
00143 };
00144 
00145 template<typename C> inline series<C>
00146 implicit_series (const routine& fun, const vector<C>& c) {
00147   series_rep<VC >* rep= (series_rep<VC >*) new implicit_series_rep<C> (fun, c);
00148   return access (solver (series<VC > (rep)), 0);
00149 }
00150 
00151 template<typename C> inline series<C>
00152 implicit_series (const routine& fun, const C& c) {
00153   return implicit_series (fun, vec<C> (c));
00154 }
00155 
00156 TMPL_DEF
00157 class implicit_vector_series_rep: public Solver_rep {
00158   routine fun;
00159   vector<vector<C> > c;
00160 public:
00161   implicit_vector_series_rep (const routine& fun2,
00162                               const vector<vector<C> >& c2):
00163     Solver_rep (N(c2[0]), c2), fun (fun2), c (c2) {}
00164   syntactic expression (const syntactic& z) const {
00165     (void) z;
00166     return apply ("implicit", flatten (fun), flatten (c)); }
00167   vector<USeries > initialize () {
00168     /*
00169     vector<USeries> mmm= this->me ();
00170     for (nat i=0; i<N(mmm); i++)
00171       for (nat j=0; j<3; j++)
00172         mmerr << i << ", " << j << " -> " << mmm[i][j] << "\n";
00173     */
00174     vector<generic> in =
00175       as<vector<generic> > (this->me ());
00176     vector<generic> out=
00177       as<vector<generic> > (fun->apply (as<generic> (in)));
00178     /*
00179     vector<USeries> rrr= as<vector<USeries > > (out);
00180     for (nat i=0; i<N(rrr); i++)
00181       for (nat j=0; j<3; j++)
00182         mmerr << i << ", " << j << " -> " << rrr[i][j] << "\n";    
00183     */
00184     return as<vector<USeries > > (out); }
00185 };
00186 
00187 template<typename C> inline vector<series<C> >
00188 implicit_vector_series (const routine& fun, const vector<vector<C> >& c) {
00189   ASSERT (N(c) > 0, "at least one initial condition required");
00190   series_rep<VC >* rep=
00191     (series_rep<VC >*) new implicit_vector_series_rep<C> (fun, c);
00192   return as_vector (solver (series<VC > (rep)));
00193 }
00194 
00195 template<typename C> inline vector<series<C> >
00196 implicit_vector_series (const routine& fun, const vector<C>& c) {
00197   return implicit_vector_series (fun, vec<vector<C> > (c));
00198 }
00199 
00200 template<typename C> inline vector<generic>
00201 gen_implicit_vector_series (const routine& fun, const vector<C>& c) {
00202   return as<vector<generic> > (implicit_vector_series (fun, c));
00203 }
00204 
00205 #undef TMPL_DEF
00206 #undef TMPL
00207 #undef Series
00208 #undef Series_rep
00209 #undef VC
00210 #undef VSeries
00211 #undef VSeries_rep
00212 #undef UC
00213 #undef UC_rep
00214 #undef Iseries
00215 #undef Iseries_rep
00216 #undef Solver_rep
00217 } // namespace mmx
00218 #endif // __MMX_SERIES_IMPLICIT_HPP
include/algebramix/series_sugar.hpp
