[GiNaC-list] Expression order changing

[Sheplyakov Alexei]

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.

-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 827 bytes
Desc: Digital signature
Url : http://www.cebix.net/pipermail/ginac-list/attachments/20070605/3612bdf8/attachment.pgp