On 2/9/07, Sheplyakov Alexei <varg@theor.jinr.ru> wrote:

On Fri, Feb 09, 2007 at 05:39:50PM +0100, Martin Sandve Alnæs wrote: This bug is unlikely to be fixed in GiNaC 1.3.X :( There are two possible work-arounds:

1) catch std::range_error exception, e.g.

ex m; try { m = s.evalm(); } catch (std::range_error& oops) { m = s; }

This is not exactly correct (one might got std::range_error for some "real" reason)...

... which is why I can't do this, since it is in a very general part of the code (see below).

2) convert the expression from power series to polynomial:

ex m = series_to_poly(s).evalm();

This one is not exaclty fool-proof: it will fail if s is not a power series object...

... which is not acceptable either.

...

Of course, it doesn't make much sense to call evalm() on this expression, but it is a problem in an application where we want to call evalm() "just in case" on all ex objects that pass a certain point in the application.

I'd say that such an "just in a case" code is buggy and/or very inefficient in first place.

Maybe I should explain my goal. This is for the swiginac python bindings, in output typemaps. Swiginac unwraps the underlying ginac object from the ex when returning an ex to python. This is necessary to make type specific functions like matrix::transpose() available in python. But when using operators (+, * etc) on matrix objects in python, this code will fail: a = matrix1 + matrix2 b = a.transpose() Since a will be an add object, not a matrix. We don't want the user to have to call evalm() everywhere, since it destroys the high level notation we strive to reach. To fix this, we need to always call evalm() on the ex output argument before determining the underlying type. If we can't do that, we can't use the operators, which will also destroy the high level notation. Maybe there's another way around this, but I don't see it right now.

...

Thus a "do nothing" behaviour for evalm() is of great importance.

I agree -- at least it should be consistent with behaviour of other eval* methods.

Thanks for the patch. But it might be problematic to require the users to have a patched or development version of ginac, so I'll have to discuss it with other developers. If this fix can't be part of ginac 1.3.x, when will 1.4 become the stable branch? martin

On 2/12/07, Chris Dams <Chris.Dams@mi.infn.it> wrote:

Dear Martin,

On Sun, 11 Feb 2007, [ISO-8859-1] Martin Sandve Alnæs wrote:

...

ex m = series_to_poly(s).evalm();

Can't you simply call your own C++-written tree transversal algorithm at this point to check whether there actually are matrix products somewhere and then only call evalm to evaluate them?

Hi Chris, that's a good idea. I'll try that. Otherwise, I think I can live with a try-catch, and use the original expression if evalm() fails. Thanks for the ideas everyone, I'll let you know when I've worked it out. martin

Hi, On Mon, Feb 12, 2007 at 10:01:14AM +0100, Chris Dams wrote:

What about the following code?

matrix m(2,2); m = exp(x), 0, 0, sin(x); ex f = m*m; f = f.series(x,2); f = evalm(f);

For series of matrix-valued functions my patch is obviously incorrect.

I'd say that it would be best to make .map() also work with power series and to make it apply the map function to the coefficients of the power series. In that case I would also like to have pseries::eval check for zero coefficients and not store these. Any objections againts this idea?

First of all, this violates principle of least surprise. For all other classes map() operates on each term, why pseries should be different? Secondly, when operating on coefficients one often needs to know corresponding degrees too. How do we provide access to degrees? Finally, a map function could transform terms in such a way that result is not a power series any more, e.g. a*x^n => a*exp(exp(...(exp(x)))) /* n times */. What then? Best regards, Alexei -- All science is either physics or stamp collecting.

after all the difference between a expansion series and a polynomial is basically that the series contains, at the end, some information about the order of the error (i.e. O(x^n) ) In general, the error term, can not be evaluated to anything meaningful (for x!=0) (but have some properties, for example can be derived respect to x, or that O(x) - O(x) = O(x) -hey, just for curiosity: can something like this be possibly ever implemented?-). Because of the existence of the error term, in my opinion evalm must either fail or evaluate to something that holds the error term in some kind of new type (e.g. error_order type). If to fail (i.e. throw) is the choice, the 'principle of less surprise' can be honored by raising a meaningful exception (with some explanation) I must admit that I have almost no experience with ginac and that I don't fully understand the philosophy of the current developers, but, in my opinion, making every method available (or evaluable) for every type is not the way to solve this kind of problems because in many cases the possible returns could be ambiguous depending on the interpretation and this is more dangerous than throwing exception and tell the user to be more specific in the syntax. Regards, Alfredo On 2/12/07, Chris Dams <Chris.Dams@mi.infn.it> wrote:

