00001
00002
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
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
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
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
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
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
00283 cols = 1;
00284 if (itr->nops() > cols)
00285 cols = itr->nops();
00286 }
00287
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)) ) {
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
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
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
00452 if (GiNaC::is_a<GiNaC::numeric>(e))
00453 {
00454 }
00455 else if (GiNaC::is_a<GiNaC::add>(e) )
00456 {
00457
00458
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
00464
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
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
00484
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
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
00505 check_visitor(e2,exlist);
00506 }
00507 }
00508
00509 else
00510 {
00511
00512
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;
00574 GiNaC::lst ret2;
00575 GiNaC::lst ret3;
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;
00606 int dof;
00607 GiNaC::ex A;
00608 GiNaC::lst basis;
00609
00610 if (nsd == 1)
00611 {
00612
00613
00614
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
00630
00631
00632
00633
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
00653
00654
00655
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;
00688 GiNaC::lst ret2;
00689 GiNaC::lst ret3;
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
00712
00713
00714
00715
00716
00717
00718
00719
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;
00731 int dof;
00732 GiNaC::ex A;
00733 GiNaC::lst basis;
00734
00735 if (nsd == 1)
00736 {
00737
00738
00739
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
00755
00756
00757
00758
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
00779
00780
00781
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
00918
00919 P=1/(pow(2,n)*GiNaC::factorial(n))*GiNaC::diff(GiNaC::pow((x*x-1),n),x,n);
00920
00921
00922
00923
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
00935 GiNaC::ex leg;
00936 GiNaC::ex A;
00937 GiNaC::lst basis;
00938 int dof;
00939
00940 GiNaC::ex b;
00941
00942
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
00957
00958
00959
00960
00961
00962
00963
00964
00965
00966
00967
00968
00969
00970
00971
00972
00973
00974
00975
00976
00977
00978
00979 else if(nsd == 2)
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
00998
00999
01000
01001
01002
01003
01004
01005
01006
01007
01008
01009
01010
01011
01012
01013
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;
01048 GiNaC::lst ret2;
01049 GiNaC::lst ret3;
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
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
01113
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
01163
01164
01165
01166
01167
01168
01169
01170
01171
01172
01173
01174
01175
01176
01177
01178
01179
01180
01181
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191
01192
01193
01194
01195
01196
01197
01198
01199
01200
01201 ex extract_symbols(const ex & e)
01202 {
01203
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
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
01254 archive a1, a2;
01255 ifstream if1(first.c_str()), if2(second.c_str());
01256 if1 >> a1;
01257 if2 >> a2;
01258
01259
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
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
01281
01282
01283 try
01284 {
01285 e2 = a2.unarchive_ex(syms, exname.c_str());
01286
01287
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
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
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
01346 es.muls += e.nops()-1;
01347 es.flops += e.nops()-1;
01348 }
01349
01350 void visit(const add & e)
01351 {
01352
01353 es.adds += e.nops()-1;
01354 es.flops += e.nops()-1;
01355 }
01356
01357 void visit(const function & e)
01358 {
01359
01360 es.functions++;
01361 }
01362 };
01363
01364 ExStats count_ops(const ex & e)
01365 {
01366
01367
01368
01369 ExStatsVisitor v;
01370 e.traverse(v);
01371 return v.es;
01372 }
01373
01374
01375
01376
01377
01378
01379
01380 ex replace_powers(const ex & ein, const list<symbol> & symbols, list<symexpair> & sel, const string & tmpsymbolprefix)
01381 {
01382 ex e = ein;
01383
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
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
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 };