[GiNaC-list] Substitution troubles

[Alexei Sheplyakov]

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

-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 835 bytes
Desc: Digital signature
URL: <http://www.cebix.net/pipermail/ginac-list/attachments/20090224/95d7b242/attachment.sig>