Hello all, I'm new to GiNaC, but so far I've been very impressed by its power. But, I seem to be having a problem with the substitutions, can anyone tell me what I'm doing wrong? 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? BTW, for anyone interesed in the point of this, I'm trying to implement the algorithm in this paper: http://www.iop.org/EJ/abstract/1367-2630/11/2/023001/ Thanks in advance! Jeremy
This is for Jeremy and not the mailing list. Can this method be used for the Jones polynomial (and other link polynomials)? Bruce Westbury -----Original Message----- From: ginac-list-bounces@ginac.de on behalf of Jeremy Jay Sent: Mon 23-Feb-09 5:37 PM To: ginac-list@ginac.de Subject: [GiNaC-list] Substitution troubles Hello all, I'm new to GiNaC, but so far I've been very impressed by its power. But, I seem to be having a problem with the substitutions, can anyone tell me what I'm doing wrong? 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? BTW, for anyone interesed in the point of this, I'm trying to implement the algorithm in this paper: http://www.iop.org/EJ/abstract/1367-2630/11/2/023001/ Thanks in advance! Jeremy _______________________________________________ GiNaC-list mailing list GiNaC-list@ginac.de https://www.cebix.net/mailman/listinfo/ginac-list
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
On Tue 24 Feb 2009 - 09:53AM, Alexei Sheplyakov wrote:
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
Hmm I guess it makes sense in that regard, but it is certainly non-intuitive. Printing the poly shows exactly what I want to replace and not the interal representation. Might I suggest someone update the tutorial with a statement like "Numeric literals are strongly discouraged in the left side of a substitution, and will not produce the desired effect." Using symbols and doing these extra replacements makes what was a pretty simple algorithm much messier... Thanks Jeremy
On Tue, 24 Feb 2009 08:50:46 -0500, Jeremy Jay <dinkumator@gmail.com> said: JJ> Hmm I guess it makes sense in that regard, but it is certainly JJ> non-intuitive.
Please consider two expressions: ((x+2)^2).subs(4==A) and (x^2+4x+4).subs(4==A). What are "intuitive" answers? Thus even if you are not referring to "internal representation" this is the only way to make subs() mathematically correct. I think Alexei made his point absolutely correct. Best wishes, Vladimir -- Vladimir V. Kisil email: kisilv@maths.leeds.ac.uk -- www: http://maths.leeds.ac.uk/~kisilv/
Hi, On Tue, Feb 24, 2009 at 08:50:46AM -0500, Jeremy Jay wrote:
Hmm I guess it makes sense in that regard, but it is certainly non-intuitive.
In every CAS I know of substitution (and pattern matching) operates on the internal representation of expressions, not on their mathematical meaning. For example, in M.........a (version 5.1): In[1]:= x^4 + 1/(16*y^4); In[2]:= InputForm[% /. a_^2 :> a] Out[3]//InputForm= x^4 + 1/(16*y^4) In[4]:= InputForm[x^4 + 1/(16*y^4) /. 1 -> q] Out[5]//InputForm= x^4 + 1/(16*y^4) Not quite intuitive, is it?
Printing the poly shows exactly what I want to replace and not the interal representation.
cout << tree << poly << endl; prints the internal representation.
Might I suggest someone update the tutorial with a statement like "Numeric literals are strongly discouraged in the left side of a substitution, and will not produce the desired effect."
That's kind of vague, because it's not quite clear what `desired effect' exactly is. Best regards, Alexei
It's not "intuitive" because I can see a plain single "1" when i print the poly, and I can call subs(1==q) on the poly with no errors, but when I print the poly again it is unchanged. If I were to print your first example as written, I would not see a "4" in the output and I would not expect anything to be substituted. If I wanted to make sure I got the 4, I would have to call expand or whatever appropriately, before the substitution. In the second example I would expect it to come out to (x^2+A*x+A). I see 4s and I'm substituting As in their place. If I don't want the 4x changed I'll have to do some extra work myself, but I don't expect GiNaC to differentiate between the 4s I see. I think any of these cases are very "intuitive" for any programmer with any experience with any pattern matching. A substitution is just that, a substitution, and it is not "well-defined" mathematically and I don't really see any reason for it to be as rigorous as you imply. If someone wants to do silly things like replace 4s with As, let him! Alexei made his point, yes, but I don't see how it is "correct." He said that the semantics of a simple substitution rely upon the internal representation of the polynomial -- although there is no documentation on this point anywhere on the website. This make it non-intuitive. If you have a reference for a rigorous mathematical definition of substitution on numbers feel free to send it to me to correct me. Jeremy On Tue 24 Feb 2009 - 01:59PM, Vladimir V. Kisil wrote:
On Tue, 24 Feb 2009 08:50:46 -0500, Jeremy Jay <dinkumator@gmail.com> said: JJ> Hmm I guess it makes sense in that regard, but it is certainly JJ> non-intuitive.
Please consider two expressions:
((x+2)^2).subs(4==A) and (x^2+4x+4).subs(4==A).
What are "intuitive" answers? Thus even if you are not referring to "internal representation" this is the only way to make subs() mathematically correct.
I think Alexei made his point absolutely correct.
Best wishes, Vladimir -- Vladimir V. Kisil email: kisilv@maths.leeds.ac.uk -- www: http://maths.leeds.ac.uk/~kisilv/
Hi! Jeremy Jay wrote:
It's not "intuitive" because I can see a plain single "1" when i print the poly, and I can call subs(1==q) on the poly with no errors, but when I print the poly again it is unchanged.
If I were to print your first example as written, I would not see a "4" in the output and I would not expect anything to be substituted. If I wanted to make sure I got the 4, I would have to call expand or whatever appropriately, before the substitution.
In the second example I would expect it to come out to (x^2+A*x+A). I see 4s and I'm substituting As in their place. If I don't want the 4x changed I'll have to do some extra work myself, but I don't expect GiNaC to differentiate between the 4s I see.
I think any of these cases are very "intuitive" for any programmer with any experience with any pattern matching. A substitution is just that, a substitution, and it is not "well-defined" mathematically and I don't really see any reason for it to be as rigorous as you imply. If someone wants to do silly things like replace 4s with As, let him!
But GiNaC is not sed! If you want to carry over a programmer's naive notion of substitution, then you could always print the expression and do the replacement in a string buffer. I recommend using boost::regex.
Alexei made his point, yes, but I don't see how it is "correct." He said that the semantics of a simple substitution rely upon the internal representation of the polynomial -- although there is no documentation on this point anywhere on the website. This make it non-intuitive.
Alexei was referring to GiNaC's actual internal representation. It could be a different representation. I'm seeing 1s all over every mathematical expression. 1+x*y could as well be written as 1^(1^1)+1^1*x^(1^1)*y^(1^1) and, upon substitution, result in q^(q^q)+q^q*x^(q^q)*y^(q^q). Or as (-1)^2+x*y and be invariant under subs(1==q). Is this really more insane than GiNaC's actual representation?
If you have a reference for a rigorous mathematical definition of substitution on numbers feel free to send it to me to correct me.
That's an interesting twist. IIRC you were the one suggesting to do this substitution on numbers. Now, we all seem to agree that it's not well-defined. Are you proposing, then, that .subs(x==y) isn't well-defined either? But every mathematician is familiar with this notion of substituting free variables and considers it sane. Is it not possible to rewrite your algorithm substituting variables instead of numbers? Cheers -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>
OK Nevermind all. I'll use a better toolkit. Jeez On Tue 24 Feb 2009 - 11:25PM, Richard B. Kreckel wrote:
Hi!
Jeremy Jay wrote:
It's not "intuitive" because I can see a plain single "1" when i print the poly, and I can call subs(1==q) on the poly with no errors, but when I print the poly again it is unchanged. If I were to print your first example as written, I would not see a "4" in the output and I would not expect anything to be substituted. If I wanted to make sure I got the 4, I would have to call expand or whatever appropriately, before the substitution. In the second example I would expect it to come out to (x^2+A*x+A). I see 4s and I'm substituting As in their place. If I don't want the 4x changed I'll have to do some extra work myself, but I don't expect GiNaC to differentiate between the 4s I see. I think any of these cases are very "intuitive" for any programmer with any experience with any pattern matching. A substitution is just that, a substitution, and it is not "well-defined" mathematically and I don't really see any reason for it to be as rigorous as you imply. If someone wants to do silly things like replace 4s with As, let him!
But GiNaC is not sed! If you want to carry over a programmer's naive notion of substitution, then you could always print the expression and do the replacement in a string buffer. I recommend using boost::regex.
Alexei made his point, yes, but I don't see how it is "correct." He said that the semantics of a simple substitution rely upon the internal representation of the polynomial -- although there is no documentation on this point anywhere on the website. This make it non-intuitive.
Alexei was referring to GiNaC's actual internal representation. It could be a different representation. I'm seeing 1s all over every mathematical expression. 1+x*y could as well be written as 1^(1^1)+1^1*x^(1^1)*y^(1^1) and, upon substitution, result in q^(q^q)+q^q*x^(q^q)*y^(q^q). Or as (-1)^2+x*y and be invariant under subs(1==q). Is this really more insane than GiNaC's actual representation?
If you have a reference for a rigorous mathematical definition of substitution on numbers feel free to send it to me to correct me.
That's an interesting twist. IIRC you were the one suggesting to do this substitution on numbers. Now, we all seem to agree that it's not well-defined. Are you proposing, then, that .subs(x==y) isn't well-defined either? But every mathematician is familiar with this notion of substituting free variables and considers it sane. Is it not possible to rewrite your algorithm substituting variables instead of numbers?
Cheers -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/> _______________________________________________ GiNaC-list mailing list GiNaC-list@ginac.de https://www.cebix.net/mailman/listinfo/ginac-list
Hello, On Tue, Feb 24, 2009 at 10:07:12AM -0500, Jeremy Jay wrote:
It's not "intuitive" because I can see a plain single "1" when i print the poly, and I can call subs(1==q) on the poly with no errors, but when I print the poly again it is unchanged.
It's never going to be "what you see is what you get" (see the example in M.........a I've given in a previous mail).
If I were to print your first example as written, I would not see a "4" in the output and I would not expect anything to be substituted. If I wanted to make sure I got the 4, I would have to call expand or whatever appropriately, before the substitution.
In the second example I would expect it to come out to (x^2+A*x+A). I see 4s and I'm substituting As in their place. If I don't want the 4x changed I'll have to do some extra work myself, but I don't expect GiNaC to differentiate between the 4s I see.
I think any of these cases are very "intuitive" for any programmer with any experience with any pattern matching.
Even for strings there are a number of different matching rules (cf. `greedy matching' versus `non-greedy matching'). Obviously, there are much more different ways to define matching rules for mathematical expressions. And for every such definition there are people who like it (for a good reason) and those who hate it (for a good reason, too!).
A substitution is just that, a substitution, and it is not "well-defined" mathematically and I don't really see any reason for it to be as rigorous as you imply.
Mathematical definitions help to avoid confusion, because typically such definitions avoid highly subjective notions (such as `intuitive'). That said, I don't know any formal definition of pattern matching (for mathematical expressions). All mathematical software I know of operates on internal representation. At best, that representation and matching rules are documented, but that's it.
If someone wants to do silly things like replace 4s with As, let him!
No problem. The code below does something like that. However, it replaces _all_ '4s with As'. #include <ginac/ginac.h> #include <iostream> using namespace std; using namespace GiNaC; struct subs_silly : public map_function { ex patt; ex repl; ex operator()(const ex& e) { if (is_a<add>(e) || is_a<mul>(e)) return e.map(*this); if (is_a<power>(e)) { return power(e.op(0).map(*this), e.op(1).map(*this)); } if (e.is_equal(patt)) return repl; else return e; } }; int main(int argc, char** argv) { symbol x("x"), y("y"), q("q"); ex e = y*pow(x, 2) + 1; subs_silly se; se.patt = ex(1); se.repl = q; cout << e << " => "; e = se(e); cout << e << endl; // prints: 1+x^2*y => q+x^2*y return 0; }
Alexei made his point, yes, but I don't see how it is "correct." He said that the semantics of a simple substitution rely upon the internal representation of the polynomial -- although there is no documentation on this point anywhere on the website.
That's not quite true. The internal documentation is somewhat documented in the tutorial (Appendix A.2, titled `Internal representation of products and sums'). Also, GiNaC provides the `tree' printing context, which outputs the internal representation of expression: #include <ginac/ginac.h> #include <iostream> using namespace std; using namespace GiNaC; int main(int argc, char** argv) { symbol x("x"), y("y"); ex e = y*pow(x, 2) + 1; cout << "Internal representation of the expression \"" << e << "\" is: " << endl; cout << tree << e << endl << dflt; return 0; } Best regards, Alexei
participants (5)
-
Alexei Sheplyakov
-
Jeremy Jay
-
Richard B. Kreckel
-
Vladimir V. Kisil
-
Westbury, Bruce