SyFi 0.3
ginac_tools.cpp
Go to the documentation of this file.
00001 // Copyright (C) 2006-2009 Kent-Andre Mardal and Simula Research Laboratory.
00002 // Licensed under the GNU GPL Version 2, or (at your option) any later version.
00003 
00004 #include "ginac_tools.h"
00005 #include "symbol_factory.h"
00006 
00007 #include <ginac/ginac.h>
00008 
00009 #include <iostream>
00010 #include <sstream>
00011 #include <fstream>
00012 #include <stdexcept>
00013 
00014 using namespace std;
00015 using namespace GiNaC;
00016 
00017 namespace SyFi
00018 {
00019 
00020         GiNaC::lst cross(GiNaC::lst& v1, GiNaC::lst& v2)
00021         {
00022                 GiNaC::lst ret;
00023                 if ( v1.nops() != v2.nops() )
00024                 {
00025                         cout <<"incompatible vectors "<<endl;
00026                         cout <<"v1.nops() "<<v1.nops();
00027                         cout <<"  v2.nops() "<<v2.nops()<<endl; ;
00028                         return GiNaC::lst();
00029                 }
00030                 ret.append(  v1.op(1)*v2.op(2) - v1.op(2)*v2.op(1));
00031                 ret.append(- v1.op(0)*v2.op(2) + v1.op(2)*v2.op(0));
00032                 ret.append(  v1.op(0)*v2.op(1) - v1.op(1)*v2.op(0));
00033                 return ret;
00034         }
00035 
00036         GiNaC::ex inner(GiNaC::ex a, GiNaC::ex b, bool transposed)
00037         {
00038                 if (GiNaC::is_a<GiNaC::matrix>(a) && GiNaC::is_a<GiNaC::matrix>(b))
00039                 {
00040                         GiNaC::matrix ma = GiNaC::ex_to<GiNaC::matrix>(a);
00041                         GiNaC::matrix mb = GiNaC::ex_to<GiNaC::matrix>(b);
00042                         if ( !transposed )
00043                         {
00044                                 if (ma.cols() != mb.cols() || ma.rows() != mb.rows() )
00045                                 {
00046                                         cout <<"Incompatible matrices "<<endl;
00047                                         cout <<"a.cols() "<<ma.cols()<<endl;
00048                                         cout <<"a.rows() "<<ma.rows()<<endl;
00049                                         cout <<"b.cols() "<<mb.cols()<<endl;
00050                                         cout <<"b.rows() "<<mb.rows()<<endl;
00051                                         cout <<"a="<<a<<endl;
00052                                         cout <<"b="<<b<<endl;
00053                                         throw std::runtime_error("Incompatible matrices.");
00054                                 }
00055 
00056                                 GiNaC::ex ret;
00057                                 for (unsigned int i=0; i<ma.rows(); i++)
00058                                 {
00059                                         for (unsigned int j=0; j<ma.cols(); j++)
00060                                         {
00061                                                 ret += ma(i,j)*mb(i,j);
00062                                         }
00063                                 }
00064                                 return ret;
00065                         }
00066                         else
00067                         {
00068                                 if (ma.cols() != mb.rows() || ma.rows() != mb.cols() )
00069                                 {
00070                                         cout <<"Incompatible matrices "<<endl;
00071                                         cout <<"a.cols() "<<ma.cols()<<endl;
00072                                         cout <<"a.rows() "<<ma.rows()<<endl;
00073                                         cout <<"b.cols() "<<mb.cols()<<endl;
00074                                         cout <<"b.rows() "<<mb.rows()<<endl;
00075                                         cout <<"a="<<a<<endl;
00076                                         cout <<"b="<<b<<endl;
00077                                         throw std::runtime_error("Incompatible matrices.");
00078                                 }
00079 
00080                                 GiNaC::ex ret;
00081                                 for (unsigned int i=0; i<ma.rows(); i++)
00082                                 {
00083                                         for (unsigned int j=0; j<ma.cols(); j++)
00084                                         {
00085                                                 ret += ma(i,j)*mb(j,i);
00086                                         }
00087                                 }
00088                                 return ret;
00089                         }
00090                 }
00091                 else if (GiNaC::is_a<GiNaC::lst>(a)
00092                         && GiNaC::is_a<GiNaC::lst>(b))
00093                 {
00094                         return inner(GiNaC::ex_to<GiNaC::lst>(a), GiNaC::ex_to<GiNaC::lst>(b));
00095                 }
00096                 else
00097                 {
00098                         return a*b;
00099                 }
00100         }
00101 
00102         GiNaC::ex inner(GiNaC::lst v1, GiNaC::lst v2)
00103         {
00104                 GiNaC::ex ret;
00105 
00106                 if ( v1.nops() != v2.nops() )
00107                 {
00108                         cout <<"incompatible vectors "<<endl;
00109                         cout <<"v1.nops() "<<v1.nops();
00110                         cout <<"  v2.nops() "<<v2.nops()<<endl; ;
00111                         return 0;
00112                 }
00113                 for (unsigned i = 0; i <= v1.nops()-1 ; ++i)
00114                 {
00115                         if ( GiNaC::is_a<GiNaC::lst>(v1.op(i)) &&
00116                                 GiNaC::is_a<GiNaC::lst>(v2.op(i)) )
00117                         {
00118                                 ret += inner(GiNaC::ex_to<GiNaC::lst>(v1.op(i)),
00119                                         GiNaC::ex_to<GiNaC::lst>(v2.op(i)));
00120                         }
00121                         else
00122                         {
00123                                 ret += v1.op(i)*v2.op(i);
00124                         }
00125                 }
00126                 return ret;
00127         }
00128 
00129         GiNaC::ex inner(GiNaC::exvector& v1, GiNaC::exvector& v2)
00130         {
00131                 GiNaC::ex ret;
00132                 for (unsigned int i=0; i< v1.size(); i++)
00133                 {
00134                         ret += v1[i]*v2[i];
00135                 }
00136                 return ret;
00137         }
00138 
00139         GiNaC::lst matvec(GiNaC::matrix& M, GiNaC::lst& x)
00140         {
00141                 GiNaC::lst ret;
00142                 int nr = M.rows();
00143                 int nc = M.cols();
00144                 for (int i = 0; i < nr; i++)
00145                 {
00146                         GiNaC::ex tmp;
00147                         for (int j = 0; j < nc; j++)
00148                         {
00149                                 tmp = tmp +  M(i,j)*(x.op(j));
00150                         }
00151                         ret.append(tmp);
00152                 }
00153                 return ret;
00154         }
00155 
00156         GiNaC::ex matvec(GiNaC::ex A, GiNaC::ex x)
00157         {
00158                 ex sol;
00159 
00160                 if (GiNaC::is_a<GiNaC::matrix>(A) && GiNaC::is_a<GiNaC::matrix>(x))
00161                 {
00162                         GiNaC::matrix AA = GiNaC::ex_to<GiNaC::matrix>(A);
00163                         GiNaC::matrix xx = GiNaC::ex_to<GiNaC::matrix>(x);
00164                         sol = AA.mul(xx);
00165                 }
00166                 else
00167                 {
00168                         throw std::runtime_error("Invalid argument types, need matrices");
00169                 }
00170                 return sol;
00171         }
00172 
00173         GiNaC::lst ex2equations(GiNaC::ex rel)
00174         {
00175                 GiNaC::ex lhs = rel.lhs();
00176                 GiNaC::ex rhs = rel.rhs();
00177 
00178                 GiNaC::ex l;
00179                 GiNaC::ex r;
00180 
00181                 GiNaC::lst eqs;
00182 
00183                 for (int i=lhs.ldegree(x); i<=lhs.degree(x); ++i)
00184                 {
00185                         for (int j=lhs.ldegree(y); j<=lhs.degree(y); ++j)
00186                         {
00187                                 for (int k=lhs.ldegree(z); k<=lhs.degree(z); ++k)
00188                                 {
00189                                         l = lhs.coeff(x,i).coeff(y, j).coeff(z,k);
00190                                         r = rhs.coeff(x,i).coeff(y, j).coeff(z,k);
00191                                         //      if (! (l == 0 && r == 0 ) )  eqs.append(l == r); OLD VERSION
00192                                         if ( (l != 0 && (r == 0 || r == 1) ) )  eqs.append(l == r);
00193                                 }
00194                         }
00195                 }
00196                 eqs.sort();
00197                 return eqs;
00198         }
00199 
00200         GiNaC::lst collapse(GiNaC::lst l)
00201         {
00202                 GiNaC::lst lc;
00203                 GiNaC::lst::const_iterator iter1, iter2;
00204 
00205                 for (iter1 = l.begin(); iter1 != l.end(); ++iter1)
00206                 {
00207                         if (GiNaC::is_a<GiNaC::lst>(*iter1))
00208                         {
00209                                 for (iter2 = GiNaC::ex_to<GiNaC::lst>(*iter1).begin(); iter2 != GiNaC::ex_to<GiNaC::lst>(*iter1).end(); ++iter2)
00210                                 {
00211                                         lc.append(*iter2);
00212                                 }
00213                         }
00214                         else
00215                         {
00216                                 lc.append(*iter1);
00217                         }
00218                 }
00219                 lc.sort();
00220                 lc.unique();
00221                 return lc;
00222         }
00223 
00224         GiNaC::matrix equations2matrix(const GiNaC::ex &eqns, const GiNaC::ex &symbols)
00225         {
00226 
00227                 GiNaC::matrix sys(eqns.nops(),symbols.nops());
00228                 GiNaC::matrix rhs(eqns.nops(),1);
00229                 GiNaC::matrix vars(symbols.nops(),1);
00230 
00231                 for (size_t r=0; r<eqns.nops(); r++)
00232                 {
00233                                                                  // lhs-rhs==0
00234                         const GiNaC::ex eq = eqns.op(r).op(0)-eqns.op(r).op(1);
00235                         GiNaC::ex linpart = eq;
00236                         for (size_t c=0; c<symbols.nops(); c++)
00237                         {
00238                                 const GiNaC::ex co = eq.coeff(GiNaC::ex_to<GiNaC::symbol>(symbols.op(c)),1);
00239                                 linpart -= co*symbols.op(c);
00240                                 sys(r,c) = co;
00241                         }
00242                         linpart = linpart.expand();
00243                         rhs(r,0) = -linpart;
00244                 }
00245                 return sys;
00246         }
00247 
00248         void matrix_from_equations(const GiNaC::ex &eqns, const GiNaC::ex &symbols, GiNaC::matrix& A, GiNaC::matrix& b)
00249         {
00250                 // build matrix from equation system
00251                 GiNaC::matrix sys(eqns.nops(),symbols.nops());
00252                 GiNaC::matrix rhs(eqns.nops(),1);
00253                 GiNaC::matrix vars(symbols.nops(),1);
00254 
00255                 for (size_t r=0; r<eqns.nops(); r++)
00256                 {
00257                                                                  // lhs-rhs==0
00258                         const GiNaC::ex eq = eqns.op(r).op(0)-eqns.op(r).op(1);
00259                         GiNaC::ex linpart = eq;
00260                         for (size_t c=0; c<symbols.nops(); c++)
00261                         {
00262                                 const GiNaC::ex co = eq.coeff(GiNaC::ex_to<GiNaC::symbol>(symbols.op(c)),1);
00263                                 linpart -= co*symbols.op(c);
00264                                 sys(r,c) = co;
00265                         }
00266                         linpart = linpart.expand();
00267                         rhs(r,0) = -linpart;
00268                 }
00269                 A = sys;
00270                 b = rhs;
00271         }
00272 
00273         GiNaC::ex lst_to_matrix2(const GiNaC::lst& l)
00274         {
00275                 GiNaC::lst::const_iterator itr, itc;
00276 
00277                 // Find number of rows and columns
00278                 size_t rows = l.nops(), cols = 0;
00279                 for (itr = l.begin(); itr != l.end(); ++itr)
00280                 {
00281                         if (!GiNaC::is_a<GiNaC::lst>(*itr))
00282                                 //              throw (std::invalid_argument("lst_to_matrix: argument must be a list of lists"));
00283                                 cols = 1;
00284                         if (itr->nops() > cols)
00285                                 cols = itr->nops();
00286                 }
00287                 // Allocate and fill matrix
00288                 GiNaC::matrix &M = *new GiNaC::matrix(rows, cols);
00289                 M.setflag(GiNaC::status_flags::dynallocated);
00290 
00291                 unsigned i;
00292                 for (itr = l.begin(), i = 0; itr != l.end(); ++itr, ++i)
00293                 {
00294                         unsigned j;
00295                         if (cols == 1)
00296                         {
00297                                 M(i, 0) = *itr;
00298                         }
00299                         else
00300                         {
00301                                 for (itc = GiNaC::ex_to<GiNaC::lst>(*itr).begin(), j = 0; itc != GiNaC::ex_to<GiNaC::lst>(*itr).end(); ++itc, ++j)
00302                                         M(i, j) = *itc;
00303                         }
00304                 }
00305                 return M;
00306         }
00307 
00308         GiNaC::lst matrix_to_lst2(const GiNaC::ex& m)
00309         {
00310                 if (GiNaC::is_a<GiNaC::matrix>(m))
00311                 {
00312                         GiNaC::matrix A = GiNaC::ex_to<GiNaC::matrix>(m);
00313                         int cols = A.cols();
00314                         int rows = A.rows();
00315 
00316                         GiNaC::lst ret;
00317                         if ( cols == 1 )
00318                         {
00319                                 for (unsigned int i=0; i<=A.rows()-1; i++)
00320                                 {
00321                                         ret.append(A(i,0));
00322                                 }
00323                         }
00324                         else if ( rows == 1 )
00325                         {
00326                                 for (unsigned int i=0; i<=A.cols()-1; i++)
00327                                 {
00328                                         ret.append(A(0,i));
00329                                 }
00330                         }
00331                         else
00332                         {
00333                                 for (unsigned int i=0; i<=A.rows()-1; i++)
00334                                 {
00335                                         GiNaC::lst rl;
00336                                         for (unsigned int j=0; j<=A.cols()-1; j++)
00337                                         {
00338                                                 rl.append(A(i,j));
00339                                         }
00340                                         ret.append(rl);
00341                                 }
00342                         }
00343                         return ret;
00344                 }
00345                 else
00346                 {
00347                         return GiNaC::lst();
00348                 }
00349         }
00350 
00351         GiNaC::lst lst_equals(GiNaC::ex a, GiNaC::ex b)
00352         {
00353                 GiNaC::lst ret;
00354                 if ( (GiNaC::is_a<GiNaC::lst>(a)) && (GiNaC::is_a<GiNaC::lst>(b)) /*&& (a.nops() == b.nops())*/ ) {
00355                 for (unsigned int i=0; i<= a.nops()-1; i++)
00356                 {
00357                         ret.append(b.op(i) == a.op(i));
00358                 }
00359         }
00360         else if ( !(GiNaC::is_a<GiNaC::lst>(a)) && !(GiNaC::is_a<GiNaC::lst>(b)))
00361         {
00362                 ret.append(b == a);
00363         }
00364         else if ( !(GiNaC::is_a<GiNaC::lst>(a)) && (GiNaC::is_a<GiNaC::lst>(b)))
00365         {
00366                 ret.append(b.op(0) == a);
00367         }
00368         else
00369         {
00370                 throw(std::invalid_argument("Make sure that the lists a and b are comparable."));
00371         }
00372         return ret;
00373 }
00374 
00375 
00376 int find(GiNaC::ex e, GiNaC::lst list)
00377 {
00378         for (unsigned int i=0; i< list.nops(); i++)
00379         {
00380                 if ( e == list.op(i) ) return i;
00381         }
00382         return -1;
00383 }
00384 
00385 
00386 void visitor_subst_pow(GiNaC::ex e, GiNaC::exmap& map, ex_int_map& intmap, string a)
00387 {
00388         static int i=0;
00389         if (map.find(e) != map.end())
00390         {
00391                 intmap[e] = intmap[e]+1;
00392                 return;
00393         }
00394         if (GiNaC::is_exactly_a<GiNaC::power>(e))
00395         {
00396                 std::ostringstream s;
00397                 s <<a<<i++;
00398                 map[e] = GiNaC::symbol(s.str());
00399                 intmap[e] = 0;
00400                 for (unsigned int i=0; i< e.nops(); i++)
00401                 {
00402                         GiNaC::ex e2 = e.op(i);
00403                         //       cout <<"power e "<<e2<<endl;
00404                         visitor_subst_pow(e2,map,intmap, a);
00405                 }
00406         }
00407         else if (GiNaC::is_a<GiNaC::function>(e))
00408         {
00409                 std::ostringstream s;
00410                 s <<a<<i++;
00411                 map[e] = GiNaC::symbol(s.str());
00412                 intmap[e] = 0;
00413                 for (unsigned int i=0; i< e.nops(); i++)
00414                 {
00415                         GiNaC::ex e2 = e.op(i);
00416                         //       cout <<"function e "<<e2<<endl;
00417                         visitor_subst_pow(e2,map,intmap, a);
00418                 }
00419         }
00420         else if (GiNaC::is_a<GiNaC::mul>(e))
00421         {
00422                 if (e.nops() > 4 && e.nops() < 10 )
00423                 {
00424                         std::ostringstream s;
00425                         s <<a<<i++;
00426                         map[e] = GiNaC::symbol(s.str());
00427                         intmap[e] = 0;
00428                 }
00429 
00430                 for (unsigned int i=0; i< e.nops(); i++)
00431                 {
00432                         GiNaC::ex e2 = e.op(i);
00433                         visitor_subst_pow(e2,map,intmap, a);
00434                 }
00435         }
00436         else if (GiNaC::is_a<GiNaC::add>(e))
00437         {
00438                 for (unsigned int i=0; i< e.nops(); i++)
00439                 {
00440                         GiNaC::ex e2 = e.op(i);
00441                         visitor_subst_pow(e2,map,intmap,a);
00442                 }
00443         }
00444 }
00445 
00446 
00447 void check_visitor(GiNaC::ex e, GiNaC::lst& exlist)
00448 {
00449         if (find(e, exlist) >= 0) return;
00450 
00451         //  cout <<"ex e "<<e<<endl;
00452         if (GiNaC::is_a<GiNaC::numeric>(e))
00453         {
00454         }
00455         else if (GiNaC::is_a<GiNaC::add>(e) )
00456         {
00457                 //    cout <<"e "<<e <<endl;
00458                 //    cout <<"e.nops() "<<e.nops() <<endl;
00459                 if (e.nops() > 4 && e.nops() < 10 ) exlist.append(e);
00460                 for (unsigned int i=0; i< e.nops(); i++)
00461                 {
00462                         GiNaC::ex e2 = e.op(i);
00463                         //       cout <<"add e "<<e2<<endl;
00464                         //       exlist.append(e2);
00465                         check_visitor(e2,exlist);
00466                 }
00467         }
00468         else if (GiNaC::is_a<GiNaC::mul>(e))
00469         {
00470                 for (unsigned int i=0; i< e.nops(); i++)
00471                 {
00472                         GiNaC::ex e2 = e.op(i);
00473                         //       cout <<"mul e "<<e2<<endl;
00474                         exlist.append(e2);
00475                         check_visitor(e2,exlist);
00476                 }
00477         }
00478         else if (GiNaC::is_a<GiNaC::lst>(e))
00479         {
00480                 for (unsigned int i=0; i< e.nops(); i++)
00481                 {
00482                         GiNaC::ex e2 = e.op(i);
00483                         //       cout <<"GiNaC::lst e "<<e2<<endl;
00484                         //       exlist.append(e2);
00485                         check_visitor(e2,exlist);
00486                 }
00487         }
00488         else if (GiNaC::is_exactly_a<GiNaC::power>(e))
00489         {
00490                 exlist.append(e);
00491                 for (unsigned int i=0; i< e.nops(); i++)
00492                 {
00493                         GiNaC::ex e2 = e.op(i);
00494                         //       cout <<"power e "<<e2<<endl;
00495                         check_visitor(e2,exlist);
00496                 }
00497         }
00498         else if (GiNaC::is_a<GiNaC::function>(e))
00499         {
00500                 exlist.append(e);
00501                 for (unsigned int i=0; i< e.nops(); i++)
00502                 {
00503                         GiNaC::ex e2 = e.op(i);
00504                         //       cout <<"function e "<<e2<<endl;
00505                         check_visitor(e2,exlist);
00506                 }
00507         }
00508 
00509         else
00510         {
00511                 //       exlist.append(e);
00512                 //    cout <<"atom e "<<e<<endl;
00513         }
00514 
00515         exlist.sort();
00516         exlist.unique();
00517 }
00518 
00519 
00520 GiNaC::ex homogenous_pol(unsigned int order, unsigned int nsd, const string a)
00521 {
00522         using SyFi::x;
00523         using SyFi::y;
00524         using SyFi::z;
00525 
00526         if ( nsd == 1)
00527         {
00528                 GiNaC::symbol a0(istr(a,0));
00529                 return GiNaC::lst(a0*pow(x,order), a0, pow(x,order));
00530         }
00531         else if ( nsd == 2 )
00532         {
00533                 GiNaC::ex variables = get_symbolic_matrix(1,order+1, a);
00534                 GiNaC::lst basis;
00535                 GiNaC::ex ret;
00536                 for (unsigned int i=0; i<= order; i++)
00537                 {
00538                         basis.append(pow(x,i)*pow(y,order-i));
00539                         ret += variables.op(i)*basis.op(i);
00540                 }
00541                 return GiNaC::lst(ret, matrix_to_lst2(variables), basis);
00542         }
00543         else if ( nsd == 3 )
00544         {
00545                 GiNaC::lst basis;
00546                 for (unsigned int i=0; i<= order; i++)
00547                 {
00548                         for (unsigned int j=0; j<= order; j++)
00549                         {
00550                                 for (unsigned int k=0; k<= order; k++)
00551                                 {
00552                                         if ( i + j + k == order )
00553                                         {
00554                                                 basis.append(pow(x,i)*pow(y,j)*pow(z,k));
00555                                         }
00556                                 }
00557                         }
00558                 }
00559                 GiNaC::ex variables = get_symbolic_matrix(1,basis.nops(), a);
00560                 GiNaC::ex ret;
00561                 for (unsigned int i=0; i<basis.nops(); i++)
00562                 {
00563                         ret += variables.op(i)*basis.op(i);
00564                 }
00565                 return GiNaC::lst(ret, matrix_to_lst2(variables), basis);
00566         }
00567         throw std::runtime_error("Homogenous polynomials only implemented in 1D, 2D and 3D");
00568 }
00569 
00570 
00571 GiNaC::lst homogenous_polv(unsigned int no_fields, unsigned int order, unsigned int nsd, const string a)
00572 {
00573         GiNaC::lst ret1;                         // contains the polynom
00574         GiNaC::lst ret2;                         // contains the coefficients
00575         GiNaC::lst ret3;                         // constains the basis functions
00576         GiNaC::lst basis_tmp;
00577         for (unsigned int i=0; i< no_fields; i++)
00578         {
00579                 GiNaC::lst basis;
00580                 std::ostringstream s;
00581                 s <<a<<""<<i<<"_";
00582                 GiNaC::ex polspace = homogenous_pol(order, nsd, s.str());
00583                 ret1.append(polspace.op(0));
00584                 ret2.append(polspace.op(1));
00585                 basis_tmp = GiNaC::ex_to<GiNaC::lst>(polspace.op(2));
00586                 for (GiNaC::lst::const_iterator basis_iterator = basis_tmp.begin();
00587                         basis_iterator != basis_tmp.end(); ++basis_iterator)
00588                 {
00589                         GiNaC::lst tmp_lst;
00590                         for (unsigned int d=1; d<=no_fields; d++) tmp_lst.append(0);
00591                         tmp_lst.let_op(i) = (*basis_iterator);
00592                         ret3.append(tmp_lst);
00593                 }
00594         }
00595         return GiNaC::lst(ret1,ret2,ret3);
00596 }
00597 
00598 
00599 GiNaC::ex pol(unsigned int order, unsigned int nsd, const string a)
00600 {
00601         using SyFi::x;
00602         using SyFi::y;
00603         using SyFi::z;
00604 
00605         GiNaC::ex ret;                           // ex to return
00606         int dof;                                         // degrees of freedom
00607         GiNaC::ex A;                             // ex holding the coefficients a_0 .. a_dof
00608         GiNaC::lst basis;
00609 
00610         if (nsd == 1)
00611         {
00612                 /* 1D:
00613                  * P^n = a_0 + a_1*x + .... + a_n*x^n
00614                  * dof : n+1
00615                  */
00616                 dof = order+1;
00617                 A = get_symbolic_matrix(1,dof, a);
00618                 int o=0;
00619                 for (GiNaC::const_iterator i = A.begin(); i != A.end(); ++i)
00620                 {
00621                         ret += (*i)*pow(x,o);
00622                         basis.append(pow(x,o));
00623                         o++;
00624                 }
00625         }
00626         else if ( nsd == 2)
00627         {
00628 
00629                 /* 2D: structure of coefficients (a_i)
00630                  * [ a_0      a_1 x     a_3 x^2     a_6 x^3
00631                  * [ a_2 y    a_4 xy    a_7 x^2y
00632                  * [ a_5 y^2  a_8 xy^2
00633                  * [ a_9 y^3
00634                  */
00635                 dof = (order+1)*(order+2)/2;
00636                 A = get_symbolic_matrix(1, dof , a);
00637 
00638                 size_t i=0;
00639                 for (unsigned int o = 0; o <= order; o++)
00640                 {
00641                         for (unsigned int d = 0; d <= o; d++)
00642                         {
00643                                 ret += A.op(i)*pow(y,d)*pow(x,o-d);
00644                                 basis.append(pow(y,d)*pow(x,o-d));
00645                                 i++;
00646                         }
00647                 }
00648         }
00649         else if (nsd == 3)
00650         {
00651 
00652                 /* Similar structure as in 2D, but
00653                  * structured as a tetraheder, i.e.,
00654                  *   a_o + a_1 x + a_2 y + a_3 z
00655                  * + a_4 x^2 + a_5 xy +
00656                  */
00657                 dof = 0;
00658                 for (unsigned int j=0; j<= order; j++)
00659                 {
00660                         dof += (j+1)*(j+2)/2;
00661                 }
00662                 A = get_symbolic_matrix(1, dof , a);
00663 
00664                 size_t i=0;
00665                 for (unsigned int o = 0; o <= order; o++)
00666                 {
00667                         for (unsigned int d = 0; d <= o; d++)
00668                         {
00669                                 for (unsigned int f = 0; f <= o; f++)
00670                                 {
00671                                         if ( int(o)-int(d)-int(f) >= 0)
00672                                         {
00673                                                 ret += A.op(i)*pow(y,f)*pow(z,d)*pow(x,o-d-f);
00674                                                 basis.append(pow(y,f)*pow(z,d)*pow(x,o-d-f));
00675                                                 i++;
00676                                         }
00677                                 }
00678                         }
00679                 }
00680         }
00681         return GiNaC::lst(ret,matrix_to_lst2(A), basis);
00682 }
00683 
00684 
00685 GiNaC::lst polv(unsigned int no_fields, unsigned int order, unsigned int nsd, const string a)
00686 {
00687         GiNaC::lst ret1;                         // contains the polynom
00688         GiNaC::lst ret2;                         // contains the coefficients
00689         GiNaC::lst ret3;                         // constains the basis functions
00690         GiNaC::lst basis_tmp;
00691         for (unsigned int i=0; i< no_fields; i++)
00692         {
00693                 GiNaC::lst basis;
00694                 std::ostringstream s;
00695                 s <<a<<""<<i<<"_";
00696                 GiNaC::ex polspace = pol(order, nsd, s.str());
00697                 ret1.append(polspace.op(0));
00698                 ret2.append(polspace.op(1));
00699                 basis_tmp = GiNaC::ex_to<GiNaC::lst>(polspace.op(2));
00700                 for (GiNaC::lst::const_iterator basis_iterator = basis_tmp.begin();
00701                         basis_iterator != basis_tmp.end(); ++basis_iterator)
00702                 {
00703                         GiNaC::lst tmp_lst;
00704                         for (unsigned int d=1; d<=no_fields; d++) tmp_lst.append(0);
00705                         tmp_lst.let_op(i) = (*basis_iterator);
00706                         ret3.append(tmp_lst);
00707                 }
00708         }
00709         return GiNaC::lst(ret1,ret2,ret3);
00710 
00711         /* Old Code:
00712            GiNaC::lst ret;
00713            for (int i=1; i<= nsd; i++) {
00714            std::ostringstream s;
00715            s <<a<<"^"<<i<<"_";
00716            GiNaC::ex p = pol(order, nsd, s.str());
00717         ret.append(p);
00718         }
00719         return ret;
00720         */
00721 }
00722 
00723 
00724 GiNaC::ex polb(unsigned int order, unsigned int nsd, const string a)
00725 {
00726         using SyFi::x;
00727         using SyFi::y;
00728         using SyFi::z;
00729 
00730         GiNaC::ex ret;                           // ex to return
00731         int dof;                                         // degrees of freedom
00732         GiNaC::ex A;                             // ex holding the coefficients a_0 .. a_dof
00733         GiNaC::lst basis;
00734 
00735         if (nsd == 1)
00736         {
00737                 /* 1D:
00738                  * P^n = a_0 + a_1*x + .... + a_n*x^n
00739                  * dof : n+1
00740                  */
00741                 dof = order+1;
00742                 A = get_symbolic_matrix(1,dof, a);
00743                 int o=0;
00744                 for (GiNaC::const_iterator i = A.begin(); i != A.end(); ++i)
00745                 {
00746                         ret += (*i)*pow(x,o);
00747                         basis.append(pow(x,o));
00748                         o++;
00749                 }
00750         }
00751         else if ( nsd == 2)
00752         {
00753 
00754                 /* 2D: structure of coefficients (a_i)
00755                  * [ a_0      a_1 x     a_3 x^2     a_6 x^3
00756                  * [ a_2 y    a_4 xy    a_7 x^2y
00757                  * [ a_5 y^2  a_8 xy^2
00758                  * [ a_9 y^3
00759                  */
00760 
00761                 dof = (order+1)*(order+1);
00762                 A = get_symbolic_matrix(1, dof , a);
00763 
00764                 size_t i=0;
00765                 for (unsigned int o = 0; o <= order; o++)
00766                 {
00767                         for (unsigned int d = 0; d <= order; d++)
00768                         {
00769                                 ret += A.op(i)*pow(y,d)*pow(x,o);
00770                                 basis.append(pow(y,d)*pow(x,o));
00771                                 i++;
00772                         }
00773                 }
00774         }
00775         else if (nsd == 3)
00776         {
00777 
00778                 /* Similar structure as in 2D, but
00779                  * structured as a tetraheder, i.e.,
00780                  *   a_o + a_1 x + a_2 y + a_3 z
00781                  * + a_4 x^2 + a_5 xy +
00782                  */
00783                 dof = (order+1)*(order+1)*(order+1);
00784                 A = get_symbolic_matrix(1, dof , a);
00785 
00786                 size_t i=0;
00787                 for (unsigned int o = 0; o <= order; o++)
00788                 {
00789                         for (unsigned int d = 0; d <= order; d++)
00790                         {
00791                                 for (unsigned int f = 0; f <= order; f++)
00792                                 {
00793                                         ret += A.op(i)*pow(y,f)*pow(z,d)*pow(x,o);
00794                                         basis.append(pow(y,f)*pow(z,d)*pow(x,o));
00795                                         i++;
00796                                 }
00797                         }
00798                 }
00799         }
00800 
00801         return GiNaC::lst(ret,matrix_to_lst2(A), basis);
00802 }
00803 
00804 
00805 GiNaC::lst coeffs(GiNaC::lst pols)
00806 {
00807         GiNaC::lst cc;
00808         GiNaC::lst tmp;
00809         for (unsigned int i=0; i<= pols.nops()-1; i++)
00810         {
00811                 tmp = coeffs(pols.op(i));
00812                 cc = collapse(GiNaC::lst(cc, tmp));
00813         }
00814         return cc;
00815 }
00816 
00817 
00818 GiNaC::lst coeffs(GiNaC::ex pol)
00819 {
00820         using SyFi::x;
00821         using SyFi::y;
00822         using SyFi::z;
00823 
00824         GiNaC::lst cc;
00825         GiNaC::ex c, b;
00826         for (int i=pol.ldegree(x); i<=pol.degree(x); ++i)
00827         {
00828                 for (int j=pol.ldegree(y); j<=pol.degree(y); ++j)
00829                 {
00830                         for (int k=pol.ldegree(z); k<=pol.degree(z); ++k)
00831                         {
00832                                 c = pol.coeff(x,i).coeff(y, j).coeff(z,k);
00833                                 if ( c != 0 ) cc.append(c);
00834                         }
00835                 }
00836         }
00837         return cc;
00838 }
00839 
00840 
00841 GiNaC::exvector coeff(GiNaC::ex pol)
00842 {
00843         using SyFi::x;
00844         using SyFi::y;
00845         using SyFi::z;
00846 
00847         GiNaC::exvector cc;
00848         GiNaC::ex c, b;
00849         for (int i=pol.ldegree(x); i<=pol.degree(x); ++i)
00850         {
00851                 for (int j=pol.ldegree(y); j<=pol.degree(y); ++j)
00852                 {
00853                         for (int k=pol.ldegree(z); k<=pol.degree(z); ++k)
00854                         {
00855                                 c = pol.coeff(x,i).coeff(y, j).coeff(z,k);
00856                                 if ( c != 0 ) cc.insert(cc.begin(),c);
00857                         }
00858                 }
00859         }
00860         return cc;
00861 }
00862 
00863 
00864 GiNaC::exmap pol2basisandcoeff(GiNaC::ex e, GiNaC::ex s)
00865 {
00866         if (GiNaC::is_a<GiNaC::symbol>(s))
00867         {
00868                 GiNaC::symbol ss = GiNaC::ex_to<GiNaC::symbol>(s);
00869                 e = expand(e);
00870                 GiNaC::ex c;
00871                 GiNaC::ex b;
00872                 GiNaC::exmap map;
00873                 for (int i=e.ldegree(ss); i<=e.degree(ss); ++i)
00874                 {
00875                         c = e.coeff(ss,i);
00876                         b = pow(ss,i);
00877                         map[b] = c;
00878                 }
00879                 return map;
00880         }
00881         else
00882         {
00883                 throw(std::invalid_argument("The second argument must be a symbol."));
00884         }
00885 }
00886 
00887 
00888 GiNaC::exmap pol2basisandcoeff(GiNaC::ex e)
00889 {
00890         using SyFi::x;
00891         using SyFi::y;
00892         using SyFi::z;
00893 
00894         e = expand(e);
00895         GiNaC::ex c;
00896         GiNaC::ex b;
00897         GiNaC::exmap map;
00898         for (int i=e.ldegree(x); i<=e.degree(x); ++i)
00899         {
00900                 for (int j=e.ldegree(y); j<=e.degree(y); ++j)
00901                 {
00902                         for (int k=e.ldegree(z); k<=e.degree(z); ++k)
00903                         {
00904                                 c = e.coeff(x,i).coeff(y, j).coeff(z,k);
00905                                 b = pow(x,i)*pow(y,j)*pow(z,k);
00906                                 map[b] = c;
00907                         }
00908                 }
00909         }
00910         return map;
00911 }
00912 
00913 
00914 GiNaC::ex legendre1D(const GiNaC::symbol x, unsigned int n)
00915 {
00916         GiNaC::ex P;
00917         // Rodrigue's formula for Legendre polynomial of 1D
00918         // The interval [-1, 1]
00919         P=1/(pow(2,n)*GiNaC::factorial(n))*GiNaC::diff(GiNaC::pow((x*x-1),n),x,n);
00920         // -----------------
00921         // The interval [0,1]
00922         //  GiNaC::ex xx = 2*x - 1;
00923         // P=1/(pow(2,2*n)*GiNaC::factorial(n))*GiNaC::diff(GiNaC::pow((xx*xx-1),n),x,n);
00924         return P;
00925 }
00926 
00927 
00928 GiNaC::ex legendre(unsigned int order, unsigned int nsd, const string s)
00929 {
00930         using SyFi::x;
00931         using SyFi::y;
00932         using SyFi::z;
00933 
00934         // The Legendre polynomials to be used in FiniteElement
00935         GiNaC::ex leg;
00936         GiNaC::ex A;
00937         GiNaC::lst basis;
00938         int dof;
00939 
00940         GiNaC::ex b;
00941 
00942         // 1D
00943         if(nsd == 1)
00944         {
00945                 dof = order+1;
00946                 A = get_symbolic_matrix(1,dof,s);
00947                 int o=0;
00948                 for(GiNaC::const_iterator i = A.begin(); i!=A.end(); ++i)
00949                 {
00950                         b= legendre1D(x,o);
00951                         leg+= (*i)*b;
00952                         basis.append(b);
00953                         o++;
00954                 }
00955         }
00956         // 2D
00957         /*
00958         else if(nsd == 2){  // NB: Only for tensor products on TRIANGLES (not boxes)
00959                 / * 2D: structure of coefficients (a_i)
00960                  * [ a_0           a_1 P_1(x)           a_3 P_2(x)        a_6 P_3(x)
00961                  * [ a_2 P_1(y)    a_4 P_1(x)*P_1(y)    a_7 P_2(x)*P_1(y)
00962                  * [ a_5 P_2(y)    a_8 P_1(x)*P_2(y)
00963         * [ a_9 P_3(y)
00964         * /
00965         dof = (order+1)*(order+2)/2;
00966         A = get_symbolic_matrix(1,dof,s);
00967         size_t i=0;
00968         for (int o = 0; o <= order; o++) {
00969         for (int d = 0; d <= o; d++) {
00970         b = legendre1D(y,d)*legendre1D(x,o-d);
00971         leg += A.op(i)*b;
00972         basis.append(b);
00973         i++;
00974 
00975         }
00976         }
00977         }
00978         */
00979         else if(nsd == 2)                        // NB: Only for tensor products on rectangles
00980         {
00981                 dof = (order+1)*(order+1);
00982                 A = get_symbolic_matrix(1,dof,s);
00983                 size_t i=0;
00984                 for (unsigned int o = 0; o <= order; o++)
00985                 {
00986                         for (unsigned int d = 0; d <= order; d++)
00987                         {
00988                                 b = legendre1D(y,d)*legendre1D(x,o);
00989                                 leg += A.op(i)*b;
00990                                 basis.append(b);
00991                                 i++;
00992 
00993                         }
00994                 }
00995         }
00996 
00997         /* tetrahedron
00998         else if(nsd==3){
00999                 dof = 0;
01000                 for (int j=0; j<= order; j++) {
01001                         dof += (j+1)*(j+2)/2;
01002                 }
01003         A = get_symbolic_matrix(1, dof , s);
01004 
01005         size_t i=0;
01006         for (int o = 0; o <= order; o++) {
01007         for (int d = 0; d <= o; d++) {
01008         for (int f = 0; f <= o; f++) {
01009         if ( o-d-f >= 0) {
01010         b = legendre1D(y,f)*legendre1D(z,d)*legendre1D(x,o-d-f);
01011         leg += A.op(i)*b;
01012         basis.append(b);
01013         i++;
01014         }
01015         }
01016         }
01017         }
01018         }
01019         */
01020 
01021         else if(nsd==3)
01022         {
01023                 dof = (order+1)*(order+1)*(order+1);
01024                 A = get_symbolic_matrix(1, dof , s);
01025 
01026                 size_t i=0;
01027                 for (unsigned int o = 0; o <= order; o++)
01028                 {
01029                         for (unsigned int d = 0; d <= order; d++)
01030                         {
01031                                 for (unsigned int f = 0; f <= order; f++)
01032                                 {
01033                                         b = legendre1D(y,f)*legendre1D(z,d)*legendre1D(x,o);
01034                                         leg += A.op(i)*b;
01035                                         basis.append(b);
01036                                         i++;
01037                                 }
01038                         }
01039                 }
01040         }
01041         return GiNaC::lst(leg,matrix_to_lst2(A), basis);
01042 }
01043 
01044 
01045 GiNaC::lst legendrev(unsigned int no_fields, unsigned int order, unsigned int nsd, const string a)
01046 {
01047         GiNaC::lst ret1;                         // contains the polynom
01048         GiNaC::lst ret2;                         // contains the coefficients
01049         GiNaC::lst ret3;                         // constains the basis functions
01050         GiNaC::lst basis_tmp;
01051         for (unsigned int i=1; i<= no_fields; i++)
01052         {
01053                 GiNaC::lst basis;
01054                 std::ostringstream s;
01055                 s <<a<<""<<i<<"_";
01056                 GiNaC::ex polspace = legendre(order, nsd, s.str());
01057                 ret1.append(polspace.op(0));
01058                 ret2.append(polspace.op(1));
01059                 basis_tmp = GiNaC::ex_to<GiNaC::lst>(polspace.op(2));
01060                 for (GiNaC::lst::const_iterator basis_iterator = basis_tmp.begin();
01061                         basis_iterator != basis_tmp.end(); ++basis_iterator)
01062                 {
01063                         GiNaC::lst tmp_lst;
01064                         for (unsigned int d=1; d<=no_fields; d++) tmp_lst.append(0);
01065                         tmp_lst.let_op(i-1) = (*basis_iterator);
01066                         ret3.append(tmp_lst);
01067                 }
01068         }
01069         return GiNaC::lst(ret1,ret2,ret3);
01070 }
01071 
01072 
01073 bool compare(const ex & e, const string & s)
01074 {
01075         ostringstream ss;
01076         ss << e;
01077         return ss.str() == s;
01078 }
01079 
01080 
01081 void EQUAL_OR_DIE(const ex & e, const string & s)
01082 {
01083         if (!compare(e, s))
01084         {
01085                 ostringstream os;
01086                 os << "ERROR: expression e: " <<e<<" is not equal to "<<s<<endl;
01087                 throw runtime_error(os.str());
01088         }
01089 }
01090 
01091 
01092 // =========== expression inspection utilities
01093 
01094 class SymbolMapBuilderVisitor:
01095 public visitor,
01096         public symbol::visitor
01097 {
01098         public:
01099                 SymbolMapBuilderVisitor(GiNaC::exmap & em): symbolmap(em) {}
01100 
01101                 GiNaC::exmap & symbolmap;
01102 
01103                 void visit(const symbol & s)
01104                 {
01105                         GiNaC::exmap::iterator it = symbolmap.find(s);
01106                         if(it != symbolmap.end())
01107                         {
01108                                 it->second++;
01109                         }
01110                         else
01111                         {
01112                                 //GiNaC::exmap::value_type p(ex(s), ex(s));
01113                                 //symbolmap.insert(p);
01114                                 symbolmap[ex(s)] = ex(s);
01115                         }
01116                 }
01117 };
01118 
01119 class SymbolCounterVisitor:
01120 public visitor,
01121         public symbol::visitor,
01122         public basic::visitor
01123 {
01124         public:
01125                 exhashmap<int> symbolcount;
01126 
01127                 void visit(const basic & s)
01128                 {
01129                         std::cout << "visiting basic " << std::endl;
01130                 }
01131 
01132                 void visit(const symbol & s)
01133                 {
01134                         ex e = s;
01135                         std::cout << "visiting symbol " << e << std::endl;
01136                         exhashmap<int>::iterator it = symbolcount.find(s);
01137                         if(it != symbolcount.end())
01138                         {
01139                                 std::cout << "found symbol " << e << std::endl;
01140                                 it->second++;
01141                         }
01142                         else
01143                         {
01144                                 std::cout << "adding symbol " << e << std::endl;
01145                                 pair<ex,int> p;
01146                                 p.first = ex(s);
01147                                 p.second = 1;
01148                                 symbolcount.insert(p);
01149                         }
01150                 }
01151 };
01152 
01153 exhashmap<int> count_symbols(const ex & e)
01154 {
01155         SymbolCounterVisitor v;
01156         e.traverse(v);
01157         return v.symbolcount;
01158 }
01159 
01160 
01161 /*
01162 GiNaC::exvector get_symbols3(const GiNaC::ex & e)
01163 {
01164         // Implemented directly to avoid copying map:
01165         SymbolCounterVisitor v;
01166         e.traverse(v);
01167 exhashmap<int> & sc = v.symbolcount;
01168 
01169 int i = 0;
01170 GiNaC::exvector l(sc.size());
01171 for(exhashmap<int>::iterator it=sc.begin(); it!=sc.end(); it++)
01172 {
01173 //l.push_back(it->first);
01174 l[i] = it->first;
01175 i++;
01176 }
01177 return l;
01178 }
01179 
01180 std::list<GiNaC::ex> get_symbols2(const ex & e)
01181 {
01182 // Implemented directly to avoid copying map:
01183 SymbolCounterVisitor v;
01184 e.traverse(v);
01185 exhashmap<int> & sc = v.symbolcount;
01186 
01187 list<ex> l;
01188 for(exhashmap<int>::iterator it=sc.begin(); it!=sc.end(); it++)
01189 {
01190 l.push_back(it->first);
01191 }
01192 return l;
01193 }
01194 
01195 void get_symbols(const ex & e, GiNaC::exmap & em)
01196 {
01197 SymbolMapBuilderVisitor v(em);
01198 e.traverse(v);
01199 }
01200 */
01201 ex extract_symbols(const ex & e)
01202 {
01203         // Implemented directly to avoid copying map:
01204         SymbolCounterVisitor v;
01205         e.traverse(v);
01206         exhashmap<int> & sc = v.symbolcount;
01207 
01208         lst l;
01209         for(exhashmap<int>::iterator it=sc.begin(); it!=sc.end(); it++)
01210         {
01211                 l.append(it->first);
01212                 std::cout << (it->first) << std::endl;
01213         }
01214         ex ret = l;
01215         return ret;
01216 }
01217 
01218 
01219 // Collect all symbols of an expression
01220 void collect_symbols(const GiNaC::ex & e, exset & v)
01221 {
01222         if (GiNaC::is_a<GiNaC::symbol>(e))
01223         {
01224                 v.insert(e);
01225         }
01226         else
01227         {
01228                 for (size_t i=0; i<e.nops(); i++)
01229                 {
01230                         collect_symbols(e.op(i), v);
01231                 }
01232         }
01233 }
01234 
01235 
01236 GiNaC::exvector collect_symbols(const GiNaC::ex & e)
01237 {
01238         exset s;
01239         collect_symbols(e, s);
01240         GiNaC::exvector v(s.size());
01241         for(exset::iterator i=s.begin(); i!= s.end(); i++)
01242         {
01243                 v.push_back(*i);
01244         }
01245         return v;
01246 }
01247 
01248 
01249 bool compare_archives(const string & first, const string & second, std::ostream & os)
01250 {
01251         bool ret = true;
01252 
01253         // read both archives
01254         archive a1, a2;
01255         ifstream if1(first.c_str()), if2(second.c_str());
01256         if1 >> a1;
01257         if2 >> a2;
01258 
01259         // compare size
01260         int n = a1.num_expressions();
01261         int n2 = a2.num_expressions();
01262         if(n != n2)
01263         {
01264                 os << "Archives " << first << " and " << second
01265                         << " has a different number of expressions, " << n << " and " << n2 << "." << endl;
01266                 os << "Comparing common expressions." << endl;
01267                 ret = false;
01268         }
01269 
01270         // iterate over all expressions in first archive
01271         ex e1,e2;
01272         for(int i=0; i<n; i++)
01273         {
01274                 lst syms;
01275                 string exname;
01276 
01277                 e1 = a1.unarchive_ex(syms, exname, i);
01278 
01279                 syms = ex_to<lst>(extract_symbols(e1));
01280                 //        os << "Comparing " << exname << " with symbols " << syms << endl;
01281 
01282                 // is this in the second archive?
01283                 try
01284                 {
01285                         e2 = a2.unarchive_ex(syms, exname.c_str());
01286 
01287                         // got it, now compare
01288                         bool isequal = is_zero(e1-e2);
01289                         if(!isequal)
01290                         {
01291                                 if(ret)
01292                                 {
01293                                         os << "Archives " << first << " and " << second
01294                                                 << " are not equal, details follow:" << endl;
01295                                 }
01296                                 os << "Expression with name " << exname << " is not equal:" << endl;
01297                                 os << "First:  " << endl << e1 << endl;
01298                                 os << "Second: " << endl << e2 << endl;
01299                                 ret = false;
01300                         }
01301                 }
01302                 catch(...)
01303                 {
01304                         os << "Expression " << exname << " is missing from " << second << "." << endl;
01305                         ret = false;
01306                 }
01307         }
01308 
01309         return ret;
01310 }
01311 
01312 
01313 class ExStatsVisitor:
01314 public visitor,
01315         public power::visitor,
01316         public mul::visitor,
01317         public add::visitor,
01318         public function::visitor
01319 {
01320         public:
01321                 ExStats es;
01322 
01323                 void visit(const power & e)
01324                 {
01325                         //cout << "pow" << endl;
01326                         ex b = e.op(1);
01327                         ex c = b.evalf();
01328                         if( is_a<numeric>(c) )
01329                         {
01330                                 numeric nu = ex_to<numeric>(c);
01331                                 int n = (int) to_double(nu);
01332                                 if( is_zero(nu-n) )
01333                                 {
01334                                         // is an integer power, interpret as a sequence of multiplications
01335                                         es.muls  += n-1;
01336                                         es.flops += n-1;
01337                                 }
01338                         }
01339 
01340                         es.pows++;
01341                 }
01342 
01343                 void visit(const mul & e)
01344                 {
01345                         //cout << "mul" << endl;
01346                         es.muls  += e.nops()-1;
01347                         es.flops += e.nops()-1;
01348                 }
01349 
01350                 void visit(const add & e)
01351                 {
01352                         //cout << "add" << endl;
01353                         es.adds  += e.nops()-1;
01354                         es.flops += e.nops()-1;
01355                 }
01356 
01357                 void visit(const function & e)
01358                 {
01359                         //cout << "func" << endl;
01360                         es.functions++;
01361                 }
01362 };
01363 
01364 ExStats count_ops(const ex & e)
01365 {
01366         //cout << "count_ops " << e << endl;
01367         //cout << "is an add: " << GiNaC::is<GiNaC::add>(e) << endl;
01368         //cout << "is a  mul: " << GiNaC::is<GiNaC::mul>(e) << endl;
01369         ExStatsVisitor v;
01370         e.traverse(v);
01371         return v.es;
01372 }
01373 
01374 
01375 // =========== expression manipulation utilities
01376 
01377 // Returns in sel a list with symbol/expression pairs for all
01378 // power replacement variables introduced. Works best on
01379 // expressions in expanded form.
01380 ex replace_powers(const ex & ein, const list<symbol> & symbols, list<symexpair> & sel, const string & tmpsymbolprefix)
01381 {
01382         ex e = ein;
01383         // build power expressions
01384         list<symbol>::const_iterator it = symbols.begin();
01385         for(; it != symbols.end(); it++)
01386         {
01387                 int deg      = e.degree(*it);
01388                 if(deg > 0)
01389                 {
01390                         symbol sym   = ex_to<symbol>(*it);
01391                         string sname = tmpsymbolprefix + sym.get_name();
01392 
01393                         // make list of new symbols
01394                         vector<symbol> symbols(deg);
01395                         symbols[0] = sym;
01396                         for(int i=1; i<deg; i++)
01397                         {
01398                                 symbols[i] = get_symbol( sname + int2string(i+1) );
01399                                 sel.push_back(make_pair(symbols[i], symbols[i-1]*sym));
01400                         }
01401 
01402                         // with highest order first, subs in e
01403                         ex prod = sym;
01404                         for(int i=deg-1; i>=1; i--)
01405                         {
01406                                 e = e.subs(power(sym,i+1) == symbols[i], subs_options::algebraic);
01407                         }
01408                 }
01409         }
01410         return e;
01411 }
01412 
01413 
01414 };                                                               // namespace SyFi
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines