small GiNaC query
keith.briggs at bt.com
keith.briggs at bt.com
Tue Jan 27 17:41:05 CET 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
More information about the GiNaC-list
mailing list