May 2010 - GiNaC-list

Term ordering and compiling C++ code
by Doug 02 Jun '10

02 Jun '10

Hi, I'm running into a term ordering problem like the one mentioned in this post. http://permalink.gmane.org/gmane.comp.mathematics.ginac.general/1107 Relevant details about my application are below, for the curious. The short of it, though, is that the C++ code I am generating with ginac is so large (50MB) that the compiler can't handle it, and moreover, the values being represented by ginac symbols are array accesses that the compiler can't optimize away. So I need to create temporary values and do sub-expression elimination on the C++ code just so the compiler can have a chance. I have successfully done this using search-and-replace on the C++ code for other applications, but I guess I was just lucky that ginac printed terms in a consistent order. Is this a problem other ginac users have confronted? From the answer to the above post, it sounds like I need to write my own print function that can ensure expressions are printed in a predictable order. It looks like I can also inspect the expression tree and possibly do a proper job of common sub-expression elimination. But if there is an easier approach, I'd be very grateful for a pointer or two. Thanks, Doug Baker ------------------- Background Details ------------------- My application is related to visual odometry. I am computing the gradient and hessian for a function of 10 variables that involves converting quaternions to rotation matrices, applying those matrices to 3D points and projecting those points onto a 2D plane. I am generating code with ginac that is initially 50MB of C++ code. (And thank you for ginac being able to handle that monstrosity with such grace.) The input symbols are all array elements (e.g. symbol x("state[t+X]")) or function calls (e.g. symbol rTi11("rTi.getValue<1,1>()")). If I take the printed output from ginac as is, the C++ compiler can't properly optimize it because it doesn't know that state[t+X] is a constant or that getValue is constant. Moreover, the expressions have huge repeated expressions, and these bog the compiler down, too. I've been able to use a simple search-and-replace approach in other problems, and on this one that same approach has gotten me from 50MB to 6MB or so, but the term ordering issue is biting me now. For example, I have one sum (a+b+c+d) that occurs 232608 times. But there are 4! = 24 permutations of that sum that I'd have to replace. It's worse than that, though, because that's not the only expression like it. Some are bigger, like this: (-2.0*sK/sKSIJ+4.0*sK*sS2/pKSIJ2+-4.0*sI*sJ*sS/pKSIJ2). This has something like 72 permutations, not even counting the permutations of pKSIJ and pKSIJ2, which are themselves expressions like (a+b+c+d). And not only are there variations on this due to term ordering, similar expressions occur with different variables (e.g. with sJ instead of sK, etc.) So when all these expressions can be output in various orders, it's no longer feasible to generate all the permutations for search-and-replace. Support NPR 20 seconds at a time. www.twentysecondsatatime.org

5 18

Re: [GiNaC-list] Term ordering and compiling C++ code
by Martin Ettl 28 May '10

28 May '10

Hello all, i have extended the GinacPrint about following functions: /** * With this function the default variable name can be changed. * By default it is "e"; Passing an empty string resets * to default. * * @param newVarName the new variable name * */ void setVarName(const std::string &newVarName); /** * This functions returns the current variable name. * By default it is "e"; * * @return the current variable name * */ const std::string & getVarName() const; /** * Clears all internal variables. This is a reset to default function. * */ void clear(); I have attached the modified version. Best regards from Germany Martin -------- Original-Nachricht -------- > Datum: Wed, 26 May 2010 17:55:19 -0700 (PDT) > Von: Doug <cape1232(a)yahoo.com> > An: GiNaC discussion list <ginac-list(a)ginac.de> > Betreff: Re: [GiNaC-list] Term ordering and compiling C++ code > Beautiful. Thanks! > > > --- On Wed, 5/26/10, Martin Ettl <ettl.martin(a)gmx.de> wrote: > > From: Martin Ettl <ettl.martin(a)gmx.de> > > ... > What do you think about my modifications (see attachment)? > > > > > -- GRATIS für alle GMX-Mitglieder: Die maxdome Movie-FLAT! Jetzt freischalten unter http://portal.gmx.net/de/go/maxdome01

3 3

Re: [GiNaC-list] Term ordering and compiling C++ code
by Martin Ettl 27 May '10

27 May '10

Hello Doug and all, thank you for sharing your code. I think this should be included to ginac release. The ability to create c++-code is very imporant. I modified your version in order to support c++-templates. This means you can define the datatype, e.g.: T, long double, some custom type,.... e.g.: for long double type use the function: gprint.setType("long double"); By default the double type is used, in order to not brake backwardcompatiblity. What do you think about my modifications (see attachment)? Best regards from Germany Martin Ettl -------- Original-Nachricht -------- > Datum: Tue, 25 May 2010 07:45:05 -0700 (PDT) > Von: Doug <cape1232(a)yahoo.com> > An: GiNaC discussion list <ginac-list(a)ginac.de> > Betreff: Re: [GiNaC-list] Term ordering and compiling C++ code > Hi All, > > Attached is the code for generating C/C++ code from a GiNaC expression > that pulls out subexpressions and declares them up front so that they only > have to be evaluated once. The .hh file has documentation and a simple usage > example. > > I hope this is useful. > > I used it to reduce 50MB of code to under 300K and improved the running > time dramatically. > > -Doug > > Support NPR 20 seconds at a time. www.twentysecondsatatime.org > > > > > -- GRATIS für alle GMX-Mitglieder: Die maxdome Movie-FLAT! Jetzt freischalten unter http://portal.gmx.net/de/go/maxdome01

2 1

Re: [GiNaC-list] Term ordering and compiling C++ code
by Doug 26 May '10

26 May '10

