GiNaC as the symbolic manipulation engine in Sage
Burcin Erocal
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
William Stein
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.
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.
> 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.
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.
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
Burcin Erocal
> > 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 <je...@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" <wst...@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
Burcin Erocal
"William Stein" <wst...@gmail.com> wrote:
>
> On Mon, Aug 11, 2008 at 1:30 AM, Burcin Erocal <bur...@erocal.org> wrote:
<snip>
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
Burcin Erocal
Jens Vollinga <je...@nikhef.nl> wrote:
> Hi from GiNaC side,
>
> just some short comments:
>
> Burcin Erocal schrieb:
> > - It takes ages to build
> > The packages above took ~25 minutes to build on my desktop
> > machine (15 for cln, 9 for ginac)
>
> Did you configure with the --disable-static option? Otherwise your build
> time is just twice the usual time (looks like it).
I didn't do anything to try to make these times better. I was more
concerned with benchmarking speed of basic manipulations, and thought
that these problems can be solved easily later.
Your suggestion reduces the compilation time of ginac to 5 minutes.
Looks like we're already below William's 8 minute limit. Now can we fix
the cln problem? :)
Thanks.
Burcin
> > - cln seems to support on only gcc, which might be a problem for the
> > windows port (which will be using ms compilers)
>
> <blame others>
> That is a real problem, but it can only be solved by the cln people.
> </blame others>
>
> > - cln duplicates functionality which is provided by different libraries
> > already distributed with Sage (mpfr, mpir, flint) at a considerable
> > speed penalty.
>
> 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).
>
> > - 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.
>
> True.
>
> > - 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.
>
> Yes, the two points above are related and I am to a certain point
> responsible for that mess. To explain: The original tinfo system in
> ginac was changed a few years ago to avoid fixed magic numbers which are
> a problem if you have a lot of user supplied ginac classes (made more
> urgent by the plan to make functions into full ginac classes). The
> current system assigns these dynamically during compilation. But it is
> not done in a smart way, so the tinfo numbers usually differ between
> different compilations. Now, these tinfo numbers also go into the hash
> calculation and the term ordering ...
>
> Personally, I would like to have this issue solved and have ginac to
> produce a fixed term ordering (as long as the declaration order of the
> symbols done by the user is the same). But at the moment I don't have
> the time to do this. Any volunteers? :-)
>
> Regards,
> Jens
> _______________________________________________
> GiNaC-list mailing list
> GiNaC...@ginac.de
> https://www.cebix.net/mailman/listinfo/ginac-list
>
Mike Hansen
I'm not going to comment on benchmarks since I don't really know how
Gary's symbolics branch does, but I will focus on what affects me the
most.
I've done a bit of work fixing bugs and such in the current calculus
code, and it's an area of the Sage codebase I am pretty familiar with.
I find it very unlikely that I would work GiNaC directly; although,
I'd probably work on the wrapper side of things. I'm concerned with
how extensible using GiNaC would be. For example, if I wanted to add
piecewise functions or support for the MeijerG or LambertW functions
or elliptic integrals, I'd have to do so within GiNaC. Maybe someone
else would, but it isn't the most active project either from looking
at the mailing lists and changelog (
http://www.ginac.de/ginac.git?p=ginac.git;a=shortlog ). For what it's
worth, switching to GiNaC would probably remove my (however small)
contribution to the bus factor.
--Mike
Pearu Peterson
Burcin Erocal wrote:
> 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.
This is just because the test above uses Python `sum` that
is not optimal. Try this (that would be equivalent to using `add`
in Maple):
from sympycore import *
x, y = Symbol("x"), Symbol("y")
a = (x+y)**2
%time t = Calculus.Add( *[a**i for i in xrange(10000)] )
In my computer, using `Calculus.Add` is 110 times faster than using
Python `sum`. Compared to Maple, sympycore would be only 10x slower.
Pearu
Burcin Erocal
I don't see why you have to implement these in GiNaC directly. What's
wrong with implementing such things on top of the (hypothetical)
wrapper in Sage? (Isn't this already done this way in
sage/functions/piecewise.py?)
I am only proposing to use GiNaC for basic arithmetic / symbolic
manipulation, people developing symbolic algorithms should be able to
interact with the data structures only using python/cython.
We should use other specialized code in Sage to cover most
functionality such as linear system solving, gcds and factorization
anyway, since Sage has access to much faster specialized libraries for
these.
Note that the ginac developers are very responsive to questions and
requests. I received a response to my initial e-mail to ginac-list in
69 minutes, which is certainly comparable to William's response time of
67 minutes. :)
Cheers,
Burcin
David Joyner
>
> Hi Mike,
>
> On Mon, 11 Aug 2008 04:21:02 -0700
> "Mike Hansen" <mha...@gmail.com> wrote:
>
>>
>> Hi Burcin,
>>
>> I'm not going to comment on benchmarks since I don't really know how
>> Gary's symbolics branch does, but I will focus on what affects me the
>> most.
>>
>> I've done a bit of work fixing bugs and such in the current calculus
>> code, and it's an area of the Sage codebase I am pretty familiar with.
>> I find it very unlikely that I would work GiNaC directly; although,
>> I'd probably work on the wrapper side of things. I'm concerned with
>> how extensible using GiNaC would be. For example, if I wanted to add
>> piecewise functions or support for the MeijerG or LambertW functions
>> or elliptic integrals, I'd have to do so within GiNaC. Maybe someone
>> else would, but it isn't the most active project either from looking
>> at the mailing lists and changelog (
>> http://www.ginac.de/ginac.git?p=ginac.git;a=shortlog ). For what it's
>> worth, switching to GiNaC would probably remove my (however small)
>> contribution to the bus factor.
>
> I don't see why you have to implement these in GiNaC directly. What's
> wrong with implementing such things on top of the (hypothetical)
> wrapper in Sage? (Isn't this already done this way in
> sage/functions/piecewise.py?)
I wrote most (or maybe all) of that module and would not relish the thought
of swapping out maxima for ginac (assuming ginac has the functionality needed).
I like the philosophy (as I understand it) William Stein has taken
with regard to
sympy. Namely, make it part of Sage and then gradually add necessary
base classes to support a possible transition, consider on a
case-by-case basis
which maxima calls might be more efficiently done with a sympy call.
Perhaps something similar could be done with ginac? Would it help your
research if a "method=maxima/sympy/ginac" option was added to some functions?
Burcin Erocal
I just read the file sage/functions/piecewise.py. The calls to maxima
are used for integration, differentiation, and trigsimplify. GiNaC does
not provide these functions. The "laplace" function seems to use some
symbolic manipulation, leaving it as it is wouldn't hurt anyone, but
changing it might make it much faster.
I suggest using GiNaC for symbolic arithmetic only. This would mean
that the pause you experience on the last line of this
sage: var("x y z")
(x, y, z)
sage: a = (x+y)
sage: b = (y+z)
sage: a*b
(y + x)*(z + y)
will go away. We will also gain the ability to write our own
simplification, integration, etc. routines as a bonus.
> I like the philosophy (as I understand it) William Stein has taken
> with regard to
> sympy. Namely, make it part of Sage and then gradually add necessary
> base classes to support a possible transition, consider on a
> case-by-case basis
> which maxima calls might be more efficiently done with a sympy call.
> Perhaps something similar could be done with ginac? Would it help your
> research if a "method=maxima/sympy/ginac" option was added to some functions?
My example patch adds an optional parameter to "var", which allows one
to create symbolic variables backed by GiNaC
I am mainly concerned with the speed of arithmetic, and support for
pattern matching.
> > I am only proposing to use GiNaC for basic arithmetic / symbolic
> > manipulation, people developing symbolic algorithms should be able to
> > interact with the data structures only using python/cython.
> >
> > We should use other specialized code in Sage to cover most
> > functionality such as linear system solving, gcds and factorization
> > anyway, since Sage has access to much faster specialized libraries for
> > these.
<snip>
Burcin
Mike Hansen
>> > wrong with implementing such things on top of the (hypothetical)
>> > wrapper in Sage? (Isn't this already done this way in
>> > sage/functions/piecewise.py?)
>>
>>
>> I wrote most (or maybe all) of that module and would not relish the thought
>> of swapping out maxima for ginac (assuming ginac has the functionality needed).
>
> I just read the file sage/functions/piecewise.py. The calls to maxima
> are used for integration, differentiation, and trigsimplify. GiNaC does
> not provide these functions. The "laplace" function seems to use some
> symbolic manipulation, leaving it as it is wouldn't hurt anyone, but
> changing it might make it much faster.
>
> I suggest using GiNaC for symbolic arithmetic only. This would mean
> that the pause you experience on the last line of this
>
> sage: var("x y z")
> (x, y, z)
> sage: a = (x+y)
> sage: b = (y+z)
> sage: a*b
> (y + x)*(z + y)
>
> will go away. We will also gain the ability to write our own
> simplification, integration, etc. routines as a bonus.
You typically want to be able do arithmetic with piecewise functions
or expressions involving LambertW. Are you saying that we should keep
the current SymbolicArithmetic hierarchy that's in Sage and just swap
out Maxima for GiNaC?
--Mike
William Stein
>
>>> > I don't see why you have to implement these in GiNaC directly. What's
>>> > wrong with implementing such things on top of the (hypothetical)
>>> > wrapper in Sage? (Isn't this already done this way in
>>> > sage/functions/piecewise.py?)
>>>
>>>
>>> I wrote most (or maybe all) of that module and would not relish the thought
>>> of swapping out maxima for ginac (assuming ginac has the functionality needed).
>>
>> I just read the file sage/functions/piecewise.py. The calls to maxima
>> are used for integration, differentiation, and trigsimplify. GiNaC does
>> not provide these functions. The "laplace" function seems to use some
>> symbolic manipulation, leaving it as it is wouldn't hurt anyone, but
>> changing it might make it much faster.
>>
>> I suggest using GiNaC for symbolic arithmetic only. This would mean
>> that the pause you experience on the last line of this
>>
>> sage: var("x y z")
>> (x, y, z)
>> sage: a = (x+y)
>> sage: b = (y+z)
>> sage: a*b
>> (y + x)*(z + y)
>>
>> will go away. We will also gain the ability to write our own
>> simplification, integration, etc. routines as a bonus.
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
William Stein
Can ginac do manipulation with symbols like:
* (sqrt(3)-3)*(sqrt(3)+1)/6
I'm wondering, because to follow up on 1 above, evidently
Maxima has a major bug according to another sage-support
thread in just manipulating an expression like the above,
which is IMHO pretty sad, given the age of maxima.
-- William
Burcin Erocal
"William Stein" <wst...@gmail.com> wrote:
<snip>
> Can ginac do manipulation with symbols like:
>
> * (sqrt(3)-3)*(sqrt(3)+1)/6
If you installed the packages, there should be a binary called "ginsh"
in SAGE_LOCAL/bin. It lets you try ginac through a simple shell
interface.
> I'm wondering, because to follow up on 1 above, evidently
> Maxima has a major bug according to another sage-support
> thread in just manipulating an expression like the above,
> which is IMHO pretty sad, given the age of maxima.
I will look at that thread more carefully, and try the example with
ginsh.
Burcin
William Stein
>
> On Mon, 11 Aug 2008 08:28:42 -0700
> "William Stein" <wst...@gmail.com> wrote:
>
> <snip>
>> Can ginac do manipulation with symbols like:
>>
>> * (sqrt(3)-3)*(sqrt(3)+1)/6
>
> If you installed the packages, there should be a binary called "ginsh"
> in SAGE_LOCAL/bin. It lets you try ginac through a simple shell
> interface.
Wow, that's nice.
By the way, anybody with an account on sage.math can try ginac
and ginsh right now by doing:
anybody@sage:~$ sage-ginac
----------------------------------------------------------------------
| SAGE Version 3.0.6.final, Release Date: 2008-07-29 |
| Type notebook() for the GUI, and license() for information. |
----------------------------------------------------------------------
sage: var("x y z",ns=1)
(x, y, z)
sage: expand((x+y+z+1)^3)
1+3*z*y^2+3*x+3*x^2*z+3*y^2+z^3+3*z+y^3+6*z*y+3*z^2+3*x^2*y+3*x^2+3*z^2*y+6*x*y+x^3+6*x*z+3*y+3*x*z^2+6*x*z*y+3*x*y^2
sage: !ginsh
ginsh - GiNaC Interactive Shell (ginac V1.4.3)
__, _______ Copyright (C) 1999-2008 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)-3)*(sqrt(3)+1)/6
> ;
1/6*(1+sqrt(3))*(-3+sqrt(3))
>
>
>> I'm wondering, because to follow up on 1 above, evidently
>> Maxima has a major bug according to another sage-support
>> thread in just manipulating an expression like the above,
>> which is IMHO pretty sad, given the age of maxima.
>
> I will look at that thread more carefully, and try the example with
> ginsh.
>
>
> Burcin
>
> >
>
--
William Stein
Associate Professor of Mathematics
University of Washington
http://wstein.org
Burcin Erocal
Here is what ginac does. Seems right to me. :)
> t1 = (sqrt(3)-3)*(sqrt(3)+1)/6;
1/6*(1+sqrt(3))*(-3+sqrt(3))
> tt1 = -1/sqrt(3);
-1/3*sqrt(3)
> t2 = sqrt(3)/6;
1/6*sqrt(3)
> t1+t2;
1/6*(1+sqrt(3))*(-3+sqrt(3))+1/6*sqrt(3)
> tt1+t2;
-1/6*sqrt(3)
> expand(t1);
-1/3*sqrt(3)
> expand(tt1);
-1/3*sqrt(3)
> expand(t2);
1/6*sqrt(3)
> expand(t1+t2);
-1/6*sqrt(3)
> expand(tt1+t2);
-1/6*sqrt(3)
Burcin Erocal
Burcin Erocal <bur...@erocal.org> wrote:
>
> On Mon, 11 Aug 2008 11:39:42 +0200
> Jens Vollinga <je...@nikhef.nl> wrote:
>
> > Hi from GiNaC side,
> >
> > just some short comments:
> >
> > Burcin Erocal schrieb:
> > > - It takes ages to build
> > > The packages above took ~25 minutes to build on my desktop
> > > machine (15 for cln, 9 for ginac)
> >
> > Did you configure with the --disable-static option? Otherwise your build
> > time is just twice the usual time (looks like it).
>
> I didn't do anything to try to make these times better. I was more
> concerned with benchmarking speed of basic manipulations, and thought
> that these problems can be solved easily later.
>
> Your suggestion reduces the compilation time of ginac to 5 minutes.
> Looks like we're already below William's 8 minute limit. Now can we fix
> the cln problem? :)
Using the --disable-static option while configuring cln brings down its
build&install time to 9 minutes.
One thing that remains in favor of cln is that it uses memory pools.
This would be hard to do with a simple rewrite. I don't know how much
this effects performance though.
Burcin
Burcin Erocal
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" <wst...@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
William Stein
cln is absolutely not an option. Full stop. No way in hell, etc.
We have a "no way in hell" policy on any code that isn't going to easily
build under MSVC in windows for at least the next year in Sage. And
there is no way I have enough patience to wait another year to get the
symbolic situation in Sage resolved.
So just completely put any thoughts of using cln in Sage out of your mind.
William
William Stein
I just can't frickin' believe that maxima gets this basic thing wrong.
That doesn't make sense. I've posted a ticket about this...
http://trac.sagemath.org/sage_trac/ticket/3805
William
David Joyner
We use maxima 5.13. It is apparently fixed in maxima 5.15. (The
current version is
5.16.)
>
> http://trac.sagemath.org/sage_trac/ticket/3805
>
> William
>
> >
>
William Stein
In my trac ticket I explain that it seems to me that doing the above
directly in maxima at the prompt works in the version of maxima
included in sage. So the problem may be either the mode sage
gets maxima into for the above, or something involving communication
back and forth with maxima. Or I'm somehow just missing the point.
Could you replicate the bug using maxima 5.13 directly?
David Joyner
Yes, using the maxima_console on an amd64 hardy heron machine.
sage: maxima_console()
Maxima 5.13.0 http://maxima.sourceforge.net
Using Lisp CLISP 2.46 (2008-07-02)
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
This is a development version of Maxima. The function bug_report()
provides bug reporting information.
(%i1) t1 : (sqrt(3)-3)*(sqrt(3)+1)/6
;
(sqrt(3) - 3) (sqrt(3) + 1)
(%o1) ---------------------------
6
(%i2) t1 : (sqrt(3)-3)*(sqrt(3)+1)/6;
(sqrt(3) - 3) (sqrt(3) + 1)
(%o2) ---------------------------
6
(%i3) tt1 : -1/sqrt(3);
1
(%o3) - -------
sqrt(3)
(%i4) t2 : sqrt(3)/6;
- 1/2
3
(%o4) ------
2
(%i5) expand(t1+t2);
- 1/2
3 3
(%o5) - --------
2 2
(%i6) expand(tt1+t2);
1
(%o6) - ---------
2 sqrt(3)
(Sorry, I know I should have turned off print2d.)
>
>>
>>
>>>
>>> http://trac.sagemath.org/sage_trac/ticket/3805
>
> >
>
Burcin Erocal
Date: Mon, 11 Aug 2008 20:29:48 +0400
From: Alexei Sheplyakov <va...@theor.jinr.ru>
To: GiNaC discussion list <ginac...@ginac.de>
Subject: Re: [GiNaC-list] Fw: [sage-devel] Re: GiNaC as the symbolic manipulation engine in Sage
On Mon, Aug 11, 2008 at 05:38:46PM +0200, Burcin Erocal wrote:
> 1. does Ginac (or the part we plan to use) have *any* known bugs at all?
Yes, although it's purely theoretical:
http://bugs.debian.org/cgi-bin/bugreport.cgi?bug=479163
> 2. does Ginac have a test suite?
Yes, have a look at the check directory.
> 3. does ginac have a bug tracker?
No.
> 4. what does Valgrind say about the quality of the code when run on
> the test suite? memory leaks?
Valgrind reports the memory allocated for static objects as a leak. AFAIK
valgrind does so for any C++ library which has static objects (with
non-trivial constructors and destructors), so this might be a misfeature
of valgrind. Still I'm not quite sure, this might be also an actual bug.
Any ideas on debugging this are wellcome.
I routinely use GiNaC in long-running computations, so the *actual* memory
leaks get noticed and fixed (see e.g. commit
dcb1d9061e4ce0a012d145530384541fc941fef4).
Best regards,
Alexei
--
All science is either physics or stamp collecting.
Robert Dodier
> 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.
In the case of the summation problem mentioned in this thread,
the central problem is that Maxima's simplification algorithm
for the result (something of the form (x + y)^2 + (x + y)^4 + ...)
seems to run in n^2 time where n is the number of terms.
There is also a trade-off in the summation algorithm:
in this case it would be faster to generate a list of all the
terms and then apply "+" to them, but if all the terms are
numbers, the calculation could be done in constant space
by adding each term to a running total; the latter is the
current implementation. Probably the trade-off should go
the other way; I 'll look at that.
I'll look at the "+" simplification algorithm to see if I can see
a way to speed it up. Interested parties may direct
followups to the Maxima mailing list.
best
Robert Dodier
rjf
1. There are fixes to the general simplifier using hash coding that
make collection and ordering of terms run much faster.
In essence, replacing sorting that probably takes O(n^3) time by O(n)
time. I have provided these, but there are some nuances to
integrating these fixes that other people should look at.
2. The fixes are, to some extent, irrelevant for the following reason:
If the problem is large and consists of polynomial arithmetic, you
should use rat() form.
If the problem is large and is NOT obviously polynomial arithmetic,
you should make some substitutions and
convert it to polynomial form if at all possible.
My guess is that code that is 1000 times slower than Mathematica will
most likely be made 100, 1000 or even 5000 times faster by the
judicious insertion of rat() somewhere.
The particular example cited above, where (presumably) this is timed
in Maxima:
a: (x+y)^2$
sum(a^i,i,1,10000)$
takes 29 seconds on my Pentium D workstation.
Perhaps you are counting the time to compute and transmit the result
to your SAGE front end.
Tests which I find revealing are those that look like this:
f(x,n):= if (n=1) then x else g(f(x,n-1)); maybe store these in an
array; Maxima makes this easy to do:
f[x,1]:=x$
f[x,n]:=g(f[x,n-1])$
sum(f(x,i), i,1,100)$
or
sum(f[x,i],i,1,100)$
or try
h[x,i]:=rat(f[x,i])
and sum that.
I assume that there are nuances in GiNac or any other system which
will cause it to appear to misbehave if you use it in ways that are
not anticipated by the authors.
RJF
mark mcclure
> The particular example cited above, where (presumably) this is timed
> in Maxima:
> a: (x+y)^2$
> sum(a^i,i,1,10000)$
> takes 29 seconds on my Pentium D workstation.
>
> Perhaps you are counting the time to compute and transmit the result
> to your SAGE front end.
Pro:
------------------ Maxima session ------------------
mark-mcclures-computer-2:~ markmcclure$ maxima
Maxima 5.13.0 http://maxima.sourceforge.net
Using Lisp CLISP 2.41 (2006-10-13)
Dedicated to the memory of William Schelter.
This is a development version of Maxima. The function bug_report()
provides bug reporting information.
(%i2) sum(a^i,i,1,10000)$
(%i3) time(%o2);
(%o3) [291.801613]
---------------------------------------------------------------
----------------- Mathematica session ------------
mark-mcclures-computer-2:~ markmcclure$ rlwrap math
Mathematica 6.0 for Mac OS X x86 (32-bit)
Copyright 1988-2008 Wolfram Research, Inc.
In[1]:= a = (x+y)^2;
Out[2]= {0.029104, Null}
-------------------------------------------------------------
Mark McClure
Ondrej Certik
Yes. One needs to postpone all such simplifications/terms sorting
until it is absolutely necessary.
As to this thread, I would like to say --- yes, please go ahead! As
you can see, Sage has enough momentum already to force the other
packages to either improve, or they will not be used. I think it is
very healthy (and necessary!) for all the packages involved.
As to SymPy, our aim is to get as fast as maple or mathematica and
also to cover all the advanced features. We are much better than ginac
on features, but much worse on speed, so the speed has became a high
priority for me now, I am playing with different cores here:
http://groups.google.com/group/sympy/browse_thread/thread/8362465daafbd4dc
Ondrej
Robert Dodier
> 1. There are fixes to the general simplifier using hash coding that
> make collection and ordering of terms run much faster.
> In essence, replacing sorting that probably takes O(n^3) time by O(n)
> time. I have provided these, but there are some nuances to
> integrating these fixes that other people should look at.
I will take a look at merging that code.
> 2. The fixes are, to some extent, irrelevant for the following reason:
> If the problem is large and consists of polynomial arithmetic, you
> should use rat() form.
> If the problem is large and is NOT obviously polynomial arithmetic,
> you should make some substitutions and
> convert it to polynomial form if at all possible.
hand.
In any event, a naive formulation should scale in a reasonable
manner, even if there are special methods to make it faster.
FWIW
Robert Dodier
Gary Furnish
my symbolics engine explicitly supports is different evaluation
mechanisms for different types. GiNaC does not do many of the things
I had in mind for symbolics, but there is no reason symbolics could
not wrap GiNaC for Complex Symbolics. Then we could merge symbolics
sooner, and we could even use maxima if there are any things GiNaC
could not handle (or native symbolics for when both can not handle
something). Furthermore symbolics would serve as a nice abstraction
layer making it easier to switch between Maxima features and GiNaC
features and native Sage features, and would also serve to prevent the
Sage calculus API from becoming too tied into the GiNaC api. Finally,
Symbolics would also provide a nice way to wrap some of the features
of GiNaC that are symbolic but not calculus based ( such as some of
the noncommutative algebras used in physics). We would then still be
able to do things like implement differential geometry in a system
independent way on top of symbolics (while making use of the typed
coercion system). Perhaps by wrapping GiNaC with symbolics we can get
get a better system sooner (with more features then either has on its
own). Thoughts?
-- Bill
rjf
On Aug 11, 9:23 pm, mark mcclure <mcmcc...@unca.edu> wrote:
.....
> I, at least, can confirm the extreme timing discrepancy on my MacBook
> Pro:
>
> ------------------ Maxima session ------------------
> mark-mcclures-computer-2:~ markmcclure$ maxima
> Distributed under the GNU Public License. See the file COPYING.
> Dedicated to the memory of William Schelter.
> This is a development version of Maxima. The function bug_report()
> provides bug reporting information.
> (%i1) a: (x+y)^2$
> (%i2) sum(a^i,i,1,10000)$
> (%i3) time(%o2);
> (%o3) [291.801613]
> ---------------------------------------------------------------
(which is what CLISP is),
or you should run Maxima compiled to machine language.
There are several compiled versions of Maxima in lisps, in particular
GCL, SBCL.
but also ECL, Allegro CL, Scieneer, ...
Also, a few pages of change to the general simplifier can make this
much faster. (100X?)
It seems that for most people the interest in Maxima is to make it do
things that it cannot do at all,
rather than making it run faster.
It seems that the he people using Sage apparently are interested in
doing things
that can generally be done (like adding polynomials?) but faster and
larger.
rjf
On Aug 12, 8:38 am, Robert Dodier <robert.dod...@gmail.com> wrote:
> Well, after trying it, it appears rat doesn't help in the case at
> hand.
for a=(x+y)^2, it would
expand out a^2, a^3, etc. and would therefore produce, , a polynomial
starting with
(x^20000*y^20000) + ....
This is quite a different problem from the one expressed by the
mathematica and maple code, which
presumably preserves
(x+y)^20000 + (x+y)^19998+ ... an expression with a mere 10,000
terms.
and one that it is hard to view as being of much interest except to
test the internal sorting program.
Oh, by the way, the answer from Maple, if you look at, is probably
not sorted. I haven't tried it, but
Maple tends to produce answers in what appears to be random order.
Maple advocates claim that
the order is not random, even though it is not ordered by degree. It
is indeed deterministic in some way.
You determine it by setting up your Maple session the same way each
time. If you set up
your Maple session differently, the order may deterministically be
different.
RJF
mark mcclure
>
> To be fair, you should either run mathematica in an interpreter
> (which is what CLISP is), or you should run Maxima compiled
> to machine language.
major source of the problem for the original poster as, I
believe, that Sage uses CLISP. At least it does on my
MacBook:
------- Start of Maxima session via Sage --------
mark-mcclures-computer-2:~ markmcclure$ /Applications/sage-3.0.2/sage -
maxima
Maxima 5.13.0 http://maxima.sourceforge.net
I think the Sage folks are trying to switch over to ECL
and, hopefully, that will bring the time discrepancies
down.
Mark McClure
mabshoff
On Aug 12, 10:28 am, mark mcclure <mcmcc...@unca.edu> wrote:
> On Aug 12, 12:02 pm, rjf <fate...@gmail.com> wrote:
> > To be fair, you should either run mathematica in an interpreter
> > (which is what CLISP is), or you should run Maxima compiled
> > to machine language.
>
> Good to know; thanks a lot. Perhaps, in fact, this is a
> major source of the problem for the original poster as, I
> believe, that Sage uses CLISP. At least it does on my
> MacBook:
>
> ------- Start of Maxima session via Sage --------
> mark-mcclures-computer-2:~ markmcclure$ /Applications/sage-3.0.2/sage -
> maxima
> -------------------------------------------------
>
> I think the Sage folks are trying to switch over to ECL
> and, hopefully, that will bring the time discrepancies
> down.
(depending on OSX, compiler, etc somewhere around 2 to 8 times as
fast), so there are still orders of magnitude in some benchmarks. But
in the end this might not be relevant for many use cases. The issue
with Maxima is plainly that we have to use pexpect for communication,
so even if Maxima were as fast as say MMA the pexpect overhead would
still kill us, especially for the common case, i.e. computing the 60th
power of some polynomial is not the case that is killing us on a
regular basis, but adding some absolute term to some parametric plot
would slow plotting down by a couple orders of magnitude due to
pexpect. That case is a million times faster with _fast_eval in Sage,
so we cannot use anything for polynomial arithmetic via pexpect.
> Mark McClure
Cheers,
Michael
mark mcclure
>
> Yes, but even with sbcl or gcl we are taking about a constant
> (depending on OSX, compiler, etc somewhere around 2 to 8 times as
> fast), so there are still orders of magnitude in some benchmarks. But
> in the end this might not be relevant for many use cases. The issue
> with Maxima is plainly that we have to use pexpect for communication,
> so even if Maxima were as fast as say MMA the pexpect overhead would
Yes, I am aware of the pexpect issue as well and that is
definitely an issue.
There's a lot of good reasons to like Sage, though, and one
is that it is an easily accessible repository of excellent
open source software. If, for example, someone asked me for
the easiest way to install Maxima on Mac OSX, I'd tell them
to install Sage and then run install_scripts from the Sage
command line. For such a person, the pexpect issue is no
problem.
Conversely, I've installed Maxima and other open source
programs on Mac OSX a number of times and generally found
it to be a complete pain. I think the ease of Sage's
installation is in itself quite an accomplishment.
So thanks,
Mark McClure
PS
> That case is a million times faster with _fast_eval in Sage,
> so we cannot use anything for polynomial arithmetic via pexpect.
index of the reference manual. It's mainly a mystery to me.
William Stein
>
> On Aug 12, 1:49 pm, mabshoff <Michael.Absh...@fsmath.mathematik.uni-
> dortmund.de> wrote:
>>
>> Yes, but even with sbcl or gcl we are taking about a constant
>> (depending on OSX, compiler, etc somewhere around 2 to 8 times as
>> fast), so there are still orders of magnitude in some benchmarks. But
>> in the end this might not be relevant for many use cases. The issue
>> with Maxima is plainly that we have to use pexpect for communication,
>> so even if Maxima were as fast as say MMA the pexpect overhead would
>> still kill us, ...
>
> Yes, I am aware of the pexpect issue as well and that is
> definitely an issue.
>
> There's a lot of good reasons to like Sage, though, and one
> is that it is an easily accessible repository of excellent
> open source software. If, for example, someone asked me for
> the easiest way to install Maxima on Mac OSX, I'd tell them
> to install Sage and then run install_scripts from the Sage
> command line. For such a person, the pexpect issue is no
> problem.
We're not at all talking about getting rid of Maxima. We're
just talking about not using Maxima to *implement* in Sage
the core symbolic manipulation routines. That's all. I'm
sure Maxima will be included in Sage for a very long time
to come.
> Conversely, I've installed Maxima and other open source
> programs on Mac OSX a number of times and generally found
> it to be a complete pain. I think the ease of Sage's
> installation is in itself quite an accomplishment.
>
> So thanks,
> Mark McClure
>
> PS
> You wrote:
>> That case is a million times faster with _fast_eval in Sage,
>> so we cannot use anything for polynomial arithmetic via pexpect.
>
> I've seen some examples of fast_eval, but it doesn't appear in the
> index of the reference manual. It's mainly a mystery to me.
Check out the source code:
http://www.sagemath.org/hg/sage-main/file/5e2426d277eb/sage/ext/fast_eval.pyx
jlapeyre
(%i2) a : (x+y)^2;
2
(%o2) (y + x)
(%i3) ans : sum(a^i,i,1,10000)$
(%i4) time(%o3);
(%o4) [26.09]
I remember in the early 90s some reasonable call to a bessel
function caused Mma to hang. I can't recall if it hung the kernel
completely. A couple years later, there was a list of slow Mma things
on the web. One was to allocate a list of zeros (maybe 1000 or 10000),
then ask for element 500 say. This took about two minutes to return.
In the early 2000's the last example had not changed.