small GiNaC query

27 Jan 2004


      How do I get the program below to compute
1/6*s.1^3+1/2*s.1*s.2+1/3*s.3
?   I don't care about ordering of summands, but I need the s.1*s.2 and s.2*s.1 terms to be merged
and s.1*s.1*s.1 to be converted to s.1^3.
Thanks,
Keith

// KMB 2004 Jan 26
// g++ -Wall tryginac1.cc -lginac && a.out

#include <iostream>
#include <stdexcept>
#include <ginac/ginac.h>

using namespace std;
using namespace GiNaC;

ex Z(int n) {
  if (n==0) return 1; 
  ex r=0;
  symbol s("s"),i_sym("i");
  idx i(i_sym,n);
  for (int j=1; j<=n; j++) {
    ex si=indexed(s,i); 
    r+=si.subs(i==j)*Z(n-j); 
  }
  return r.expand().collect(s)/n;          // -> 1/6*s.2*s.1+1/3*s.1*s.2+1/3*s.3+1/6*s.1*s.1*s.1
  return simplify_indexed(r/n);            // -> 1/6*s.2*s.1+1/3*s.1*s.2+1/3*s.3+1/6*s.1*s.1*s.1
  return collect_common_factors(r/n); // -> 1/6*s.1*(s.2+s.1*s.1)+1/3*s.3+1/3*s.1*s.2
  return r.expand()/n;                         // -> 1/6*s.1*s.1*s.1+1/6*s.1*s.2+1/3*s.3+1/3*s.1*s.2
}

int main(void) {
  try {
    cout<<Z(3)<<endl;
  } catch (exception &p) {
    cerr << p.what() << endl;
    return 1;
  }
  return 0;
}


	Dr. Keith M. Briggs
	Senior Mathematician, Complexity Research, BT Exact
	http://more.btexact.com/people/briggsk2/ (internet)
	http://research.btexact.com/teralab/keithbriggs.html (BT intranet)
	phone: +44(0)1473  work: 641 911 home: 610 517  fax: 642 161
	mail: Keith Briggs, Polaris 134, Adastral Park, Martlesham, Suffolk IP5 3RE, UK

keith.briggs＠bt.com

varg＠thsun1.jinr.ru

tags

participants (2)