Dear Alexei and others,

On Mon, 12 Feb 2007, Sheplyakov Alexei wrote:

...

First of all, this violates principle of least surprise. For all other classes map() operates on each term, why pseries should be different?

Okay, then I think it is best to create a method pseries::evalm and have evalm act on every coefficient in the power series. The variable and expansion piont shoulnt' be matrix value anyway, as far as I know. Any objections against that?

Best wishes, Chris

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

Dear Alexei, On Mon, 12 Feb 2007, Sheplyakov Alexei wrote:

So, here is a second attempt:

[PATCH] Declare and implement pseries::evalm() which evalm()s each coefficient.

I committed this patch now after adding a check that the coefficients aren't zero. We don't want to have zero as a coefficient in a power series. Best wishes, Chris

my question is this: do you know, if it is possible to implement the above procedure in C++? and how? maybe using some macros? Because it is certainly non standard, that you construct an instance of some class Add(x,x) and it returns an instance of a completely different class Mul(2,x).

Well, one solution is to define all the C++ classes as Add_, Mul_, etc. and then have a macro (or a function) Add, Mul, that will construct an instance of Add_ and then return Add_.eval(). But it's not very elegant.... Ondrej

I gather that you can already tell us something about the speed of the pure-Python SymPy compared to GiNaC. We are so curious!

Well, it's slow. For example expanding (x+y)**10 takes 2s, (x+y)**20 takes 20s etc. In GiNaC if I remember correctly, you can easily expand (x+y)**300 or more in a moment. Python generally is around 2000 times slower than C/C++, but I think I am using some bad algorithm in expand(). On the other hand, I have a code for calculating a Riemann and Ricci tensor from a given metric in 4D, I tried that on the Schwarzschild solution of Einstein equations and it produces the correct differential equations in about 2s (and it can be done even faster). Doing this by hand would take me 3 hours. The metric tensor contains 2 unknown functions, I don't know if GiNaC can do it as well.

In the end, it all boils down to the question how much you are willing to do statically (at compile time) as opposed to dynamically (at run time). If you want the degree of flexibility of GiNaC or a language like Python that features duck typing, where x-x can evaluate to a simple number at run time, then I don't see how macros, templates, covariant return types, and all the tricks you can read about in books on C++ can ever be helpful. Think about the eval() trick as GiNaC's way of doing duck typing in the C++ language! Please drop me a note if you happen to know a more straightforward way.

I think all the simplification is happening at runtime in both GiNaC and SymPy. But in GiNaC you are concerned with speed, that's why it's not easy to say which solution is the best. For example I believe that GiNaC is slower than FORM, but GiNaC is certainly more flexible. In SymPy I decided I want to have the simplest code possible.

I may be very biased, but I see two ways for you: Stick with Python, live happy, and don't be concerned about an occasional performance drop. Or reimplement (a possibly better) GiNaC and do Python bindings (again). But implementing the core of SymPy in C++ just for the sake of speed is probably not worth worrying about.

That is of course correct. I also implemented symbolic limits, and SymPy can do almost all of them below 1s. Currently we are just playing with the interface and algorithms, so it certainly doesn't make sense to rewrite it in C++ now. I want to do it in the future though, but only in exactly the same simple way as in SymPy (otherwise I could use GiNaC directly) if it is possible and also because that would answer your question, if there is a more straightforward way than the GiNaC's way. I think there is - using the solution from my last email + a garbage collector in C++. This should make the code as simple as SymPy is - but it's probably going to be slower (but who knows).

But back to the original topic: I do hope that the sum-of-matrices evaluation thingie is no reason for despair. Way back, that was a deliberte decision. But maybe it was a mistake to place matrices in the type hierarchy at all. Maybe they would be less confusing if they were not subtypes of basic, managed by ex, but classes standing apart from the rest of the hierarchy. That would duplicate some code, however. From past discussions, my impression is that this idea would have some followers on this list.

Exactly. I decided not to subclass a Matrix class from the Basic from precisely these reasons. Cheers, Ondrej

At the moment, I think my own alternative would be to write such a matrix class myself. I might also handle general rank r tensors as well if I do that.

Or, since you use Python anyway, you can write in Python everything and you'll get SymPy. Ondrej

