Hello! On Mon, Jun 04, 2007 at 10:33:21AM +0200, Martin Sandve Alnæs wrote:
Hi, when generating code from ginac expressions (with csrc), the order of the printed expressions will vary depending on previous operations in the program, in particular the previous construction of symbols. I understand this is because ginac uses an internal symbol id counter of some kind.
Is it possible to "reset" this counter, or somehow make sure a set of newly created symbols keep the same relative order independent of previous operations in the program?
You could create custom printing function (or print_context) which does not care of symbol ordering. Here is what I use: #include <iostream> #include <ginac/ginac.h> using namespace GiNaC; #define unlikely(cond) __builtin_expect((cond), 0) /** * @brief transform multivariate polynomial into an exmap * @param e exporession to operate on * @param l list of variables * @param do_expand -- expand expression before transforming it. * Set this to false if you know in advance that e is expanded w.r.t. * all variables in l. Useful for processing polynomials of the form * x_1^j_1...x_k^j_k(a_i+b_i+c_i...)^n_i. * BEWARE: if the initial expression is not expanded in variables * listed in l, the result is complete nonsense. */ const exmap collect_vargs(const ex& e, const lst& l, const bool do_expand) { static const ex _ex1(1); static const ex _ex0(0); ex x; exmap cmap = exmap(); if (unlikely(!do_expand)) x = e; else x = e.expand(); if (unlikely(! is_a<add>(x))) { ex key = _ex1; ex pre_coeff = x; for (lst::const_iterator li=l.begin(); li!=l.end(); ++li) { int cexp = pre_coeff.degree(*li); pre_coeff = pre_coeff.coeff(*li, cexp); key *= pow(*li, cexp); } cmap.insert(exmap::value_type(key, pre_coeff)); } else { for (const_iterator xi=x.begin(); xi!=x.end(); ++xi) { ex key = _ex1; ex pre_coeff = *xi; for (lst::const_iterator li=l.begin(); li!=l.end(); ++li) { int cexp = pre_coeff.degree(*li); pre_coeff = pre_coeff.coeff(*li, cexp); key *= pow(*li, cexp); } exmap::iterator ci = cmap.find(key); if (ci != cmap.end()) ci->second += pre_coeff; else cmap.insert(exmap::value_type(key, pre_coeff)); } } return cmap; } #define MARK \ (std::cout << std::endl << __FUNCTION__ << ":" << __LINE__ << std::endl) void example0() { MARK; symbol x("x"), y("y"), a("a"), b("b"), c("c"); ex e = pow(x, 2)*a + x*y*pow(b+c, 10); exmap em = collect_vargs(e, lst(x, y), false); // N.B. e is already expanded in x, y for (exmap::const_iterator ii=em.begin(); ii!=em.end(); ++ ii) std::cout << ii->first << " => " << ii->second << std::endl; // prints // x^2 => a // x*y => (a+b)^10 } /** * @brief Print monomial in a human-readable form * @param l user-specified ordering, a list of expressions * @param e expression to print, should be a monomial (w.r.t. l) * */ inline void print_monomial_ordered(std::ostream& s, const ex& e, const lst& l) { static const ex _ex1(1); if (! (is_a<mul>(e) || is_a<power>(e))) { s << e; return; } // split the monomial into a vector... exvector ev; ex cmp = _ex1; for (lst::const_iterator li=l.begin(); li!=l.end(); ++li) { int cexp = e.degree(*li); ex tmp = pow(*li, cexp); if (tmp.is_equal(_ex1)) continue; ev.push_back(tmp); cmp *= tmp; // N.B.: vector will be automagically sorted. if (cmp.is_equal(e)) break; } // ...and print out the vector size_t nops = ev.size(); size_t count = 0; for (exvector::const_iterator vi=ev.begin(); vi!=ev.end(); ++vi) { bool need_parens = false; if (is_a<add>(*vi)) need_parens = true; if (need_parens) s << '('; s << *vi; if (need_parens) s << ')'; ++count; if (count != nops) s << '*'; } } void example1() { MARK; symbol y("y"),z("z"); // shuffle symbols for (unsigned i=0; i < 100; i++) { symbol* xxx = new symbol("xxx"); // do something stupid xxx->set_name("foo"); delete xxx; } symbol x("x"); ex e = pow(x, 2)*pow(z,3)*y; print_monomial_ordered(std::cout, e, lst(x, y, z)); // prints // x^2*y*z^3 // not y*z^3*x^2, z^3*x^2*y, etc. std::cout << std::endl; } /** * @brief print collected @var{add} object one term per line */ inline void print_by_one(std::ostream& s, const exmap& e, const lst& l) { for (exmap::const_iterator ei=e.begin(); ei != e.end(); ++ei) { s << "+("; print_monomial_ordered(s, ei->first, l); s << ")*(" << ei->second << ")" << std::endl; } } void example2() { MARK; symbol x("x"), y("y"), a("a"), b("b"), c("c"); ex e = pow(x, 2)*a + x*y*pow(b+c, 10); lst l(x, y); exmap em = collect_vargs(e, l, false); // N.B. e is already expanded in x, y, hense `false' argument // to collect_vargs print_by_one(std::cout, em, l); // prints // // +(x^2)*(a) // +(x*y)*((b+c)^10) } /** * @brief print multivariate polynomial one term per line * * The order of terms is still random, but expressions are *look* much * better. Pipe through sort(1) to get even better result. * * @param s destination * @param e expression to print * @param l list of variables * @param need_expand (default: true) -- expand expression before processing. * This is the safe choice, but sometimes it spoils nice properties of * the input expression (e.g., partial factorization). Set to false only if * you know in advance that the expression (e) is already expanded (w.r.t. l). * You have been warned. * * @todo: find out how to determine this in a fast and reliable way. */ void print_add_by_one(std::ostream& s, const ex& e, const lst& l, const bool need_expand=true) { exmap em = collect_vargs(e, l, need_expand); // N.B.: .normal() makes coefficients much "simpler" (I know in // advance that some denominators in my expressions *must* cancel). // But in general case this might be just a waste of time. So don't // hesitate to comment out these 3 lines: for (exmap::iterator i=em.begin(); i!=em.end(); ++i) i->second = collect_common_factors(i->second.normal()); print_by_one(s, em, l); } void example3() { MARK; symbol x("x"), y("y"), a("a"), b("b"), c("c"); ex e = pow(x, 2)*a + x*y*pow(b+c, 10); // N.B. It is simple to see that e is already expanded in x, y, // hense `false' argument to collect_vargs is safe. print_add_by_one(std::cout, e, lst(x,y), false); // prints // +(x^2)*(a) // +(x*y)*((b+c)^10) } int main(int argc, char** argv) { example0(); example1(); example2(); example3(); return 0; } Best regards, Alexei -- All science is either physics or stamp collecting.