Hi All, Attached is the code for generating C/C++ code from a GiNaC expression that pulls out subexpressions and declares them up front so that they only have to be evaluated once. The .hh file has documentation and a simple usage example. I hope this is useful. I used it to reduce 50MB of code to under 300K and improved the running time dramatically. -Doug Support NPR 20 seconds at a time. www.twentysecondsatatime.org

2 1

input compiling sample
by Martin Ettl 26 May '10

26 May '10

Hello, i failed to compile the sample from :http://www.ginac.de/tutorial/Input_002foutput.html#Input_002foutput on Ubuntu linux with gcc-4.4.0. #include <iostream> #include <string> #include <stdexcept> #include <ginac/ginac.h> using namespace std; using namespace GiNaC; int main() { cout << "Enter an expression containing 'x': " << flush; parser reader; try { ex e = reader(cin); symtab table = reader.get_syms(); symbol x = table.find("x") != table.end() ? ex_to<symbol>(table["x"]) : symbol("x"); cout << "The derivative of " << e << " with respect to x is "; cout << e.diff(x) << "." << endl; } catch (exception &p) { cerr << p.what() << endl; } } Do you have any suggestions, i used the following command line to compile the sample: g++ -o test test1.cpp -lginac $ g++ -o test test1.cpp -lginac test1.cpp: In function »int main()«: test1.cpp:11: Error: »parser« not declared in this scope test1.cpp:11: Error: expected »;« before »reader« test1.cpp:14: Error: »reader« not declared in this scope test1.cpp:15: Error: »symtab« not declared in this scope test1.cpp:15: Error: expected »;« before »table« test1.cpp:16: Error: »table« not declared in this scopedefiniert It looks like i missed some #include <ginac/parser.h> ??? Thanks in advance Martin -- GRATIS für alle GMX-Mitglieder: Die maxdome Movie-FLAT! Jetzt freischalten unter http://portal.gmx.net/de/go/maxdome01

1 1

Re: [GiNaC-list] Term ordering and compiling C++ code
by Doug 24 May '10

24 May '10

kreckel wrote: ...and second, sharing is currently not pursued aggressively! Rather, it is exploited whenever it is trivial to do so. (The product rule of differentiation is an example where exactly the same terms pop up again and again so exploiting sharing comes at no cost.) I think the benefit would still be large, even if coverage isn't 100%. See below. Well, I think that if the size of generated code is so prohibitively large and compiler CSE doesn't help you may be better off writing your own algorithm collecting similar terms in an associative array. You could then artificially introduce temporaries, in order to help the compiler. This would boil down to a more aggressive variant of GiNaC's subexpression sharing. What do you think? This is exactly what I have done. As mentioned, I will post that solution here in the next few days. This works well, but is very slow, as it takes something that is sub-linear in the size of the expression tree and makes it O(n^2). If GiNaC expression tree nodes had a unique identifier, this approach could achieve sub-linear performance. Specifically, it would be linear in the size of the unique subexpressions, not the size of the actual expression. (If you have a lot of sharing, your expression can be large, but the number of unique expressions can be small.) Even if GiNaC doesn't identify every single instance of a shared sub-expression, the benefits could still be large. In fact, I tried an experiment to test this. I treated ex::gethash() as a unique identifier, even though it isn't. The code I generated was wrong, of course, but it was very fast to generate it, and very fast to compile and execute. Thus, whatever number of unique expressions were not actually identified as such was small enough that this approach would still be quite useful. Also, FYI, when using the hash value as a unique identifier, I only got a few collisions out of thousands of sub-expressions. That is, I only found a few expressions with the same hash even though they were different expressions. Thanks, Doug

2 1

Can GiNaC set info_flag?
by Seungyun Cho 19 May '10

19 May '10

For example, I want to set a symbol r to be rational, or real, but without specifying a value. I know one can use info(info_flags::real) to determine if an object is real, but I want to do the opposite. I found a type 'realsymbol' to do this, but I want many other scopes such as rational, prime, positive integer, pure imaginary, etc. Is there any solution?

2 1

[PATCH] found a memory leak in /ginac/excompiler.cpp
by Martin Ettl 17 May '10

17 May '10

Hello, during a longer debugging session i identified a memory leak in ginac/excompiler.cpp at line 116 and 121. Please refer the attached patch that fixes the issue. Best regards Ettl Martin -- GRATIS für alle GMX-Mitglieder: Die maxdome Movie-FLAT! Jetzt freischalten unter http://portal.gmx.net/de/go/maxdome01

3 4

conjugates of power objects
by Burcin Erocal 06 May '10

06 May '10

On Mon, 26 Apr 2010 18:02:40 -0700 (PDT) kcrisman <kcrisman(a)gmail.com> wrote: > On Apr 26, 4:09 pm, John Cremona <john.crem...(a)gmail.com> wrote: > > This is certainly a bug: > > > > sage: a = sqrt(-3) > > sage: a > > sqrt(-3) > > sage: a.conjugate() > > sqrt(-3) > > > > sage: bool(a==a.conjugate()) > > True > > > > Yeah, this is bad, and according to a.conjugate?? it seems to be > something wrong with our Ginac/Pynac use. This is a bug in GiNaC: ginsh - GiNaC Interactive Shell (ginac V1.5.7) __, _______ Copyright (C) 1999-2010 Johannes Gutenberg University Mainz, (__) * | Germany. This is free software with ABSOLUTELY NO WARRANTY. ._) i N a C | You are welcome to redistribute it under certain conditions. <-------------' For details type `warranty;'. Type ?? for a list of help topics. > sqrt(-3); sqrt(-3) > conjugate(sqrt(-3)); sqrt(-3) For conjugation, power objects just compute the conjugate of the basis and the exponent, and construct a new power object from these. Here is the relevant function: http://pynac.sagemath.org/hg/file/3ece9ba22005/ginac/power.cpp#l805 Cheers, Burcin

4 9