Hello, On Sun, Mar 25, 2007 at 10:30:52PM +0200, Ondrej Certik wrote:

...

I agree with you completely: it's much better to not provide evalm() for anything but matrix (frankly, I think all those eval* things are plain evil). But with current [sloppy] type system such kind of interface is hardly ever possible. :(

Hi, this is the reason why GiNaC is using ex - so that it automatically calls eval(). We had this discussion already - it's because of C++, memory handling, speed and automatic evaluation.

But in Python I recently found a very elegant solution to this (in SymPy): I also have eval() but it's called automatically right after the construtor and if it returns a different instance, then that instance is returned. Example:

e = Add(x,x)

however "e" contains Mul(2,x), because in python you can make the Add.eval() get called right after the Add constructor. Thus Add.eval() evaluates x+x to 2*x and returns Mul(2,x). So in the end Add(x,x) is in fact the same as Mul(2,x).

Anyway, Add(A,Mul(-1,A)) is going to be zero (a number). Sometimes this is just plain wrong, e.g. think about power series: (x+O(x))-(x+O(x)) should be O(x), not zero. I don't think the problem is some C++ limitation[s]. I think the problem is the sloppy, mathematically inconsistent type system of GiNaC. What is add, exactly? Is it a polynomial? A rational expression? An algebraic one? A tensor? Or what?

And in the end, you don't have to think about eval() at all. And you don't need ex (python is handling the garbage collection automatically).

Typically those general-purpose garbage collectors suck badly. Otherwise I would be using some [Free] CAS (such as Axiom or Maxima) instead of fiddling with C++ (to put it mildly, it is certainly not my favourite programming language).

But soon I would like to implement the core of SymPy in C++ for speed and I am curious, if my system is going to be faster or slower than GiNaC.

my question is this: do you know, if it is possible to implement the above procedure in C++? and how? maybe using some macros?

IFAIK this is what eval() for: some kind of "virtual constructor". The constructor can not return an object of arbitrary type, but eval() can return anything.

Because it is certainly non standard, that you construct an instance of some class Add(x,x) and it returns an instance of a completely different class Mul(2,x).

This is certainly impossible in C++. But the question is: may be "virtual constructors" are not necessary at all? Suppose there is the Polynomial<Ring> type, and Polynomial<Ring3> operator+(const Polynomial<Ring2>&, const Polynomial<Ring1>&); Polynomial<Ring3> operator*(const Polynomial<Ring2>&, const Polynomial<Ring1>&); Polynomial<Ring3> operator*(const Polynomial<Ring2>&, const Ring1&); etc. Thus the return types are known at the compile time: if x is polynomial, then x + x is polynomial too. Thus there is no need for "virtual constructors". And compiler can catch more mistakes: Polynomial<Z> A = x; // OK Polynomial<Z> B = x + 0.5; // ERROR Best regards, Alexei -- All science is either physics or stamp collecting.

Hi! Sheplyakov Alexei wrote:

...

Because it is certainly non standard, that you construct an instance of some class Add(x,x) and it returns an instance of a completely different class Mul(2,x).

This is certainly impossible in C++. But the question is: may be "virtual constructors" are not necessary at all? Suppose there is the Polynomial<Ring> type, and

Polynomial<Ring3> operator+(const Polynomial<Ring2>&, const Polynomial<Ring1>&); Polynomial<Ring3> operator*(const Polynomial<Ring2>&, const Polynomial<Ring1>&); Polynomial<Ring3> operator*(const Polynomial<Ring2>&, const Ring1&); etc.

I suppose you meant Ring1 == Ring2 == Ring3.

Thus the return types are known at the compile time: if x is polynomial, then x + x is polynomial too. Thus there is no need for "virtual constructors". And compiler can catch more mistakes:

Polynomial<Z> A = x; // OK Polynomial<Z> B = x + 0.5; // ERROR

Right, but all this is not too helpful when there is no way of knowing at compile time what the ring is. After all, modular rings are different types if they have different moduli and knowing the modulus at compile time is not always possible. Think of public-key crypto algorithms where moduli have typically around 1000 bits and are determined dynamically by RNGs followed by a probabilistic primality test. Did you know that CLN's cl_univpoly_ring class is doing more or less what you suggest? It is delegating consistency checks like the above to runtime, however. And it is somwehat too complicated, in my opinion. Alas, there is no way out of this: Parametric types in C++ are determined at compile time. Just like all types in C++. Best -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>