- GiNaC-list - ginac.de

GiNaC as the symbolic manipulation engine in Sage
by Burcin Erocal 21 Aug '08

21 Aug '08

Hello, This e-mail is going to two lists (sage-devel and ginac-list), so there will be some redundant information. I hope you still take the time to read it to the end, where you can see some cool benchmarks. :) William often states that Sage is now capable of fully supporting his research. Unfortunately, for my research in symbolic summation & integration Sage lacks some of the very basic facilities. If Sage is to become a "viable alternative to Maple and Mathematica", it needs to be able to support research in symbolic computation, and development of new algorithms. At a bare minimum this means providing support for fast symbolic manipulation and arithmetic as well as some form of pattern matching. Most people reading this e-mail are probably aware of the shortcomings of the symbolics code in Sage at the moment, but here are some timings to make the situation clear. MMA: In[1]:= a = (x+y)^2; In[2]:= Timing[Sum[a^i,{i,1,10000}];] Out[2]= {0.028001, Null} Maple: > a:=(x+y)^2: > st:=time(): t:=add(a^i,i=1..10000): time()-st; 0.024 Sage via the maxima interface didn't finish this task in < 5 minutes. Timings for some of the alternatives that were considered as a backend for symbolics: Sympy (version 0.6.0 in Sage-3.0.6) didn't complete the above operation in < 5 minutes. sympycore (SVN revision 1047): sage: from sympycore import Symbol sage: x, y = Symbol("x"), Symbol("y") sage: a = (x+y)^2 sage: %time t = sum( [a^i for i in xrange(10000)] ) CPU times: user 3.45 s, sys: 0.78 s, total: 4.23 s Wall time: 4.23 s Note that even with its highly optimized structures, sympycore is far from reaching Maple/MMA performance from pure Python. At the ACA'08 conference, I've met some physicists using GiNaC (http://www.ginac.de/) a C++ library that provides symbolic manipulation and pattern matching facilities. Its lastest release was on Apr 04 2008, it also has extensive documentation, and active mailing lists. I wrote a basic Cython wrapper to benchmark its performance. Here are the numbers: sage: var("x y z",ns=1) (x, y, z) sage: a = (x+y)^2 sage: %time t = sum( [a^i for i in xrange(10000)] ) CPU times: user 1.83 s, sys: 0.00 s, total: 1.83 s Wall time: 1.83 s And some more: sage: b = (y+z)^2 sage: %time t = sum( [a^i*b^j for i in xrange(1,200,5) for j in xrange (1, 300, 3)] ) CPU times: user 0.35 s, sys: 0.00 s, total: 0.35 s Wall time: 0.36 s Maple: > a := (x+y)^2: > b := (y+z)^2: > st := time(): S := 0: for i from 1 to \ > 200 by 5 do for j from 1 to 300 by 3 d\ > o S := S+a^i*b^j end do end do: time() - st; 3.464 MMA: In[1]:= a=(x+y)^2; In[2]:= b=(y+z)^2; In[3]:= S = 0; In[4]:= Timing[ Do[ Do[ S=S+a^i*b^j, {j, 1, 300, 3}], {i, 1, 200, 5}]; ] Out[4]= {11.0527, Null} Thinking that these timings might be off because I don't know the proper way to do this in Maple or MMA, here is a different python timing: sage: s = a.parent(0) sage: st = cputime() sage: for i in xrange(1,200,5): ....: for j in xrange(1,300,3): ....: s += a^i*b^j ....: sage: et = cputime() sage: et-st 0.42002600000000001 Like many other subsystems in Sage, the symbolic arithmetic / manipulation backend should consist of a c/c++ library with a wrapper which allows the user to implement more advanced algorithms such as substitutions/tree traversals in Sage/Cython. One can of course implement the arithmetic in pure C from Cython, but this goes very much against the Sage philosophy of not reinventing the wheel. I suggest GiNaC should be the library we use in Sage for this backend. Its performance is very impressive and it can get us to a competitive position in the short term. Once we have a stable symbolics interface, we can work on more advanced algorithms such as limits / integration and improving the backend at the same time. For those who want to try my wrapper, you need to install these two spkgs (in order): http://sage.math.washington.edu/home/burcin/spkg/cln-1.2.2.spkg http://sage.math.washington.edu/home/burcin/spkg/ginac-1.4.3.spkg And apply this patch: http://www.risc.jku.at/people/berocal/sage/ginac-wrapper.patch While GiNaC seems to be the best option for Sage, it has some problems (apart from the awkward capitalization of the name which is a pain to type). - It takes ages to build The packages above took ~25 minutes to build on my desktop machine (15 for cln, 9 for ginac) - cln seems to support on only gcc, which might be a problem for the windows port (which will be using ms compilers) - cln duplicates functionality which is provided by different libraries already distributed with Sage (mpfr, mpir, flint) at a considerable speed penalty. - Creating symbolic functions at runtime is not supported. I can probably work around this in the wrapper. I've done something similar in my SCrypt package which provides symbolic computation capabilities over characteristic 2 for experimenting in crypto. It relies on PolyBoRi for arithmetic, which is just a polynomial library, i.e. knows nothing about symbolics. It would still be better to have support for runtime creation of symbolic functions in GiNaC. - IIRC, GiNaC documentation suggests the printing order for the variables is stable between different sessions, regardless of creation order. However, running the doctests in the experimental wrapper reveals different results. We may need to write a custom pretty printer. - The gethash() function of GiNaC objects doesn't return stable results, which is probably the cause of the problem above. We need a stablehash() function. I think all these problems can be resolved with some effort, and even the speed of basic arithmetic in GiNaC can be improved with modifying it's data structures. Admittedly some of these tasks are much easier than others, but supporting an already stable and mature library like GiNaC is a much better approach than trying to write our own. Thoughts? Comments? Cheers, Burcin

4 12

Fw: [sage-devel] Re: GiNaC as the symbolic manipulation engine in Sage
by Burcin Erocal 11 Aug '08

11 Aug '08

Can I get some help with these questions? I don't see a link to a bug tracker on the web site. Is there a list of known bugs? I guess the test code is placed under the directory check/ in the source tree? Thanks. Burcin Begin forwarded message: Date: Mon, 11 Aug 2008 08:24:52 -0700 From: "William Stein" <wstein(a)gmail.com> To: sage-devel(a)googlegroups.com Subject: [sage-devel] Re: GiNaC as the symbolic manipulation engine in Sage Burcin, 1. does Ginac (or the part we plan to use) have *any* known bugs at all? 2. does Ginac have a test suite? 3. does ginac have a bug tracker? 4. what does Valgrind say about the quality of the code when run on the test suite? memory leaks? -- William

2 1

Re: [GiNaC-list] [sage-devel] Re: GiNaC as the symbolic manipulation engine in Sage
by Burcin Erocal 11 Aug '08

11 Aug '08

Hi, I forwarded your questions to the ginac list. I'll try to answer as much as I know. On Mon, 11 Aug 2008 08:24:52 -0700 "William Stein" <wstein(a)gmail.com> wrote: <snip> > Burcin, > > 1. does Ginac (or the part we plan to use) have *any* known bugs at all? I can't find a public bug tracker on the web site. They suggest sending bugs to the mailing list. > 2. does Ginac have a test suite? Some of the files under the check/ directory seems to verify correctness. > 3. does ginac have a bug tracker? > 4. what does Valgrind say about the quality of the code when run on > the test suite? memory leaks? I haven't run it through valgrind at all. Maybe mabshoff can lend a helping hand here? Burcin

1 0

Re: [GiNaC-list] [sage-devel] Re: GiNaC as the symbolic manipulation engine in Sage
by Burcin Erocal 11 Aug '08

11 Aug '08

On Mon, 11 Aug 2008 02:37:40 -0700 "William Stein" <wstein(a)gmail.com> wrote: > > On Mon, Aug 11, 2008 at 1:30 AM, Burcin Erocal <burcin(a)erocal.org> wrote: <snip> > > > > At the ACA'08 conference, I've met some physicists using GiNaC > > (http://www.ginac.de/) a C++ library that provides symbolic > > manipulation and pattern matching facilities. Its lastest release was > > on Apr 04 2008, it also has extensive documentation, and active mailing > > lists. I wrote a basic Cython wrapper to benchmark its performance. > > Here are the numbers: > > > > sage: var("x y z",ns=1) > > (x, y, z) > > sage: a = (x+y)^2 > > sage: %time t = sum( [a^i for i in xrange(10000)] ) > > CPU times: user 1.83 s, sys: 0.00 s, total: 1.83 s > > Wall time: 1.83 s > > This benchmark might be lame. Even the current "insanely slow" > maxima-based symbolics in Sage seem to do better than what you > just wrote above: > > sage: var('x,y') > (x, y) > sage: a = (x+y)^2 > sage: time t = sum(a^i for i in xrange(10000)) > CPU times: user 0.51 s, sys: 0.03 s, total: 0.54 s > Wall time: 0.55 s > > and > > sage: var('y,z'); b = (y+z)^2 > sage: time t = sum([a^i*b^j for i in xrange(1,200,5) for j in > xrange(1, 300, 3)] ) > CPU times: user 0.30 s, sys: 0.01 s, total: 0.31 s > Wall time: 0.32 s > > That's because the benchmark is maybe testing nothing > but memory allocation. Your benchmark doesn't call maxima at all. When you quit Sage after these lines, there is no Exiting spawned Maxima process. notice printed. > I have seen a lot of *real* benchmarks involving symbolic calculus > when sad users and students come to me and ask why their code > is 1000 times slower than Mathematica. This happened a *few* times > today at Sage Days 9 today. It also comes up in sage-support all > the time. I really need to come up with a large list of specific examples > like this, since I'm getting way too confused by benchmarks. I couldn't find any standard benchmarks, and had to come up with the one above. I think it is reflects some aspect of performance, as the power operator is very similar to a symbolic function application, and the system has to do a basic simplification like ((x+y)^2)^i = (x+y)^(2*i) for each step. Cheers, Burcin

1 0

Re: [GiNaC-list] [sage-devel] Re: GiNaC as the symbolic manipulation engine in Sage
by Burcin Erocal 11 Aug '08

11 Aug '08

> On Mon, Aug 11, 2008 at 1:30 AM, Burcin Erocal <burcin(a)erocal.org> wrote: <snip> > > While GiNaC seems to be the best option for Sage, it has some problems > > (apart from the awkward capitalization of the name which is a pain to > > type). > > > > - It takes ages to build > > The packages above took ~25 minutes to build on my desktop > > machine (15 for cln, 9 for ginac) > > > > - cln seems to support on only gcc, which might be a problem for the > > windows port (which will be using ms compilers) > > > > - cln duplicates functionality which is provided by different libraries > > already distributed with Sage (mpfr, mpir, flint) at a considerable > > speed penalty. On Mon, 11 Aug 2008 11:39:42 +0200 Jens Vollinga <jensv(a)nikhef.nl> wrote: > It would be quite useful for ginac if the numeric class that wraps the > cln library would allow for other libraries as well, even pure C double > precision ones. The problem in realizing this idea was always that other > libraries were missing certain features (some functions, modular > arithmetic,...) and to make up for these would incur a lot of extra > work. It didn't seem like it was worth the hassle (yet). I also think that using other libraries instead of cln in the (subclasses of) numeric class is the way to go. With William's remarks below, this becomes essential for Sage. :) On Mon, 11 Aug 2008 02:37:40 -0700 "William Stein" <wstein(a)gmail.com> wrote: > > The single big question that needs to be answered that you don't > answer in your post is the following: > > Is it likely in your opinion that with a month of hard work by you (say), > would it be possible to create something based on Ginac that would > actually satisfy the following: > > 1. Does not depend on cln. > 2. Builds in less than 8 minutes. > > I mean this as a 100% technical questions (not social, etc.). > I just had a look at the source code directory for ginac, and it is about > 20,000 lines of pretty well documented C++ code total. I have > the strong impression that the use of cln is nearly completely > encapsulated in the file numeric.cpp, which is about 2000 lines. > I'm guessing a Sage version of Ginac could simply replace numeric.cpp > by our own implementation. Of course this is possible. I am willing to put in most of the work to make this happen. However, support from the GiNaC developers for such an effort is very important. They know their code base best, with their suggestions this could be accomplished much faster. > > cln will not be included standard in Sage ever. > I could definitely see a version of ginac that doesn't depend on cln > as possible for Sage... > > -- William> Cheers, Burcin

1 0

website down?
by Burcin Erocal 11 Aug '08

11 Aug '08

Hello, The GiNaC web site seems to be down at the moment. Does anyone know when it will be back up? Thanks. Burcin

2 1

bug with is_polynomial?
by Jonathan Cross 09 Aug '08

09 Aug '08

Hello, I'm using GiNaC 1.4.3 built on Mac OS X. I think I've found a problem with is_polynomial when powers are involved. When there is an expression involving a non-integer power---even when the expression is independent of the dummy variable of the polynomial--- then the expression is not considered a polynomial by GiNaC. The following code illustrates the problem: (or perhaps I don't fully understand "is_polynomial"?) ------------------------- #include <iostream> #include <ginac/ginac.h> using namespace std; int main() { GiNaC::symbol x("x"); GiNaC::symbol s("s"); // here are two expressions that are both polynomials with respect to s. GiNaC::ex expr1 = sin(x) + 2*s; GiNaC::ex expr2 = pow(2,x) + 2*s; // works as expected if (expr1.is_polynomial(s)) cout << expr1 << " is a polynomial with respect to " << s << endl; else cout << expr1 << " is NOT a polynomial with respect to " << s << endl; // reports 2^x + 2s is not a polynomial in s? if (expr2.is_polynomial(s)) cout << expr2 << " is a polynomial with respect to " << s << endl; else cout << expr2 << " is NOT a polynomial with respect to " << s << endl; return 0; } ------------------------ I'm not familiar enough with the internal workings of GiNaC to make a patch, but this seems like the changes required would be fairly small. Thanks!

2 4

Re: [GiNaC-list] cl_F_idecode.cc line 30 ?
by Chris Bouchard 30 Jul '08

30 Jul '08

I've since been able to archive an expression whose length is ~1200 lines... I'm at a loss. The only difference between the expressions that generate the error and the ones that don't seems irrelevant, but here it is. Based on the code fragment: ... ex c = TwoLoop2Pt(...); // c.nops() = 2 ex numeric = c.op(0); //returned for numerical integration ex analytic = c.op(1); //returned for analytic use, f(eps,...) ex ana_series = analytic.series(eps==0,2); ex ana_0 = ana_series.coeff(eps,0); ex ana_1 = ana_series.coeff(eps,-1); ex ana_2 = ana_series.coeff(eps,-2); ... numeric doesn't work (~1200 lines) ana_0 doesn't work (~1200 lines) ana_1 works (~1200 lines) ana_2 works (~600 lines) I guess there's either something fundamentally different between ana_0 and ana_1 (seems unlikely since they come from the same expression, c.op(1)) or there's a critical expression length right around 1200 lines. Perhaps some knowledge about cl_F_idecode.cc would help? Sincerely, Chris Bouchard ---- Original message ---- >Date: Tue, 29 Jul 2008 16:08:22 -0500 (CDT) >From: Chris Bouchard <cbouchrd(a)illinois.edu> >Subject: [GiNaC-list] cl_F_idecode.cc line 30 ? >To: GiNaC discussion list <ginac-list(a)ginac.de> > >GiNaC mailing list, > >I've run across an error with a notification request. The >error occurs while trying to archive an expression. The >error and request are: > >"Internal error: statement in file >./float/misc/cl_F_idecode.cc, line 30 has been reached!! >Please send the authors of the program a description how you >produced this error!" > >The file cl_F_idecode.cc is in the cln library and a snippet >is given here, with a comment indicating line 30: > >... >namespace cln { > >const cl_idecoded_float integer_decode_float (const cl_F& x) >{ > floatcase(x > , return integer_decode_float(x); > , return integer_decode_float(x); > , return integer_decode_float(x); > , return integer_decode_float(x); > ); // this is line 30 >} >} > >After playing around with my archive code I'm beginning to >think the error occurs when the expression being archived >becomes too long. I haven't nailed down a critical >expression length but 27 lines works while ~1200 lines >generated the error. > >I can't think of anything else that's changing. > >The expressions are individual terms from series expansions >of expressions returned from xloops function calls (eg. >TwoLoop2Pt(...)) and don't contain any noticeably offensive >pieces. > >Sincerely, >Chris Bouchard >_______________________________________________ >GiNaC-list mailing list >GiNaC-list(a)ginac.de >https://www.cebix.net/mailman/listinfo/ginac-list

2 2

cl_F_idecode.cc line 30 ?
by Chris Bouchard 29 Jul '08

29 Jul '08

GiNaC mailing list, I've run across an error with a notification request. The error occurs while trying to archive an expression. The error and request are: "Internal error: statement in file ./float/misc/cl_F_idecode.cc, line 30 has been reached!! Please send the authors of the program a description how you produced this error!" The file cl_F_idecode.cc is in the cln library and a snippet is given here, with a comment indicating line 30: ... namespace cln { const cl_idecoded_float integer_decode_float (const cl_F& x) { floatcase(x , return integer_decode_float(x); , return integer_decode_float(x); , return integer_decode_float(x); , return integer_decode_float(x); ); // this is line 30 } } After playing around with my archive code I'm beginning to think the error occurs when the expression being archived becomes too long. I haven't nailed down a critical expression length but 27 lines works while ~1200 lines generated the error. I can't think of anything else that's changing. The expressions are individual terms from series expansions of expressions returned from xloops function calls (eg. TwoLoop2Pt(...)) and don't contain any noticeably offensive pieces. Sincerely, Chris Bouchard

1 0

Re: [GiNaC-list] comparing symbols
by Chris Bouchard 25 Jul '08

25 Jul '08

Hi, >try the code > >cout << tree << "x = " << x << "\n"; > >then we can see what is going on. x = lst @0x8a08c70, hash=0xc, flags=0x13, nops=2 0 (symbol) @0x8a3df40, serial=62607, hash=0x771ff4, flags=0x7, domain=0 1 (symbol) @0x8a2bde8, serial=62608, hash=0x771ff4, flags=0x7, domain=0 >Maybe r_sum is called with something like > symbol f1("0"); > symbol f2("1"); > r_sum(...,lst(f1,f2),...); Wouldn't that eliminate functional dependence on the variables represented by f1 & f2, making the return value a function of "0" and "1". The end user could handle this by manually simplifying the output, but that wouldn't be optimal. Sincerely, Chris Bouchard

2 1