Hi, On Mon, Feb 23, 2009 at 12:37:39PM -0500, Jeremy Jay wrote:
A short example is the following code:
----
#include <iostream> #include <ginac/ginac.h>
using namespace std; using namespace GiNaC;
//////////////////////////////////////////////////////////////
// delta( c_i, c_j ) = 1 if c_i==c_j otherwise 0 DECLARE_FUNCTION_2P(delta); REGISTER_FUNCTION(delta, dummy());
//////////////////////////////////////////////////////////////
int main(int argc, char **argv) { symbol q("q"), c0("c0"), c1("c1"), c2("c2");
ex poly( (1-delta(c0,c1))*(1-delta(c0,c2)) ); poly=poly.expand();
cout << "START : " << poly << endl;
poly = poly.subs(1 == q);
cout << "sub 1 => q : " << poly << endl;
return 0; }
----
which gives the output:
----
START : 1-delta(c0,c2)-delta(c0,c1)+delta(c0,c2)*delta(c0,c1) sub 1 => q : 1-delta(c0,c2)-delta(c0,c1)+delta(c0,c2)*delta(c0,c1)
----
which is not what I would expect. I've also tried using the 'numeric' type instead of the literal 1's, but get the same effect. Am I missing something?
Although there's nothing wrong about your code, such substitution won't work. Rewrite the `poly' expression as symbol A("A"); ex poly((A-delta(c0,c1))*(A-delta(c0,c2)) ); and than substitute A with whatever you need (1, q, etc). In general, subs() and friends (match(), has()) operate on internal representation of expressions (as opposed to their mathematical properties). GiNaC stores sums (and products) as (overall_coeff . ((coeff_0 . rest_0) (coeff_1 . rest_1) ... )) where overall_coeff, coeff_0, coeff_1, ... are numbers, and rest_0, rest_1 are arbitrary expressions. subs() operates only on `rest' parts (there are some tricks for products, but that's a different story): ginac/expairseq.cpp 1710 // Substitute only in the "rest" part of the pairs 1711 epvector::const_iterator cit = seq.begin(), last = seq.end(); 1712 while (cit != last) { 1713 1714 const ex &subsed_ex = cit->rest.subs(m, options); 1715 if (!are_ex_trivially_equal(cit->rest, subsed_ex)) { Making substitution in "coeff" part(s) would give a (very) surprising results, i.e. (1 + x*y).subs(1 == q) => q + x^q*y^q (1/2*x^2/z^2).subs(2 == q) => 1/2*x^q/z^2 Best regards, Alexei