[PATCH] clifford_unit: fix possible bugs due to wrong operator[!=]= usage.

[Alexei Sheplyakov-2]

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.



_______________________________________________
GiNaC-devel mailing list
GiNaC-devel@...
https://www.cebix.net/mailman/listinfo/ginac-devel

[Jens Vollinga]

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
_______________________________________________
GiNaC-devel mailing list
GiNaC-devel@...
https://www.cebix.net/mailman/listinfo/ginac-devel

[Alexei Sheplyakov-2]

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.



_______________________________________________
GiNaC-devel mailing list
GiNaC-devel@...
https://www.cebix.net/mailman/listinfo/ginac-devel

[Jens Vollinga]

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
_______________________________________________
GiNaC-devel mailing list
GiNaC-devel@...
https://www.cebix.net/mailman/listinfo/ginac-devel

[Alexei Sheplyakov-2]

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.



_______________________________________________
GiNaC-devel mailing list
GiNaC-devel@...
https://www.cebix.net/mailman/listinfo/ginac-devel

[Jens Vollinga]

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

_______________________________________________
GiNaC-devel mailing list
GiNaC-devel@...
https://www.cebix.net/mailman/listinfo/ginac-devel

[Alexei Sheplyakov-2]

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.



_______________________________________________
GiNaC-devel mailing list
GiNaC-devel@...
https://www.cebix.net/mailman/listinfo/ginac-devel

[Richard B. Kreckel]

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/>
_______________________________________________
GiNaC-devel mailing list
GiNaC-devel@...
https://www.cebix.net/mailman/listinfo/ginac-devel

[Richard B. 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/>
_______________________________________________
GiNaC-devel mailing list
GiNaC-devel@...
https://www.cebix.net/mailman/listinfo/ginac-devel

[Alexei Sheplyakov-2]

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.



_______________________________________________
GiNaC-devel mailing list
GiNaC-devel@...
https://www.cebix.net/mailman/listinfo/ginac-devel

[Richard B. Kreckel]

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/>
_______________________________________________
GiNaC-devel mailing list
GiNaC-devel@...
https://www.cebix.net/mailman/listinfo/ginac-devel