[PATCH] clifford_unit: fix possible bugs due to wrong operator[!=]= usage.
The condition like this one unsigned n /* = something */; if (n == ex_to<idx>(mu).get_dim) { // do something } is *always* true, because == creates an object of type `relational'. Likewise, this condition matrix M /* = something */; if (M(i, j) != M(j, i)) { // more code } is *always* true too. Frankly, I've got annoyed by that evil `relational' class. It doesn't do anything useful and forces everyone to write the code in a very unnatural, counter-intuitive way. I plan to remove it in the next version of GiNaC (1.5). --- ginac/clifford.cpp | 17 ++++++++++++++--- 1 files changed, 14 insertions(+), 3 deletions(-) diff --git a/ginac/clifford.cpp b/ginac/clifford.cpp index 100954d..4d73d50 100644 --- a/ginac/clifford.cpp +++ b/ginac/clifford.cpp @@ -716,6 +716,17 @@ ex dirac_ONE(unsigned char rl) return clifford(ONE, rl); } +static unsigned get_dim_uint(const ex& e) +{ + if (!is_a<idx>(e)) + throw std::invalid_argument("get_dim_uint: argument is not an index"); + ex dim = ex_to<idx>(e).get_dim(); + if (!dim.info(info_flags::posint)) + throw std::invalid_argument("get_dim_uint: dimension of index should be a positive integer"); + unsigned d = ex_to<numeric>(dim).to_int(); + return d; +} + ex clifford_unit(const ex & mu, const ex & metr, unsigned char rl) { static ex unit = (new cliffordunit)->setflag(status_flags::dynallocated); @@ -734,10 +745,10 @@ ex clifford_unit(const ex & mu, const ex & metr, unsigned char rl) static idx xi((new symbol)->setflag(status_flags::dynallocated), n), chi((new symbol)->setflag(status_flags::dynallocated), n); - if ((n == M.cols()) && (n == ex_to<idx>(mu).get_dim())) { + if ((n == M.cols()) && (n == get_dim_uint(mu))) { for (unsigned i = 0; i < n; i++) { for (unsigned j = i+1; j < n; j++) { - if (M(i, j) != M(j, i)) { + if (!M(i, j).is_equal(M(j, i))) { symmetric = false; } } @@ -1170,7 +1181,7 @@ ex lst_to_clifford(const ex & v, const ex & e) { ex mu = e.op(1); ex mu_toggle = is_a<varidx>(mu) ? ex_to<varidx>(mu).toggle_variance() : mu; - unsigned dim = (ex_to<numeric>(ex_to<idx>(mu).get_dim())).to_int(); + unsigned dim = get_dim_uint(mu); if (is_a<matrix>(v)) { if (ex_to<matrix>(v).cols() > ex_to<matrix>(v).rows()) { -- 1.5.6 Best regards, Alexei -- All science is either physics or stamp collecting.
Hi, Alexei Sheplyakov schrieb:
Frankly, I've got annoyed by that evil `relational' class. It doesn't do anything useful and forces everyone to write the code in a very unnatural, counter-intuitive way. I plan to remove it in the next version of GiNaC (1.5).
could you elaborate a little bit more on your plans to remove the evil from the world? It sounds like a big API change with no apparent benefits for me. Regards, Jens
Dear Jens, On Tue, Jul 15, 2008 at 09:01:01AM +0200, Jens Vollinga wrote:
It sounds like a big API change with no apparent benefits for me.
The benefits are 1. One can compare expressions in a natural way (with operator==, operator==, etc) instead of GiNaC-specific is_equal(). 2. (a consequence of 1) No bugs like this: ex e1, e2; // do some calculation if (e1 == e2) { // gotcha: this is always true. // do something } 3. (a consequence of 1) STL containers and algorithms "just work" with GiNaC classes. The changes to achive this are not really big. Best regards, Alexei -- All science is either physics or stamp collecting.
Dear Alexei, Alexei Sheplyakov schrieb:
1. One can compare expressions in a natural way (with operator==, operator==, etc) instead of GiNaC-specific is_equal(). 2. (a consequence of 1) No bugs like this: 3. (a consequence of 1) STL containers and algorithms "just work" with GiNaC classes.
okay, so far the advantages. I like them. But then, how would one write the equalities to be given to solve(), series(), subs(), ...? I can only imagine giving two arguments to each function instead of one, but how does it work with lists? Double nested lists? Also, the answers returned by GiNaC like x==0 need to be given differently.
The changes to achive this are not really big.
Mmmh, probably not. Regards, Jens
On Tue, Jul 15, 2008 at 10:51:18AM +0200, Jens Vollinga wrote:
Dear Alexei,
Alexei Sheplyakov schrieb:
1. One can compare expressions in a natural way (with operator==, operator==, etc) instead of GiNaC-specific is_equal(). 2. (a consequence of 1) No bugs like this: 3. (a consequence of 1) STL containers and algorithms "just work" with GiNaC classes.
okay, so far the advantages. I like them. But then, how would one write the equalities to be given to solve(),
What about associative arrays? I.e. exmap eqns; eqns[a*x + b*y] = 5; eqns[c*x + d*y] = 7; exvector vars = boost::list_of(x)(y); exmap result = solve(eqns, vars);
subs(), ...?
This one is even more simple, because subs() already accepts an exmap.
series(),
Just pass 2 arguments (the variable and the expansion point).
Also, the answers returned by GiNaC like x==0 need to be given differently.
I think an exmap is a natural form for such answers. Best regards, Alexei -- All science is either physics or stamp collecting.
Dear Alexei, Alexei Sheplyakov schrieb:
What about associative arrays? I.e.
that's a very good idea. But I would wrap them in some kind of GiNaC class then like lst does with STL::list.
exmap eqns; eqns[a*x + b*y] = 5; eqns[c*x + d*y] = 7; exvector vars = boost::list_of(x)(y); exmap result = solve(eqns, vars);
But from your example (using exvector and the boost lib) I get the impression you don't like that kind of wrapping ... ? Anyway, it looks like a good improvement, at least as long as we don't care about inequalities ;-). Regards, Jens
Dear Jens, On Tue, Jul 15, 2008 at 01:16:00PM +0200, Jens Vollinga wrote:
What about associative arrays? I.e.
that's a very good idea. But I would wrap them in some kind of GiNaC class then like lst does with STL::list.
I don't quite understand why it's necessary to wrap every data type into a GiNaC class (the same applies to lst, actually). Such wrapping results in a nonsense like {x, y} + {x, y} 2*{x, y} and doesn't seem to do anything useful.
But from your example (using exvector and the boost lib) I get the impression you don't like that kind of wrapping ... ?
Exactly. I think "everything is an expression" is a bad idea. Best regards, Alexei -- All science is either physics or stamp collecting.
Dear Alexei, Alexei Sheplyakov wrote:
2. (a consequence of 1) No bugs like this:
ex e1, e2; // do some calculation if (e1 == e2) { // gotcha: this is always true. // do something }
Well, over here I get: $ cat relationaltest.cc #include <iostream> #include <ginac/ginac.h> using namespace std; using namespace GiNaC; int main() { symbol x("x"); ex e1, e2; // do some calculation e1 = series(tgamma(x),x==0,3).op(3); e2 = series(log(x),x==1,6).op(4); cout << e1 << endl; cout << e2 << endl; if (e1 == e2) { // gotcha: what happens? cout << "e1 and e2 appear to be equal." << endl; } else { cout << "e1 and e2 appear to differ." << endl; } } $ g++ -O relationaltest.cc -lginac $ ./a.out -1/12*x^2*(2*Euler^3+Euler*Pi^2+4*zeta(3)) 1/5*(-1+x)^5 e1 and e2 appear to differ. What is your point? Best wishes -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>
Dear Alexei, I wrote:
Alexei Sheplyakov wrote:
2. (a consequence of 1) No bugs like this:
ex e1, e2; // do some calculation if (e1 == e2) { // gotcha: this is always true. // do something }
Well, over here I get:
$ cat relationaltest.cc #include <iostream> #include <ginac/ginac.h> using namespace std; using namespace GiNaC;
int main() { symbol x("x"); ex e1, e2; // do some calculation e1 = series(tgamma(x),x==0,3).op(3); e2 = series(log(x),x==1,6).op(4); cout << e1 << endl; cout << e2 << endl; if (e1 == e2) { // gotcha: what happens? cout << "e1 and e2 appear to be equal." << endl; } else { cout << "e1 and e2 appear to differ." << endl; } } $ g++ -O relationaltest.cc -lginac $ ./a.out -1/12*x^2*(2*Euler^3+Euler*Pi^2+4*zeta(3)) 1/5*(-1+x)^5 e1 and e2 appear to differ.
What is your point?
I'm reading my batch of email backwards. I've now seen <http://www.ginac.de/pipermail/ginac-devel/2008-July/001370.html>, but when I try to complete these examples, they don't work as you describe: $ cat relationaltest1.cc #include <iostream> #include <ginac/ginac.h> using namespace std; using namespace GiNaC; int main() { unsigned n = 2; symbol A("A"); symbol mu_sym("mu"); idx mu(mu_sym,4); if (n == ex_to<idx>(mu).get_dim()) { cout << "doing something" << endl; } else { cout << "doing nothing" << endl; } } $ g++ -O relationaltest1.cc -lginac $ ./a.out doing nothing $ cat relationaltest2.cc #include <iostream> #include <ginac/ginac.h> using namespace std; using namespace GiNaC; int main() { matrix M(4,4); for (unsigned r=0; r<4; ++r) { for (unsigned c=0; c<4; ++c) { M.set(r,c,numeric(1,r+c+1)); } } if (M(1,2) != M(2,1)) { cout << "M(1,2) and M(2,1) appear to differ." << endl; } else { cout << "M(1,2) and M(2,1) appear to be equal." << endl; } if (M(1,3) != M(0,0)) { cout << "M(1,3) and M(0,0) appear to differ." << endl; } else { cout << "M(1,3) and M(0,0) appear to be equal." << endl; } } $ g++ -O relationaltest2.cc -lginac $ ./a.out M(1,2) and M(2,1) appear to be equal. M(1,3) and M(0,0) appear to differ. Cheers -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>
Dear Richard, On Tue, Jul 15, 2008 at 09:41:06PM +0200, Richard B. Kreckel wrote:
Well, over here I get:
$ cat relationaltest.cc #include <iostream> #include <ginac/ginac.h> using namespace std; using namespace GiNaC;
int main() { symbol x("x"); ex e1, e2; // do some calculation e1 = series(tgamma(x),x==0,3).op(3); e2 = series(log(x),x==1,6).op(4); cout << e1 << endl; cout << e2 << endl; if (e1 == e2) { // gotcha: what happens? cout << "e1 and e2 appear to be equal." << endl; } else { cout << "e1 and e2 appear to differ." << endl; } } $ g++ -O relationaltest.cc -lginac $ ./a.out -1/12*x^2*(2*Euler^3+Euler*Pi^2+4*zeta(3)) 1/5*(-1+x)^5 e1 and e2 appear to differ.
I'm a bit confused. The program fails to compile if I comment out declarations of operator[=!<>][=]*, like this diff --git a/ginac/operators.h b/ginac/operators.h index c08713f..a821aed 100644 --- a/ginac/operators.h +++ b/ginac/operators.h @@ -74,12 +74,14 @@ const numeric operator++(numeric & lh, int); const numeric operator--(numeric & lh, int); // binary relational operators ex with ex +/* const relational operator==(const ex & lh, const ex & rh); const relational operator!=(const ex & lh, const ex & rh); const relational operator<(const ex & lh, const ex & rh); const relational operator<=(const ex & lh, const ex & rh); const relational operator>(const ex & lh, const ex & rh); const relational operator>=(const ex & lh, const ex & rh); +*/ // input/output stream operators std::ostream & operator<<(std::ostream & os, const ex & e); So, e1 == e2 definitely creates the `relational' object. But it seems to give correct results anyway. I don't quite understand what's going on.
What is your point?
I've noticed that my program 1. *sometimes* gives incorrect result 2. *always* consumes more memory if I use operator== to compare expressions instead of .is_equal() method. Unfortunately, I'm unable to provide a simple test case. Also some (most?) templates libraries don't work with GiNaC::ex "out of box" because comparision operators return `relational' instead of bool. Best regards, Alexei -- All science is either physics or stamp collecting.
Dear Alexei, Alexei Sheplyakov wrote:
So, e1 == e2 definitely creates the `relational' object. But it seems to give correct results anyway. I don't quite understand what's going on.
Initially, class relational had an operator bool(), which evaluated to true or false when such an object was used in a boolean context. With GiNaC-1.1.0, that was replaced by operator safe_bool(), a technique inspired by Boost. Best wishes -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>
participants (3)
-
Alexei Sheplyakov
-
Jens Vollinga
-
Richard B. Kreckel