programming: define a new function

23 views
Skip to first unread message

Maurizio

unread,
Apr 14, 2009, 5:35:37 PM4/14/09
to sage-support, sage-...@googlegroups.com
Hi all.

I'm willing to invest some of my time to understand if I can be able
to do a step ahead with symbolic functions.

How are special symbolic functions supposed to be defined? I am
willing to experiment with delta of dirac function. This has some
special properties (see http://en.wikipedia.org/wiki/Dirac_delta_function
), some of them are really useful but I don't know how to define them
in a CAS like maxima or SAGE.

I'm aware that it is already present in maxima, even though I don't
think it is recognized by SAGE. I am wondering whether a viable
approach could be to add to calculus.py a section similar to the one
of Function_gamma, so that SAGE simply interfaces to maxima. I don't
know if this is useful or not.

Otherwise, I would be interested in knowing if this could be done with
the new symbolic package. Burcin proved to be very helpful in showing
me a simple way to define delta function by means of its values and he
assigned it as being the derivative of heaviside function (defined in
a sort of piecewise function):

sage: heaviside(x).diff(x)
dirac(x)

is there a way to implement the other properties? I am willing to know
if is there any documentation about that, because I am not able to
find that!

I am willing to learn something about pynac, but please feel free to
discourage me if you think it is too far away from being ready. Is
there any integration or derivation capability ready? Is it possible
to start testing it using maxima's integration capabilities? (I don't
think so...)
I was browsing the todo page ( http://wiki.sagemath.org/symbolics/pynac_todo
) but it seems that many action items went away... are they already
accomplished? (what about the TODO showed in http://wiki.sagemath.org/symbolics
?)

I have to say I find the actual SAGE documentation seems pretty hard
to browse, but could be my fault. Yesterday I lost half an hour
looking for a "numerical solve" or something like that, before finding
the "find_roots" function. Today, I spent half an hour looking for
"differential equation solve" with no success. I am sure I could do
some DE solution in the past (something with maxima, something with
SymPy, I think, all through SAGE), but I think that I found the way to
do it pretty easily browsing the old reference manual...

Thanks a lot

Maurizio

PS: delta of dirac is already in SymPy (
http://code.google.com/p/sympy/issues/detail?id=672&can=1&q=dirac ),
am I correct in thinking that the current function definition is
different than SAGE? In this case, I assume this could be some good
reading, but not necessarily a source of inspiration, right?

Burcin Erocal

unread,
Apr 15, 2009, 4:32:46 AM4/15/09
to sage-...@googlegroups.com
Hi Maurizio,

On Tue, 14 Apr 2009 14:35:37 -0700 (PDT)
Maurizio <maurizio...@gmail.com> wrote:

>
> Hi all.
>
> I'm willing to invest some of my time to understand if I can be able
> to do a step ahead with symbolic functions.

Great to hear that.

> How are special symbolic functions supposed to be defined?

They will be defined based on the SFunction class in
sage/symbolics/function.pyx. There is no real documentation in that
file, just some examples to test if the wrapper is working as intended.

There is some more documentation in the GiNaC tutorial here:

http://www.ginac.de/tutorial/Symbolic-functions.html

You can do (almost) all that you see in the tutorial, without using any
C++. I say almost, because I still did not set up an equivalent
of .hold(). That should come soon given the current interest in pynac.

Mike Hansen is working on switching Sage's symbolics over to pynac
completely. We talked about changing the design of SFunction on IRC, he
might have some code for you to build on.

Even if you start now, rebasing your code to the new design should be
trivial. You'll just have to put the various functions you define in a
class, which subclasses SFunction.

> I am
> willing to experiment with delta of dirac function. This has some
> special properties (see
> http://en.wikipedia.org/wiki/Dirac_delta_function ), some of them are
> really useful but I don't know how to define them in a CAS like
> maxima or SAGE.
>
> I'm aware that it is already present in maxima, even though I don't
> think it is recognized by SAGE. I am wondering whether a viable
> approach could be to add to calculus.py a section similar to the one
> of Function_gamma, so that SAGE simply interfaces to maxima. I don't
> know if this is useful or not.

With more and more people getting interested in the switch, I think it
has a real chance of happening very soon. I refrained from saying this
before, but at this moment, I don't think it's wise to invest more time
in the old symbolics code.

> Otherwise, I would be interested in knowing if this could be done with
> the new symbolic package. Burcin proved to be very helpful in showing
> me a simple way to define delta function by means of its values and he
> assigned it as being the derivative of heaviside function (defined in
> a sort of piecewise function):
>
> sage: heaviside(x).diff(x)
> dirac(x)
>
> is there a way to implement the other properties? I am willing to know
> if is there any documentation about that, because I am not able to
> find that!

For reference, the code I sent Maurizio with an example containing
heaviside and dirac functions is here:

http://sage.math.washington.edu/home/burcin/pynac/dirac.py

Implementing heaviside is #2452 on trac.

> I am willing to learn something about pynac, but please feel free to
> discourage me if you think it is too far away from being ready. Is
> there any integration or derivation capability ready? Is it possible
> to start testing it using maxima's integration capabilities? (I don't
> think so...)

Derivation is easy, and it's already done by pynac. We will be using
maxima for integration for a while now, though I plan to start
submitting native integration (and summation) code to Sage once the
switch is complete.


For people interested in helping out, a few words about integration is
in order. As we discussed on this list before, many problems in
integration can be handled with heuristics before calling more advanced
algorithms. The heuristics also give more user friendly answers in most
cases.

Integration code in Sage should run through some preprocessing before
calling anything else, be it maxima, or a native implementation of the
Risch algorithm.

There are some links to descriptions of heuristics used by maxima and maple at this page:

http://wiki.sagemath.org/symbolics

Pattern matching and rewriting capabilities of the new symbolics will
be very useful for this. Anybody interested in improving the
integration capabilities of Sage can help, no high level mathematical
knowledge is necessary.

> I was browsing the todo page
> ( http://wiki.sagemath.org/symbolics/pynac_todo ) but it seems that
> many action items went away... are they already accomplished?

The items in the todo list didn't go away, they were copied to the
"done" section, with a link to the trac ticket containing the changes.

There are two pynac packages waiting for review on trac. Since they
both contain substantial changes, it would be good if somebody other
than the developer took them for a ride. :)

I only know two people who tried to use pynac for anything real until
now, Alex Raichev and Nick Alexander. Their comments have been really
helpful and motivating.

> (what
> about the TODO showed in http://wiki.sagemath.org/symbolics ?)


This is a more long term todo list, both to remember what still needs
to be done and for reference if anybody is interested in helping out.
Most tasks on this list are highly nontrivial, and need substantial
work to be complete.

Feel free to edit the list and add items you think are important. For
example, you can add z-transforms and laplace transforms as we
discussed before.


Problems related to special functions used to be classified under the
component calculus on trac. I am not sure if I should add special
functions stuff to that page, or start a new one. We have quite a few
special functions related bugs on trac. There was also some effort to
get Sage on the list of software in the DLMF, but the wiki page for
this seems to be deleted. (mabshoff?)

> I have to say I find the actual SAGE documentation seems pretty hard
> to browse, but could be my fault. Yesterday I lost half an hour
> looking for a "numerical solve" or something like that, before finding
> the "find_roots" function. Today, I spent half an hour looking for
> "differential equation solve" with no success. I am sure I could do
> some DE solution in the past (something with maxima, something with
> SymPy, I think, all through SAGE), but I think that I found the way to
> do it pretty easily browsing the old reference manual...

sage: deso<tab>
desolve desolve_laplace desolve_system desolvers

Especially "desolvers?" seems helpful.

> Thanks a lot

Thank you.


Cheers,
Burcin

mabshoff

unread,
Apr 15, 2009, 4:39:15 AM4/15/09
to sage-devel


On Apr 15, 1:32 am, Burcin Erocal <bur...@erocal.org> wrote:

<SNIP>

Hi Burcin,

> There was also some effort to
> get Sage on the list of software in the DLMF, but the wiki page for
> this seems to be deleted. (mabshoff?)

I checked old backups and the directory I keep spam pages around and I
cannot find it. So either the conversion a while ago caused the
trouble or something went terribly wrong. Sorry.

Cheers,

Michael

Burcin Erocal

unread,
Apr 15, 2009, 4:58:39 AM4/15/09
to sage-...@googlegroups.com
Hi Michael,

No problem. We didn't lose anything much, since the "effort" didn't get
anywhere. The page Robert Bradshaw started contained a simple header,
and a few section headings without any real content. It won't be hard
to reproduce it, if there is any time that is...

BTW, could we get the subscriptions for the wiki working again?

Thanks.

Burcin

mabshoff

unread,
Apr 15, 2009, 5:23:22 AM4/15/09
to sage-devel


On Apr 15, 1:58 am, Burcin Erocal <bur...@erocal.org> wrote:
> Hi Michael,

<SNIP>

Hi Burcin,

> > I checked old backups and the directory I keep spam pages around and I
> > cannot find it. So either the conversion a while ago caused the
> > trouble or something went terribly wrong. Sorry.
>
> No problem. We didn't lose anything much, since the "effort" didn't get
> anywhere. The page Robert Bradshaw started contained a simple header,
> and a few section headings without any real content. It won't be hard
> to reproduce it, if there is any time that is...

Ok, I understand :)

> BTW, could we get the subscriptions for the wiki working again?

Well, that has annoyed me for a long, long time, but I have never
found the time to fix this. I will raise the issue tomorrow (well,
today) at UW's Sage status meeting and maybe someone will get
motivated to look into it.

> Thanks.
>
> Burcin

Cheers,

Michael

Maurizio

unread,
Apr 16, 2009, 6:43:41 PM4/16/09
to sage-devel
> For people interested in helping out, a few words about integration is
> in order. As we discussed on this list before, many problems in
> integration can be handled with heuristics before calling more advanced
> algorithms. The heuristics also give more user friendly answers in most
> cases.
>
> Integration code in Sage should run through some preprocessing before
> calling anything else, be it maxima, or a native implementation of the
> Risch algorithm.
>
> There are some links to descriptions of heuristics used by maxima and maple at this page:
>
> http://wiki.sagemath.org/symbolics
>
> Pattern matching and rewriting capabilities of the new symbolics will
> be very useful for this. Anybody interested in improving the
> integration capabilities of Sage can help, no high level mathematical
> knowledge is necessary.
>

This is very interesting, indeed I have no high level mathematical
knowledge :)

Just for fun, I was trying to implement the most basic examples (1.4
and 1.5) of this page, linked in the wiki:
http://www.cs.uwaterloo.ca/~kogeddes/papers/Integration/IntSurvey.html

Let me make this straight: I have not understood what is intended to
be a resultant() and I have absolutely no experience with rings nor
with multivariate polynomials. Nonetheless, I tried this code in SAGE:

reset()
P.<x,y,z> = PolynomialRing(QQ)
# P.<x,y,z> = GF(5)[]

A = -1
B = x^3 + x
tores = A - z*B.derivative(x)
res = tores.resultant(B,x); factor(res)

I get:
(-1) * (z + 1) * (2*z - 1)^2

That is somehow close to the result shown in the link, but why is it
not identical?

Moreover, why if I use the second statement (the one in the comment)
to define the ring, do I get a completely different answer? E.g.:
(z + 1) * (z + 2)^2

Which is the correct way to define the ring in this case?

Finally, even assuming that I can get the right answer from this,
which is the recommended way to get the roots of an equation given by
a "univariate polynomials == 0"? This is supposed to be the next step
of the algorithm.

The worst way I can think of is to use repr(), so that I can use solve
(), but hopefully there's a better way, at least to convert a
polynomials to a symbolic expression... Is there?

I know this is not leading anywhere... but thank you for your time! :)

Regards

Maurizio

Burcin Erocal

unread,
Apr 17, 2009, 8:40:44 AM4/17/09
to sage-...@googlegroups.com

A is 1 in example 1.5.

> B = x^3 + x
> tores = A - z*B.derivative(x)
> res = tores.resultant(B,x); factor(res)
>
> I get:
> (-1) * (z + 1) * (2*z - 1)^2

With A = 1:

sage: P.<x,z> = QQ[]
sage: b = x^3+x
sage: b.resultant(1-z*b.derivative(x),x)
-4*z^3 + 3*z + 1
sage: factor(_)
(-1) * (z - 1) * (2*z + 1)^2

which matches the maple result.

> That is somehow close to the result shown in the link, but why is it
> not identical?
>
> Moreover, why if I use the second statement (the one in the comment)
> to define the ring, do I get a completely different answer? E.g.:
> (z + 1) * (z + 2)^2
>
> Which is the correct way to define the ring in this case?

The theory only works over characteristic 0, i.e., your fields should contain QQ. Also note that,

sage: P.<x,z> = GF(5)[]
sage: (x+x^5).derivative(x)
1

What is the integral of 1 w.r.t x in this case?

> Finally, even assuming that I can get the right answer from this,
> which is the recommended way to get the roots of an equation given by
> a "univariate polynomials == 0"? This is supposed to be the next step
> of the algorithm.
>
> The worst way I can think of is to use repr(), so that I can use solve
> (), but hopefully there's a better way, at least to convert a
> polynomials to a symbolic expression... Is there?

You can get this information from the factorization you compute above.


Cheers,
Burcin

Carl Witty

unread,
Apr 18, 2009, 10:44:56 AM4/18/09
to sage-...@googlegroups.com
On Thu, Apr 16, 2009 at 3:43 PM, Maurizio <maurizio...@gmail.com> wrote:
> Finally, even assuming that I can get the right answer from this,
> which is the recommended way to get the roots of an equation given by
> a "univariate polynomials == 0"? This is supposed to be the next step
> of the algorithm.

Taking a quick look at that page, it looks like they want the exact
roots in CC of a polynomial with algebraic coefficients. In Sage, we
can get this with QQbar:

sage: K.<x> = QQbar[]
sage: (x^5-x-1).roots(ring=QQbar)

[(1.167303978261419?, 1),
(-0.7648844336005847? - 0.3524715460317263?*I, 1),
(-0.7648844336005847? + 0.3524715460317263?*I, 1),
(0.1812324444698754? - 1.083954101317711?*I, 1),
(0.1812324444698754? + 1.083954101317711?*I, 1)]

(AFAIK, maxima can't do this; I don't think maxima can handle general
algebraic numbers. Since Sage's solve() is implemented using maxima,
solve() won't work for this problem.)

Carl

Maurizio

unread,
Apr 18, 2009, 5:46:19 PM4/18/09
to sage-devel
Carl, Burcin,

thank you very much for your support.

Burcin, I'm sorry for the trivial mistake. Thank you for pointing it
out.
Unfortunately, I don't understand this:
"
The theory only works over characteristic 0, i.e., your fields should
contain QQ. Also note that,

sage: P.<x,z> = GF(5)[]
sage: (x+x^5).derivative(x)
1

What is the integral of 1 w.r.t x in this case?
"

Could you be clearer? As I told, I'm not familiar with rings. I don't
even know the meaning of the argument of GF (I took the number 5 from
an example I see in sage-support group, I think). Do you think that QQ
[] could fit in this case? Moreover, what's the difference between QQ
and QQbar?

About this:
"
You can get this information from the factorization you compute above.
"
It seems that factorization over GF(p) is broken (see
http://groups.google.com/group/sage-devel/tree/browse_frm/thread/527e426dbfa04eda/2641f1509d1295db?rnum=1&_done=%2Fgroup%2Fsage-devel%2Fbrowse_frm%2Fthread%2F527e426dbfa04eda%3F#doc_2641f1509d1295db
), does it work for QQ? I don't seriously need factor, as long as I
need the root finding next.

Now let's go to Carl's help...

> Taking a quick look at that page, it looks like they want the exact
> roots in CC of a polynomial with algebraic coefficients.  In Sage, we
> can get this with QQbar:
>
> sage: K.<x> = QQbar[]
> sage: (x^5-x-1).roots(ring=QQbar)
>

First problem with QQbar: it seems that resultant() doesn't like it,
because it is not able to convert it to a Singular ring (this is the
error, I'm not attaching all the output, tell me if you need it)

TypeError: no conversion of this ring to a Singular ring defined

On the contrary, QQ[] seems to work fine with resultant (but it
doesn't have roots() )

Moreover, it seems that QQbar roots() is not working for multivariate
polynomials ring... is it true or am I just missing something else? In
that case, is possible to let it work in multivariate polynomials? As
you can imagine, I would like to think about this as a method of
solving integrals, so it is very likely to have a symbolic expression
with more than just a single symbolic variable.

> [(1.167303978261419?, 1),
>  (-0.7648844336005847? - 0.3524715460317263?*I, 1),
>  (-0.7648844336005847? + 0.3524715460317263?*I, 1),
>  (0.1812324444698754? - 1.083954101317711?*I, 1),
>  (0.1812324444698754? + 1.083954101317711?*I, 1)]
>
> (AFAIK, maxima can't do this; I don't think maxima can handle general
> algebraic numbers.  Since Sage's solve() is implemented using maxima,
> solve() won't work for this problem.)
>
> Carl

Apart from this, is there another way to solve an equation (with more
than a single symbolic variable) obtaining exact roots? It seems that
maxima would do the work (with algebraic numbers...), is it possible
that it is the only symbolic equation solver within SAGE? What about
SymPy or anything else?

Finally, I still would like to know which is the best way to translate
the output of a calculation with polynomial rings into a symbolic
expression, that can be carried on with maxima or pynac. Can you help
me?

Thank you very much

Regards

Maurizio

Carl Witty

unread,
Apr 18, 2009, 8:27:27 PM4/18/09
to sage-...@googlegroups.com
On Sat, Apr 18, 2009 at 2:46 PM, Maurizio <maurizio...@gmail.com> wrote:
> Could you be clearer? As I told, I'm not familiar with rings. I don't
> even know the meaning of the argument of GF (I took the number 5 from
> an example I see in sage-support group, I think). Do you think that QQ
> [] could fit in this case? Moreover, what's the difference between QQ
> and QQbar?

GF(5) means the Galois field of characteristic 5, a.k.a. the integers
modulo 5. (So in GF(5), you have 2*3 = 1, 3+4 = 2, etc.) It's
probably quite irrelevant for computing integrals.

QQ is the rational numbers (fractions). QQbar is the algebraic
closure of QQ; this means it includes every complex number which is
the root of a polynomial with rational coefficients. So it includes
things like sqrt(2) (which is a root of x^2-2), and sqrt(-1) (a root
of x^2+1), as well as more exotic numbers like the roots of x^5-x-1,
which can't be expressed using radicals (roots). (QQbar does not
include all complex numbers, though; for instance, it does not include
pi or e, which are transcendental rather than algebraic.)

> Now let's go to Carl's help...
>
>> Taking a quick look at that page, it looks like they want the exact
>> roots in CC of a polynomial with algebraic coefficients.  In Sage, we
>> can get this with QQbar:
>>
>> sage: K.<x> = QQbar[]
>> sage: (x^5-x-1).roots(ring=QQbar)
>>
>
> First problem with QQbar: it seems that resultant() doesn't like it,
> because it is not able to convert it to a Singular ring (this is the
> error, I'm not attaching all the output, tell me if you need it)
>
> TypeError: no conversion of this ring to a Singular ring defined

Looks like this hasn't been implemented yet.

> On the contrary, QQ[] seems to work fine with resultant (but it
> doesn't have roots() )

Univariate polynomials over QQ definitely have roots(); were you using
a multivariate polynomial ring?

> Moreover, it seems that QQbar roots() is not working for multivariate
> polynomials ring... is it true or am I just missing something else? In
> that case, is possible to let it work in multivariate polynomials? As
> you can imagine, I would like to think about this as a method of
> solving integrals, so it is very likely to have a symbolic expression
> with more than just a single symbolic variable.
>

> Apart from this, is there another way to solve an equation (with more
> than a single symbolic variable) obtaining exact roots? It seems that
> maxima would do the work (with algebraic numbers...), is it possible
> that it is the only symbolic equation solver within SAGE? What about
> SymPy or anything else?

What does this even mean? .roots() gives a list of all the solutions
of a univariate polynomial equation. But a multivariate polynomial
equation will usually have an infinite number of solutions; for
instance, x^2+y^2-1=0 has an infinite number of solutions over the
rationals (or the reals, or the algebraic reals, etc.)

If you have a system of multivariate polynomial equations, then the
system might have only finitely many solutions.

> Finally, I still would like to know which is the best way to translate
> the output of a calculation with polynomial rings into a symbolic
> expression, that can be carried on with maxima or pynac. Can you help
> me?

If p is a polynomial, then SR(p) is a symbolic expression.

Carl

Maurizio

unread,
Apr 18, 2009, 9:04:18 PM4/18/09
to sage-devel
Thanks for the answer.

As the time goes, I get more understanding of the complexity of the
problem (much more than I expected at first).

On 19 Apr, 02:27, Carl Witty <carl.wi...@gmail.com> wrote:
> On Sat, Apr 18, 2009 at 2:46 PM, Maurizio <maurizio.gran...@gmail.com> wrote:
> > Could you be clearer? As I told, I'm not familiar with rings. I don't
> > even know the meaning of the argument of GF (I took the number 5 from
> > an example I see in sage-support group, I think). Do you think that QQ
> > [] could fit in this case? Moreover, what's the difference between QQ
> > and QQbar?
>
> GF(5) means the Galois field of characteristic 5, a.k.a. the integers
> modulo 5.  (So in GF(5), you have 2*3 = 1, 3+4 = 2, etc.)  It's
> probably quite irrelevant for computing integrals.
>

Thank you, I should probably have looked for it on the internet rather
than annoying you :)

> QQ is the rational numbers (fractions).  QQbar is the algebraic
> closure of QQ; this means it includes every complex number which is
> the root of a polynomial with rational coefficients.  So it includes
> things like sqrt(2) (which is a root of x^2-2), and sqrt(-1) (a root
> of x^2+1), as well as more exotic numbers like the roots of x^5-x-1,
> which can't be expressed using radicals (roots).  (QQbar does not
> include all complex numbers, though; for instance, it does not include
> pi or e, which are transcendental rather than algebraic.)
>

I understand. This means (as I was kind of suspecting) that QQbar
could have been a good candidate to try to work with (I am just trying
to solve examples 1.4 and 1.5 of that web page). Unfortunately, QQbar
doesn't provide resultant() yet, so it should jump in after the
resultant has been calculated over QQ[]. I think that is not really an
elegant solution (and probably very suboptimal...).

> > Now let's go to Carl's help...
>
> >> Taking a quick look at that page, it looks like they want the exact
> >> roots in CC of a polynomial with algebraic coefficients.  In Sage, we
> >> can get this with QQbar:
>
> >> sage: K.<x> = QQbar[]
> >> sage: (x^5-x-1).roots(ring=QQbar)
>
> > First problem with QQbar: it seems that resultant() doesn't like it,
> > because it is not able to convert it to a Singular ring (this is the
> > error, I'm not attaching all the output, tell me if you need it)
>
> > TypeError: no conversion of this ring to a Singular ring defined
>
> Looks like this hasn't been implemented yet.

Unfortunately...

>
> > On the contrary, QQ[] seems to work fine with resultant (but it
> > doesn't have roots() )
>
> Univariate polynomials over QQ definitely have roots(); were you using
> a multivariate polynomial ring?
>

Actually, I was using a multivariate polynomial ring to let SAGE
understand the expression:

tores = A - z*B.derivative(x)

(by the way...wasn't there a trac opened to let derivative syntax
unify to diff or something like that? it seems that singular doesn't
like B.diff(x) yet... but that's very off-topic)

if I don't define a ring containing x and z, how do I inject z into
the workspace? certainly not as a symbolic variable, so what else? I
know that at the very end this is a univariate polynomial ring root
finding problem, I'm just wondering how to correctly manage it.

> > Moreover, it seems that QQbar roots() is not working for multivariate
> > polynomials ring... is it true or am I just missing something else? In
> > that case, is possible to let it work in multivariate polynomials? As
> > you can imagine, I would like to think about this as a method of
> > solving integrals, so it is very likely to have a symbolic expression
> > with more than just a single symbolic variable.
>
> > Apart from this, is there another way to solve an equation (with more
> > than a single symbolic variable) obtaining exact roots? It seems that
> > maxima would do the work (with algebraic numbers...), is it possible
> > that it is the only symbolic equation solver within SAGE? What about
> > SymPy or anything else?
>
> What does this even mean?  .roots() gives a list of all the solutions
> of a univariate polynomial equation.  But a multivariate polynomial
> equation will usually have an infinite number of solutions; for
> instance, x^2+y^2-1=0 has an infinite number of solutions over the
> rationals (or the reals, or the algebraic reals, etc.)
>
> If you have a system of multivariate polynomial equations, then the
> system might have only finitely many solutions.
>

I understand, but let's consider the problem I'm coming from:
indefinite integral. I'm not using polynomials because I'm well aware
of the implications, rather because I found resultant() implemented
there :)
What I intended with my sentences was: imagine a polynomial like: x^2
+ a^2. The solutions of x^2 - a^2 = 0 with respect to x over any field
wouldn't always coincide with +/- a? That is the kind of reasoning I
would make to solve the indefinite integral, like when I work with
symbolic variables that have a precise meaning (i.e. they are not
really variables in the problem, it's just that I want to keep that
degree of freedom with respect to the solution, especially if I want
to see the effect of that value over the solution) into my problem.

> > Finally, I still would like to know which is the best way to translate
> > the output of a calculation with polynomial rings into a symbolic
> > expression, that can be carried on with maxima or pynac. Can you help
> > me?
>
> If p is a polynomial, then SR(p) is a symbolic expression.
>

Funny, the next minutes I spent, I recognized how to do that in this
same way :) By the way, is there a SR equivalent for the new pynac
symbolic? Something like NSR (new symbolic ring)? I can see that pynac
already has gcd implemented, although I can already (unfortunately)
see:
RuntimeError: gcd: arguments must be polynomials over the rationals

when I try to do a gcd over quantities containing square roots and
complex numbers... hopefully pynac gcd could be ported to work with
QQbar :) Jokes apart, I have seen a very long thread (which took place
around april 2008) about multivariate gcd, and it seems a VERY big
issue, in the sense that it is not trivial to get an algorithm for gcd
robust, fast, and working over fields like QQbar (and possibly over
functionals)...

Basically, my aim was to see whether we did already have the tools to
implement the most basic indefinite integrals of rational functions,
but it seems that I already need: resultant(), solve(), gcd() and
these three are not that trivial to get... (probably I wasn't paying
attention to how important this basic functions are, when I took them
for granted, working with MatLab [Maple]...)

> Carl

Thanks and keep up the good work!

Maurizio

Carl Witty

unread,
Apr 18, 2009, 9:54:27 PM4/18/09
to sage-...@googlegroups.com
On Sat, Apr 18, 2009 at 6:04 PM, Maurizio <maurizio...@gmail.com> wrote:
>> QQ is the rational numbers (fractions).  QQbar is the algebraic
>> closure of QQ; this means it includes every complex number which is
>> the root of a polynomial with rational coefficients.  So it includes
>> things like sqrt(2) (which is a root of x^2-2), and sqrt(-1) (a root
>> of x^2+1), as well as more exotic numbers like the roots of x^5-x-1,
>> which can't be expressed using radicals (roots).  (QQbar does not
>> include all complex numbers, though; for instance, it does not include
>> pi or e, which are transcendental rather than algebraic.)
>>
>
> I understand. This means (as I was kind of suspecting) that QQbar
> could have been a good candidate to try to work with (I am just trying
> to solve examples 1.4 and 1.5 of that web page). Unfortunately, QQbar
> doesn't provide resultant() yet, so it should jump in after the
> resultant has been calculated over QQ[]. I think that is not really an
> elegant solution (and probably very suboptimal...).

It's not elegant, but it's probably much more efficient... QQbar is
vastly less efficient than QQ, so it's definitely a good idea to stick
with QQ as long as possible before jumping to QQbar. You can use
.roots(ring=QQbar) on a polynomial over QQ.

> if I don't define a ring containing x and z, how do I inject z into
> the workspace? certainly not as a symbolic variable, so what else? I
> know that at the very end this is a univariate polynomial ring root
> finding problem, I'm just wondering how to correctly manage it.

If you have a polynomial in a multivariate ring, but your actual
polynomial actually contains only one variable, you can get the
corresponding univariate polynomial with .univariate_polynomial().

> I understand, but let's consider the problem I'm coming from:
> indefinite integral. I'm not using polynomials because I'm well aware
> of the implications, rather because I found resultant() implemented
> there :)
> What I intended with my sentences was: imagine a polynomial like: x^2
> + a^2. The solutions of x^2 - a^2 = 0 with respect to x over any field
> wouldn't always coincide with +/- a? That is the kind of reasoning I
> would make to solve the indefinite integral, like when I work with
> symbolic variables that have a precise meaning (i.e. they are not
> really variables in the problem, it's just that I want to keep that
> degree of freedom with respect to the solution, especially if I want
> to see the effect of that value over the solution) into my problem.

I think it's going to be almost impossible to carry through the
algorithm with parameters... it's definitely not where you should be
starting :)

Consider x^5-x-a. In general I don't know of a better description of
the roots of this polynomial as a function of a than "the roots of
x^5-x-a as a function of a" :) (And if you want the distinct roots,
that's even worse. The number of distinct roots will depend on the
values of the parameters. Fortunately it probably doesn't matter
whether the roots are distinct, if you carefully adapt the
algorithm...)

>> > Finally, I still would like to know which is the best way to translate
>> > the output of a calculation with polynomial rings into a symbolic
>> > expression, that can be carried on with maxima or pynac. Can you help
>> > me?
>>
>> If p is a polynomial, then SR(p) is a symbolic expression.
>>
>
> Funny, the next minutes I spent, I recognized how to do that in this
> same way :) By the way, is there a SR equivalent for the new pynac
> symbolic? Something like NSR (new symbolic ring)? I can see that pynac
> already has gcd implemented, although I can already (unfortunately)
> see:
> RuntimeError: gcd: arguments must be polynomials over the rationals

There is an NSR, but it's not available on the command line by
default; you can get it with:

from sage.symbolic.ring import NSR

Carl

William Stein

unread,
Apr 18, 2009, 9:55:01 PM4/18/09
to sage-...@googlegroups.com
I think Singular doesn't have QQbar as a base ring, so I don't think
this is likely to be implemented for a while.

>
>> On the contrary, QQ[] seems to work fine with resultant (but it
>> doesn't have roots() )
>
> Univariate polynomials over QQ definitely have roots(); were you using
> a multivariate polynomial ring?
>
>> Moreover, it seems that QQbar roots() is not working for multivariate
>> polynomials ring... is it true or am I just missing something else? In
>> that case, is possible to let it work in multivariate polynomials? As
>> you can imagine, I would like to think about this as a method of
>> solving integrals, so it is very likely to have a symbolic expression
>> with more than just a single symbolic variable.
>>
>> Apart from this, is there another way to solve an equation (with more
>> than a single symbolic variable) obtaining exact roots? It seems that
>> maxima would do the work (with algebraic numbers...), is it possible
>> that it is the only symbolic equation solver within SAGE? What about
>> SymPy or anything else?
>
> What does this even mean?  .roots() gives a list of all the solutions
> of a univariate polynomial equation.  But a multivariate polynomial
> equation will usually have an infinite number of solutions; for
> instance, x^2+y^2-1=0 has an infinite number of solutions over the
> rationals (or the reals, or the algebraic reals, etc.)
>
> If you have a system of multivariate polynomial equations, then the
> system might have only finitely many solutions.

The "variety" command is relevant here.



>
>> Finally, I still would like to know which is the best way to translate
>> the output of a calculation with polynomial rings into a symbolic
>> expression, that can be carried on with maxima or pynac. Can you help
>> me?
>
> If p is a polynomial, then SR(p) is a symbolic expression.
>
> Carl
>
> >
>



--
William Stein
Associate Professor of Mathematics
University of Washington
http://wstein.org

Carl Witty

unread,
Apr 18, 2009, 10:03:48 PM4/18/09
to sage-...@googlegroups.com
On Sat, Apr 18, 2009 at 6:55 PM, William Stein <wst...@gmail.com> wrote:
>>> First problem with QQbar: it seems that resultant() doesn't like it,
>>> because it is not able to convert it to a Singular ring (this is the
>>> error, I'm not attaching all the output, tell me if you need it)
>>>
>>> TypeError: no conversion of this ring to a Singular ring defined
>>
>> Looks like this hasn't been implemented yet.
>
> I think Singular doesn't have QQbar as a base ring, so I don't think
> this is likely to be implemented for a while.

Well, we just need a resultant algorithm that doesn't go through
Singular. I'm planning to write such a thing as part of my
cylindrical algebraic decomposition implementation sometime in the
next few months.

Carl

Maurizio

unread,
Apr 19, 2009, 10:44:51 AM4/19/09
to sage-devel
Hi all

> Well, we just need a resultant algorithm that doesn't go through
> Singular.  I'm planning to write such a thing as part of my
> cylindrical algebraic decomposition implementation sometime in the
> next few months.
>
> Carl

yes, I agree with that.

William, unfortunately I can't understand what the function "variety"
will help for. Basically, I just wanted to see if I could implement
the most basic algorithm for rational functions integration, and I
followed the link in the Pynac wiki. Probably, my mistake was to
introduce polynomials in this path, but that's the only place where I
found the "resultant()" function. Do you have any alternative
suggestion?

Carl, I took advantage of your suggestion, even though I assume I
can't still go through the whole process with the current gcd
capabilities in Pynac. But before than that, I'd like to point out
something strange I did notice, and maybe also Burcin can help with
that:

reset()
# P.<x,y,z> = PolynomialRing(QQ)
# P.<x,y,z> = GF(5)[]
P.<x,z> = QQ[]

A = 1
B = x^3 + x

tores = A - z*diff(B,x)
res = tores.resultant(B,x); factor(res)
res1 = res.univariate_polynomial()
sol1 = res1.roots(ring = QQbar)

with this code, I get the roots over QQbar, which is useful. Then I'd
like to move to the symbolic field and I do this:

var('x, zs', ns = 1)
from sage.symbolic.ring import NSR
As = NSR(A)
Bs = NSR(B)
Bs
x^3 + x
Bs.diff(x)
0

So, the derivative is not working. Which is the cause? It seems that
the "x" in Bs is not the "x" I declared, so the derivative gets 0 as a
result. Which is the reason?

Assuming to go on (manually for the moment), I do:
c1 = QQbar(sol1[0][0])
v1a = A - c1*(3*x^2 + 1)
Bs.gcd(v1a)

Traceback (click to the left for traceback)
...
RuntimeError: gcd: arguments must be polynomials over the rationals

Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "/home/notebook/sage_notebook/worksheets/admin/4/code/135.py",
line 9, in <module>
Bs.gcd(v1a)
File "/usr/local/sage/local/lib/python2.5/site-packages/
zope.interface-3.3.0-py2.5-linux-i686.egg/", line 1, in <module>

File "expression.pyx", line 1624, in
sage.symbolic.expression.Expression.gcd (sage/symbolic/expression.cpp:
8608)
RuntimeError: gcd: arguments must be polynomials over the rationals

My point is: even though I could get the roots in QQbar (which are
exact), it seems that Pynac is not happy to work with QQbar
quantities, the only supported seems to be QQ pure rationals.

Moreover, I don't see this being the right way to do this, because
(for this particular problem: integration) I don't like having the
numerical representation of things like sqrt(5), even if the result is
still correct, so that

temp = QQbar(sqrt(5)); temp
2.236067977499790?

temp^2
5.000000000000000?

So, please tell me. Which should be the right way to try to approach
this indefinite integration problem? You can see that I'm not that
good in deep mathematical theory, but approaching the simplest problem
(that could be different from this one I'm looking at right now) is
fun :)

Regards

Maurizio
Message has been deleted

Maurizio

unread,
Apr 19, 2009, 11:38:09 AM4/19/09
to sage-devel
Challenging? Just because about the problem of integration in finite
terms Hardy in 1916 stated that “there is reason to suppose that no
such method can be given” ? :)

I want to add to this discussion that I found a lot of useful
information in this thread from SymPy list:
http://groups.google.com/group/sympy/browse_frm/thread/47259e49ad1cfd13/7a98521ffb13e311?lnk=gst&q=risch#7a98521ffb13e311

I'm wondering whether it could be a good idea to add to the Pynac wiki
a link to this article from Manuel Bronstein:
http://www-sop.inria.fr/cafe/Manuel.Bronstein/publications/issac98.pdf

This seems a good reading!

Moreover, in the same thread, I read that a lot of work has been
accomplished in SymPy. I know I can sometimes sound like a broken
record (always repeating the same crap...), but I more and more see
like a great advance in taking advantage of SymPy's python code. I am
willing to go and read it, and hopefully I can catch at least a single
bit from it :P I want to point out that this (in my opinion) is just
to give pynac a fastest bootstrap, so that we can start working on
improving something already there.

Unfortunately, porting some of SymPy's code sounds more like a
Computer Science task, which is not really my field.

If I look at my first aim (Laplace and Fourier transforms...) I can
see a long way to cover, hopefully in the shortest time frame!

Regards

Maurizio

On 19 Apr, 17:10, Martin Musatov <marty.musa...@gmail.com> wrote:
> Aha! Quite the challenge is it not?

Carl Witty

unread,
Apr 19, 2009, 1:34:21 PM4/19/09
to sage-...@googlegroups.com
On Sun, Apr 19, 2009 at 7:44 AM, Maurizio <maurizio...@gmail.com> wrote:
> Carl, I took advantage of your suggestion, even though I assume I
> can't still go through the whole process with the current gcd
> capabilities in Pynac. But before than that, I'd like to point out
> something strange I did notice, and maybe also Burcin can help with
> that:
>
> reset()
> P.<x,z> = QQ[]

>
> B = x^3 + x
>
> var('x, zs', ns = 1)
> from sage.symbolic.ring import NSR
> Bs = NSR(B)
> Bs
>     x^3 + x
> Bs.diff(x)
>     0
>
> So, the derivative is not working. Which is the cause? It seems that
> the "x" in Bs is not the "x" I declared, so the derivative gets 0 as a
> result. Which is the reason?

Looks like a bug to me.

Burcin, any comments?

> Moreover, I don't see this being the right way to do this, because
> (for this particular problem: integration) I don't like having the
> numerical representation of things like sqrt(5), even if the result is
> still correct, so that
>
> temp = QQbar(sqrt(5)); temp
> 2.236067977499790?
>
> temp^2
> 5.000000000000000?

Note that it's only the representation that's numerical; internally
these are actually exact numbers. (For instance, if you do "temp^2 ==
5" and Sage prints "True", then that means that temp really is exactly
the square root of 5.) I use a numeric representation because not all
algebraic numbers can be represented with radicals, and solving by
radicals is much more difficult than the symbolic-numeric algorithms I
used in QQbar.

> So, please tell me. Which should be the right way to try to approach
> this indefinite integration problem? You can see that I'm not that
> good in deep mathematical theory, but approaching the simplest problem
> (that could be different from this one I'm looking at right now) is
> fun :)

Unfortunately, I think the right way is a lot of work.

1) Teach yourself enough math that you can mostly understand the
algorithm in the paper. (Not just mechanically follow the steps, but
understand what each step does, and have some idea why the algorithm
as a whole works.)

For some parts of this task, wikipedia and mathworld.wolfram.com will
be useful resources (planetmath.org is also interersting); but you'll
probably also need to read some books, etc.

2) Implement each missing sub-algorithm in Sage.

3) Implement the entire algorithm in Sage.

This sounds glib, but it really is possible... I'm implementing
cylindrical algebraic decomposition in Sage following this program.
I've been working on it in my spare time for several years; I probably
spent a couple of years on step 1, then I started step 2 a couple of
years ago, and now I'm working on step 2 and step 3.

The problem is that the Risch algorithm just is not simple.

A good integration algorithm probably has a couple of phases -- first
you try some heuristics, and pattern-match against a database of known
integrals; if those fail, then you bring out the big guns and apply
the decision procedure. If you're interested in working on
integration, maybe you'd rather work on the first phase? That
probably requires more computer science and less math than the second
phase.

Carl

William Stein

unread,
Apr 19, 2009, 4:14:08 PM4/19/09
to sage-...@googlegroups.com
Indeed. I finally read up on it a bit, and it involves some ideas
from algebra/differential field theory that are familiar to me because
of my background in number theory and arithmetic geometry. The basic
idea is to generalize the partial fraction decomposition algorithm as
much as possible, first to algebraic extensions of QQ(t), then to
differential and log extensions.

Wikipedia also has a few interesting remarks, e.g., that the Risch
algorithm isn't an algorithm, because it depends on being able to
check equality of general elementary functions, which is evidently an
open problem in general (so in practice you just fake it by evaluating
numerically at lots of points to decide if something is probably equal
to something else). It's also evidently not implemented anywhere,
e.g., a nice example on the Wikipedia page, is that if you let

f = (x^2 + 2*x + 1 +
(3*x+1)*sqrt(x+log(x)))/(x*sqrt(x+log(x))*(x+sqrt(x+log(x))))

then it has the antiderivative

g = 2*(sqrt(x+log(x)) + log(x+sqrt(x+log(x))))

since

sage: h = g.derivative() - f
sage: h.full_simplify()
0

However, Sage, Maple, and Mathematica, all simply give "integral(f)"
back when asked to integrate f. (I just checked this with the latest
versions.)

Fortunately, if one writes a program, using any number of heuristics,
to compute an antiderivative g, one can compute the derivative of g
and test whether it is *probably* equal to f by evaluating at lots of
random points. I wonder if some programs do some of the integration
by working numerically, getting a g with floating point coefficients
to high precision, then using LLL to recognize them as exact algebraic
numbers, then later deduce with high probability (and sometimes with
certainty) that g is correct?

Anyway, I'm personally definitely not at all scared of implementing
some heuristic variant of the Risch procedure, and I don't think other
people should be either. I've certainly encountered more complicated
"algorithms"...


>
> A good integration algorithm probably has a couple of phases -- first
> you try some heuristics, and pattern-match against a database of known
> integrals; if those fail, then you bring out the big guns and apply
> the decision procedure.  If you're interested in working on
> integration, maybe you'd rather work on the first phase?  That
> probably requires more computer science and less math than the second
> phase.
>
> Carl
>
> >
>



Fredrik Johansson

unread,
Apr 19, 2009, 4:23:14 PM4/19/09
to sage-...@googlegroups.com
On Sun, Apr 19, 2009 at 10:14 PM, William Stein <wst...@gmail.com> wrote:
> Wikipedia also has a few interesting remarks, e.g., that the Risch
> algorithm isn't an algorithm, because it depends on being able to
> check equality of general elementary functions, which is evidently an
> open problem in general (so in practice you just fake it by evaluating
> numerically at lots of points to decide if something is probably equal
> to something else).   It's also evidently not implemented anywhere,
> e.g., a nice example on the Wikipedia page, is that if you let
>
> f = (x^2 + 2*x + 1 +
> (3*x+1)*sqrt(x+log(x)))/(x*sqrt(x+log(x))*(x+sqrt(x+log(x))))
>
> then it has the antiderivative
>
> g = 2*(sqrt(x+log(x)) + log(x+sqrt(x+log(x))))
>
> since
>
> sage: h = g.derivative() - f
> sage: h.full_simplify()
> 0
>
> However, Sage, Maple, and Mathematica, all simply give "integral(f)"
> back when asked to integrate f.  (I just checked this with the latest
> versions.)

Curiously though, SymPy knows this particular integral.

>>> from sympy import *
>>> x=Symbol('x')
>>> f=(x**2+2*x+1+(3*x+1)*sqrt(x+log(x)))/(x*sqrt(x+log(x))*(x+sqrt(x+log(x))))
>>> integrate(f,x)
2*log(x + (x + log(x))**(1/2)) + 2*(x + log(x))**(1/2)

Fredrik

root

unread,
Apr 19, 2009, 6:33:12 PM4/19/09
to sage-...@googlegroups.com, sage-...@googlegroups.com

A much shorter example is:

integrate(sqrt(x+log(x)),x)

to which Axiom replies:

integrate: implementation incomplete (constant residues)

Tim Daly

root

unread,
Apr 19, 2009, 6:41:42 PM4/19/09
to sage-...@googlegroups.com, sage-...@googlegroups.com

A much shorter example is:

integrate(sqrt(x+log(x)),x)

to which Axiom replies:

integrate: implementation incomplete (constant residues)

For further information, Bronstein deals with constant residues
(written after his Axiom work) in his book [1] after page 147.

Tim Daly

[1] Bronstein, Manuel "Symbolic Integration I", Second Edition
Springer ISBN 3-540-21493-3

William Stein

unread,
Apr 19, 2009, 6:06:59 PM4/19/09
to sage-...@googlegroups.com

What is f(x) = sqrt(x+log(x)) supposed to be an example of? Does f
has an antiderivative that can be expressed in terms of elementary
functions? If so, what is it?

-- William

Ondrej Certik

unread,
Apr 19, 2009, 6:25:59 PM4/19/09
to sage-...@googlegroups.com

I was just about to point it out. It even works in Sage:

In [2]: f = (x^2 + 2*x + 1 +
...: (3*x+1)*sqrt(x+log(x)))/(x*sqrt(x+log(x))*(x+sqrt(x+log(x))))
In [5]: import sympy
In [6]: sympy.integrate(f, x)
Out[6]: 2*log(x + (x + log(x))**(1/2)) + 2*(x + log(x))**(1/2)


So we have a good start to implement the Risch algorithm in sympy already.

Ondrej

William Stein

unread,
Apr 19, 2009, 6:28:54 PM4/19/09
to sage-...@googlegroups.com

Ondrej,

How does sympy compute the integral of

(x^2 + 2*x + 1 + (3*x+1)*sqrt(x+log(x)))/(x*sqrt(x+log(x))*(x+sqrt(x+log(x))))?

William

root

unread,
Apr 19, 2009, 7:20:00 PM4/19/09
to sage-...@googlegroups.com, sage-...@googlegroups.com
> So we have a good start to implement the Risch algorithm in sympy already.

Ondrej, what result do you get for:

integrate(sqrt(x+log(x)),x)

Tim

root

unread,
Apr 19, 2009, 7:17:28 PM4/19/09
to sage-...@googlegroups.com, sage-...@googlegroups.com
> > A much shorter example is:
> >
> > integrate(sqrt(x+log(x)),x)
> >
> > to which Axiom replies:
> >
> > integrate: implementation incomplete (constant residues)
> >
>
> What is f(x) = sqrt(x+log(x)) supposed to be an example of? Does f
> has an antiderivative that can be expressed in terms of elementary
> functions? If so, what is it?

Axiom has one of the most complete Risch implementations. Many of
the original authors of the theory worked on the implementation.

Axiom can almost do the Risch integral of the formula you posted.
It fails with:

integrate: implementation incomplete (constant residues)

The failing subexpression sqrt(x+log(x)) falls into the unimplemented
case. Replacing this subexpression by 'a' succeeds.

Most systems do pattern matching and heuristics. Fateman had a server
available that would "integrate" using these techniques.

Tim


Ondrej Certik

unread,
Apr 19, 2009, 6:57:35 PM4/19/09
to sage-...@googlegroups.com

In [1]: integrate(sqrt(x+log(x)),x)
Out[1]:

⎮ ⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽
⎮ ╲╱ x + log(x) dx


So we can't do it yet.

Ondrej

Ondrej Certik

unread,
Apr 20, 2009, 1:18:28 AM4/20/09
to sage-...@googlegroups.com
On Sun, Apr 19, 2009 at 3:28 PM, William Stein <wst...@gmail.com> wrote:

> Ondrej,
>
> How does sympy compute the integral of
>
> (x^2 + 2*x + 1 + (3*x+1)*sqrt(x+log(x)))/(x*sqrt(x+log(x))*(x+sqrt(x+log(x))))?

See this thread in our list:

http://groups.google.com/group/sympy/browse_thread/thread/9a61d681e6f967e8/

Ondrej

Maurizio

unread,
Apr 20, 2009, 8:31:59 AM4/20/09
to sage-devel
Kudos to SymPy!

I'm wondering why the python integration algorithms implemented there
aren't in the short term adopted by SAGE.
At least, they are already aware of their shortcomings (ie: cannot
compute the integral of log(x)/x ).

I'm sure SAGE people could give big contribute to those, send patches
upstream, and move forward if needed.

Regards and many thanks

Maurizio

On Apr 20, 12:57 am, Ondrej Certik <ond...@certik.cz> wrote:

Burcin Erocal

unread,
Apr 20, 2009, 3:42:01 PM4/20/09
to sage-...@googlegroups.com
On Sun, 19 Apr 2009 10:34:21 -0700
Carl Witty <carl....@gmail.com> wrote:

>
> On Sun, Apr 19, 2009 at 7:44 AM, Maurizio
> <maurizio...@gmail.com> wrote:
> > Carl, I took advantage of your suggestion, even though I assume I
> > can't still go through the whole process with the current gcd
> > capabilities in Pynac. But before than that, I'd like to point out
> > something strange I did notice, and maybe also Burcin can help with
> > that:
> >
> > reset()
> > P.<x,z> = QQ[]
> >
> > B = x^3 + x
> >
> > var('x, zs', ns = 1)
> > from sage.symbolic.ring import NSR
> > Bs = NSR(B)
> > Bs
> >     x^3 + x
> > Bs.diff(x)
> >     0
> >
> > So, the derivative is not working. Which is the cause? It seems that
> > the "x" in Bs is not the "x" I declared, so the derivative gets 0
> > as a result. Which is the reason?
>
> Looks like a bug to me.
>
> Burcin, any comments?

I agree that it's confusing, but it's not a bug.

The command

sage: Bs = NSR(B)

converts the polynomial B = x^3 + x in QQ[x] to a symbolic expression,
with one numeric coefficient, namely B.

sage: Bs.coeff(x)
0
sage: Bs.is_constant()
True
sage: Bs.pyobject()
x^3 + x
sage: Bs.pyobject().parent()
Multivariate Polynomial Ring in x, z over Rational Field

Since Bs is a constant w.r.t. the symbolic x, the derivative of Bs
w.r.t. the symbolic variable x is 0.


In order to convert the polynomial to a symbolic expression, I suggest
substituting the symbolic variables in the polynomial.

sage: B.subs(x=x)
x^3 + x
sage: t = B.subs(x=x)
sage: t.coeff(x)
1
sage: t.coeff(x^3)
1
sage: t.derivative(x)
3*x^2 + 1

<snip>


> A good integration algorithm probably has a couple of phases -- first
> you try some heuristics, and pattern-match against a database of known
> integrals; if those fail, then you bring out the big guns and apply
> the decision procedure. If you're interested in working on
> integration, maybe you'd rather work on the first phase? That
> probably requires more computer science and less math than the second
> phase.


Indeed, this was what I tried to say in those comments about
integration that got Maurizio started on this track.

There are descriptions of heuristics we can use in the maxima manual [1]
and the book by Geddes, Czapor and Labahn [2].

[1]
http://maxima.sourceforge.net/docs/manual/en/maxima_20.html#Item_003a-integrate

[2] http://books.google.at/books?id=9fOUwkkRxT4C&pg=PR2&hl=en#PPA473,M1

It would be great if someone went through these links, and made them
more explicit/precise. As I said before, implementing these should be
within reach with the pattern matching/rewriting capabilities of pynac.


Here is an example, pointed out to me by Clemens Raab, which maxima can
do and MMA can't:

sage: g=log(1-x)/(polylog(2,x)^2+log(1-x)*polylog(3,x))
sage: f = g.derivative(x)
sage: integral(f,x)
log(1 - x)/(log(1 - x)*polylog(3, x) + polylog(2, x)^2)


Cheers,
Burcin

Jason Grout

unread,
Apr 20, 2009, 3:53:52 PM4/20/09
to sage-...@googlegroups.com


So is this what happens anytime you do NSR(sage_object)? It converts it
to a constant (as far as symbolics is concerned) with a coefficient of
sage_object? Or are there any exceptions?

Jason

--
Jason Grout

Maurizio

unread,
Apr 20, 2009, 4:12:18 PM4/20/09
to sage-devel
Hi Burcin, thanks for replying!

> I agree that it's confusing, but it's not a bug.
>
> The command
>
> sage: Bs = NSR(B)
>
> converts the polynomial B = x^3 + x in QQ[x] to a symbolic expression,
> with one numeric coefficient, namely B.
>

Excuse me, but I don't understand the reason for this. When does this
come handful, if you can't directly translate that polynomial into a
symbolic expression?

> sage: Bs.coeff(x)
> 0
> sage: Bs.is_constant()
> True
> sage: Bs.pyobject()
> x^3 + x
> sage: Bs.pyobject().parent()
> Multivariate Polynomial Ring in x, z over Rational Field
>
> Since Bs is a constant w.r.t. the symbolic x, the derivative of Bs
> w.r.t. the symbolic variable x is 0.
>
> In order to convert the polynomial to a symbolic expression, I suggest
> substituting the symbolic variables in the polynomial.
>
> sage: B.subs(x=x)
> x^3 + x
> sage: t = B.subs(x=x)

Is really this the correct way? Are you satisfied with this? Even if
its user friendliness could be questionable, I find it rather
inelegant... what do you think? I would consider switching the normal
behavior, if possible.

> sage: t.coeff(x)
> 1
> sage: t.coeff(x^3)
> 1
> sage: t.derivative(x)
> 3*x^2 + 1
>
> <snip>
>
> > A good integration algorithm probably has a couple of phases -- first
> > you try some heuristics, and pattern-match against a database of known
> > integrals; if those fail, then you bring out the big guns and apply
> > the decision procedure.  If you're interested in working on
> > integration, maybe you'd rather work on the first phase?  That
> > probably requires more computer science and less math than the second
> > phase.
>
> Indeed, this was what I tried to say in those comments about
> integration that got Maurizio started on this track.
>
> There are descriptions of heuristics we can use in the maxima manual [1]
> and the book by Geddes, Czapor and Labahn [2].
>
> [1]http://maxima.sourceforge.net/docs/manual/en/maxima_20.html#Item_003a...
>
> [2]http://books.google.at/books?id=9fOUwkkRxT4C&pg=PR2&hl=en#PPA473,M1
>
> It would be great if someone went through these links, and made them
> more explicit/precise. As I said before, implementing these should be
> within reach with the pattern matching/rewriting capabilities of pynac.
>
> Here is an example, pointed out to me by Clemens Raab, which maxima can
> do and MMA can't:
>
> sage: g=log(1-x)/(polylog(2,x)^2+log(1-x)*polylog(3,x))
> sage: f = g.derivative(x)
> sage: integral(f,x)
> log(1 - x)/(log(1 - x)*polylog(3, x) + polylog(2, x)^2)
>
> Cheers,
> Burcin

I don't know what about those algorithms, but it seems to me that
SymPy already implements some good heuristics, which can solve
integrals that Mathematica can't. So can we take this as a starting
point? I see that it is certainly possible to implement everything
from the beginning, but a bootstrapped start seems better to me in the
short term, because this could provide the user the functionalities
needed. This could be useful to speed up the switch to the new
symbolic system (I think you wouldn't do that without a good
integration engine), and once there, you got a great exposure to bug
fixing by users.

Thanks

Maurizio

Robert Bradshaw

unread,
Apr 20, 2009, 4:27:37 PM4/20/09
to sage-...@googlegroups.com

I consider this bad behavior as well. Perhaps it should at least try
to call a _symbolic_ method.

- Robert


mabshoff

unread,
Apr 20, 2009, 7:12:30 PM4/20/09
to sage-devel


On Apr 20, 1:12 pm, Maurizio <maurizio.gran...@gmail.com> wrote:
> Hi Burcin, thanks for replying!

<SNIP>

> I don't know what about those algorithms, but it seems to me that
> SymPy already implements some good heuristics, which can solve
> integrals that Mathematica can't.

Well, there are many, many integrals that MMA can do that Sympy
cannot. So it doesn't matter that in this specific case SymPy solves
it.

> So can we take this as a starting
> point? I see that it is certainly possible to implement everything
> from the beginning, but a bootstrapped start seems better to me in the
> short term, because this could provide the user the functionalities
> needed. This could be useful to speed up the switch to the new
> symbolic system (I think you wouldn't do that without a good
> integration engine), and once there, you got a great exposure to bug
> fixing by users.

Well, talking about it won't make it happen :). People working on
Symbolics in Sage are well aware of Sympy and its capabilities, but it
is BSD (i.e. some people won't work on it) and abstract mathematical
capabilities aren't as well developed as in Sage and often quite a bit
slower. Pynac is building Symbolics from the ground up and it has
taken a while to get to where we are and another couple months more or
less won't make a difference at this point.

> Thanks
>
> Maurizio

Cheers,

Michael

Maurizio

unread,
Apr 21, 2009, 2:44:31 AM4/21/09
to sage-devel
Hi Michael,

Actually, I thought that this discussion (especially people much more
expert than me) has clarified the point that implementing integrals is
not really just matter of a couple of months... but I would be glad to
see this happen!

I know there are some license issues with SymPy (not really issues,
just differences probably) but I think there's no problem in taking
inspiration and some pieces of code from it, right?

I'm saying this, because I can see this new symbolic stuff exciting,
but without the right amount of interest on it, its development will
always be a little slow

Regards

Maurizio

Robert Bradshaw

unread,
Apr 21, 2009, 3:32:55 AM4/21/09
to sage-...@googlegroups.com
On Apr 20, 2009, at 11:44 PM, Maurizio wrote:

> Hi Michael,
>
> Actually, I thought that this discussion (especially people much more
> expert than me) has clarified the point that implementing integrals is
> not really just matter of a couple of months... but I would be glad to
> see this happen!

Implementing something that can handle a huge number of integrals is
a reasonable goal for a couple of months. Implementing something that
can handle everything that we know how to handle in theory...well,
that hasn't ever happened yet.

> I know there are some license issues with SymPy (not really issues,
> just differences probably) but I think there's no problem in taking
> inspiration and some pieces of code from it, right?

There is a problem with just lifting code--but we can and do ship
SymPy as part of Sage.

> I'm saying this, because I can see this new symbolic stuff exciting,
> but without the right amount of interest on it, its development will
> always be a little slow

Given that we ship SymPy, so anything it can handle we should be able
to handle. I imagine when you want to integrate something, it will
try to do it natively first, and that failing ask SymPy and/or maxima
before returning a failure. Over time we'll depend less and less on
the other two.

- Robert

mabshoff

unread,
Apr 21, 2009, 3:44:52 AM4/21/09
to sage-devel


On Apr 21, 12:32 am, Robert Bradshaw <rober...@math.washington.edu>
wrote:
> On Apr 20, 2009, at 11:44 PM, Maurizio wrote:
>
> > Hi Michael,
>
> > Actually, I thought that this discussion (especially people much more
> > expert than me) has clarified the point that implementing integrals is
> > not really just matter of a couple of months... but I would be glad to
> > see this happen!
>
> Implementing something that can handle a huge number of integrals is  
> a reasonable goal for a couple of months. Implementing something that  
> can handle everything that we know how to handle in theory...well,  
> that hasn't ever happened yet.

Yep, certainly implementation are certainly better than other
implementation here on average, but one could claim that neither Sympy
nor Maxima are at the head of the class. Maxima can do many things,
but often a little massaging is required by the user to get optimum
results comparable to MMA for example. And I believe plainly and
simply that this is not the role of the user of a CAS to be an expert
at term maniplulation. If I write intgerate($FOO) I expect to get a
correct answer without being required to transform the expression due
to knowledge one has about the underlying implementation.

> > I know there are some license issues with SymPy (not really issues,
> > just differences probably) but I think there's no problem in taking
> > inspiration and some pieces of code from it, right?
>
> There is a problem with just lifting code--but we can and do ship  
> SymPy as part of Sage.

Yes. One could take code from Sympy (the BSD license allows this) and
make GPL ed changes on top of it. The main issue is just that our
pattern matching engine via pynac will be more powerful, our
arithmetic is faster and the other abstract math bits in Sage are way
more powerful than what is in Sympy. And the goal of Sympy is to be
self contained and depending on certain operations in Sage is not an
option and not wanted, so contributing such code back to Sympy is not
an option . There is the goal for Sympy to optionally depend on Pynac
and use it when available, but no one is working on this, so this is
not something we ought to be waiting for.

Plainly and simply: Waiting on someone else to fix the problem for us
has not worked, so we will do it ourselves.

> > I'm saying this, because I can see this new symbolic stuff exciting,
> > but without the right amount of interest on it, its development will
> > always be a little slow
>
> Given that we ship SymPy, so anything it can handle we should be able  
> to handle. I imagine when you want to integrate something, it will  
> try to do it natively first, and that failing ask SymPy and/or maxima  
> before returning a failure. Over time we'll depend less and less on  
> the other two.

Yep, that is the only way in my opinion. We control the Sage library
and can coordinate releases with improvements in the Sage library -
Sympy has not released as frequently as it used to be and unless
someone steps up and helps Ondrej out I doubt that is going to change
any time soon.

And by the way: the latest Symbolics code that is not yet available
for review can already user either Sympy and/or Maxima to do
integration. That code ought to be in Sage 4.0.

> - Robert

Cheers,

Michael

Ondrej Certik

unread,
Apr 21, 2009, 4:06:56 PM4/21/09
to sage-...@googlegroups.com
On Mon, Apr 20, 2009 at 5:31 AM, Maurizio <maurizio...@gmail.com> wrote:
>
> Kudos to SymPy!
>
> I'm wondering why the python integration algorithms implemented there
> aren't in the short term adopted by SAGE.

They are --- you can use them from sympy inside Sage. It's my goal
that all sympy features are nicely integrated in Sage. I work on this
as time permits.

> At least, they are already aware of their shortcomings (ie: cannot
> compute the integral of log(x)/x ).
>
> I'm sure SAGE people could give big contribute to those, send patches
> upstream, and move forward if needed.

We got 5 summer of code students this year, so we'll move sympy a lot
forward this summer. The most important project:

http://socghop.appspot.com/student_project/show/google/gsoc2009/portland_state/t124024247737

will disentangle the assumptions, so that we should then be able to
easily use any other core, like pynac, or our own cython core. Besides
that when my school ends in less than a month, I plan to work hard on
SPD:

http://code.google.com/p/spdproject/

so that one can easily create Sage like projects with his own
packages. I expect this to make it easier for people to create
libraries in a way, so that they are easy to use with Sage, e.g. they
work with plotting, etc., but they could distribute their own version
of all in one distribution.

In any case, this will be an exciting summer.

Ondrej

William Stein

unread,
Apr 21, 2009, 4:10:47 PM4/21/09
to sage-...@googlegroups.com, Mike Hansen, Dan Shumow
On Tue, Apr 21, 2009 at 1:06 PM, Ondrej Certik <ond...@certik.cz> wrote:
>
> On Mon, Apr 20, 2009 at 5:31 AM, Maurizio <maurizio...@gmail.com> wrote:
>>
>> Kudos to SymPy!
>>
>> I'm wondering why the python integration algorithms implemented there
>> aren't in the short term adopted by SAGE.
>
> They are --- you can use them from sympy inside Sage. It's my goal
> that all sympy features are nicely integrated in Sage. I work on this
> as time permits.

Also, in the pynac-based symbolics that Mike Hansen has been polishing
up for full inclusion in Sage (to replace the maxima based symbolics),
one can just do

f.integrate(algorithm="sympy")

and sage will compute the integral using sympy. He's already
implemented this and I've seen it work well when I tried it out.

>> At least, they are already aware of their shortcomings (ie: cannot
>> compute the integral of log(x)/x ).
>>
>> I'm sure SAGE people could give big contribute to those, send patches
>> upstream, and move forward if needed.
>
> We got 5 summer of code students this year, so we'll move sympy a lot
> forward this summer. The most important project:
>
> http://socghop.appspot.com/student_project/show/google/gsoc2009/portland_state/t124024247737
>
> will disentangle the assumptions, so that we should then be able to
> easily use any other core, like pynac, or our own cython core. Besides
> that when my school ends in less than a month, I plan to work hard on
> SPD:
>
> http://code.google.com/p/spdproject/
>
> so that one can easily create Sage like projects with his own
> packages. I expect this to make it easier for people to create
> libraries in a way, so that they are easy to use with Sage, e.g. they
> work with plotting, etc., but they could distribute their own version
> of all in one distribution.
>
> In any case, this will be an exciting summer.

Cool. Are you going to keep Microsoft Windows in mind wrt SPD? I
think the work Dan Shumow (and others) have been
doing on the Sage windows port could, in a sense, be seen as a base
for a completely open SPD for Windows. See
http://windows.sagemath.org/


-- William

Ondrej Certik

unread,
Apr 21, 2009, 4:25:50 PM4/21/09
to sage-...@googlegroups.com, Mike Hansen, Dan Shumow
On Tue, Apr 21, 2009 at 1:10 PM, William Stein <wst...@gmail.com> wrote:
>
> On Tue, Apr 21, 2009 at 1:06 PM, Ondrej Certik <ond...@certik.cz> wrote:
>>
>> On Mon, Apr 20, 2009 at 5:31 AM, Maurizio <maurizio...@gmail.com> wrote:
>>>
>>> Kudos to SymPy!
>>>
>>> I'm wondering why the python integration algorithms implemented there
>>> aren't in the short term adopted by SAGE.
>>
>> They are --- you can use them from sympy inside Sage. It's my goal
>> that all sympy features are nicely integrated in Sage. I work on this
>> as time permits.
>
> Also, in the pynac-based symbolics that Mike Hansen has been polishing
> up for full inclusion in Sage (to replace the maxima based symbolics),
> one can just do
>
>   f.integrate(algorithm="sympy")
>
> and sage will compute the integral using sympy.   He's already
> implemented this and I've seen it work well when I tried it out.

Yes, that way Sage allows to call any of those libraries very easily
in Sage and at the same time uses just one library by default, that
currently is the best at integrating, I guess still Maxima.

>
>>> At least, they are already aware of their shortcomings (ie: cannot
>>> compute the integral of log(x)/x ).
>>>
>>> I'm sure SAGE people could give big contribute to those, send patches
>>> upstream, and move forward if needed.
>>
>> We got 5 summer of code students this year, so we'll move sympy a lot
>> forward this summer. The most important project:
>>
>> http://socghop.appspot.com/student_project/show/google/gsoc2009/portland_state/t124024247737
>>
>> will disentangle the assumptions, so that we should then be able to
>> easily use any other core, like pynac, or our own cython core. Besides
>> that when my school ends in less than a month, I plan to work hard on
>> SPD:
>>
>> http://code.google.com/p/spdproject/
>>
>> so that one can easily create Sage like projects with his own
>> packages. I expect this to make it easier for people to create
>> libraries in a way, so that they are easy to use with Sage, e.g. they
>> work with plotting, etc., but they could distribute their own version
>> of all in one distribution.
>>
>> In any case, this will be an exciting summer.
>
> Cool.  Are you going to keep Microsoft Windows in mind wrt SPD?    I

Yes, definitely, working on windows is one of the goals.

> think the work Dan Shumow (and others) have been
> doing on the Sage windows port could, in a sense, be seen as a base
> for a completely open SPD for Windows. See
>                 http://windows.sagemath.org/

Excellent, I just joined the list and will work closely with him.

Ondrej

Jason Grout

unread,
Apr 21, 2009, 4:58:00 PM4/21/09
to sage-...@googlegroups.com
Ondrej Certik wrote:
> On Tue, Apr 21, 2009 at 1:10 PM, William Stein <wst...@gmail.com> wrote:
>> On Tue, Apr 21, 2009 at 1:06 PM, Ondrej Certik <ond...@certik.cz> wrote:
>>> On Mon, Apr 20, 2009 at 5:31 AM, Maurizio <maurizio...@gmail.com> wrote:
>>>> Kudos to SymPy!
>>>>
>>>> I'm wondering why the python integration algorithms implemented there
>>>> aren't in the short term adopted by SAGE.
>>> They are --- you can use them from sympy inside Sage. It's my goal
>>> that all sympy features are nicely integrated in Sage. I work on this
>>> as time permits.
>> Also, in the pynac-based symbolics that Mike Hansen has been polishing
>> up for full inclusion in Sage (to replace the maxima based symbolics),
>> one can just do
>>
>> f.integrate(algorithm="sympy")
>>
>> and sage will compute the integral using sympy. He's already
>> implemented this and I've seen it work well when I tried it out.
>
> Yes, that way Sage allows to call any of those libraries very easily
> in Sage and at the same time uses just one library by default, that
> currently is the best at integrating, I guess still Maxima.

unless maxima doesn't produce an answer (i.e., returns your expression
unevaluated or throws an exception because of an assumption, etc.) Then
it could automatically try sympy.

-Jason

--
Jason Grout

Maurizio

unread,
Apr 21, 2009, 5:33:06 PM4/21/09
to sage-devel
thank you for clarifying this, I didn't know that sympy was already
pretty well working with the new symbolics

I hope this has at least given some information to the community as
well! I always learn a lot from these discussions

Regards

Maurizio

On 21 Apr, 22:58, Jason Grout <jason-s...@creativetrax.com> wrote:
> Ondrej Certik wrote:
> > On Tue, Apr 21, 2009 at 1:10 PM, William Stein <wst...@gmail.com> wrote:
> >> On Tue, Apr 21, 2009 at 1:06 PM, Ondrej Certik <ond...@certik.cz> wrote:

Ondrej Certik

unread,
Apr 21, 2009, 6:48:57 PM4/21/09
to sage-...@googlegroups.com
On Tue, Apr 21, 2009 at 2:33 PM, Maurizio <maurizio...@gmail.com> wrote:
>
> thank you for clarifying this, I didn't know that sympy was already
> pretty well working with the new symbolics
>
> I hope this has at least given some information to the community as
> well! I always learn a lot from these discussions

Excellent. Please report any bugs that you discover in the sympy<->
sage conversion. I believe there still may be some.

Thanks,
Ondrej

William Stein

unread,
Apr 22, 2009, 12:56:52 AM4/22/09
to sage-...@googlegroups.com

Please don't report any bugs you discover in the sympy<-->sage conversion
right now, since Mike Hansen literally wrote a bunch of it a few days ago,
and it exists only on his laptop.

Once Sage-4.0 appears with the new symbolics, then we'll definitely
want to hear about bugs.

William


--
William Stein
Associate Professor of Mathematics
University of Washington
http://wstein.org

Ondrej Certik

unread,
Apr 22, 2009, 12:08:51 PM4/22/09
to sage-...@googlegroups.com
On Tue, Apr 21, 2009 at 12:44 AM, mabshoff <mabs...@googlemail.com> wrote:
>
>
>
> On Apr 21, 12:32 am, Robert Bradshaw <rober...@math.washington.edu>
> wrote:
>> On Apr 20, 2009, at 11:44 PM, Maurizio wrote:
>>
>> > Hi Michael,
>>
>> > Actually, I thought that this discussion (especially people much more
>> > expert than me) has clarified the point that implementing integrals is
>> > not really just matter of a couple of months... but I would be glad to
>> > see this happen!
>>
>> Implementing something that can handle a huge number of integrals is
>> a reasonable goal for a couple of months. Implementing something that
>> can handle everything that we know how to handle in theory...well,
>> that hasn't ever happened yet.
>
> Yep, certainly implementation are certainly better than other
> implementation here on average, but one could claim that neither Sympy
> nor Maxima are at the head of the class. Maxima can do many things,
> but often a little massaging is required by the user to get optimum
> results comparable to MMA for example. And I believe plainly and
> simply that this is not the role of the user of a CAS to be an exper
> at term maniplulation. If I write intgerate($FOO) I expect to get a
> correct answer without being required to transform the expression due
> to knowledge one has about the underlying implementation.

I believe that too.

>
>> > I know there are some license issues with SymPy (not really issues,
>> > just differences probably) but I think there's no problem in taking
>> > inspiration and some pieces of code from it, right?
>>
>> There is a problem with just lifting code--but we can and do ship
>> SymPy as part of Sage.
>
> Yes. One could take code from Sympy (the BSD license allows this) and
> make GPL ed changes on top of it. The main issue is just that our
> pattern matching engine via pynac will be more powerful, our

Could you give some specific examples where the pynac pattern engine
is more powerful?

> arithmetic is faster and the other abstract math bits in Sage are way
> more powerful than what is in Sympy. And the goal of Sympy is to be
> self contained and depending on certain operations in Sage is not an
> option and not wanted, so contributing such code back to Sympy is not

We want sympy to run without Sage. But I have absolutely no problems
calling Sage for certain things if it's available. In fact, I do want
to call Sage if it's available.

> an option . There is the goal for Sympy to optionally depend on Pynac
> and use it when available, but no one is working on this, so this is
> not something we ought to be waiting for.
>
> Plainly and simply: Waiting on someone else to fix the problem for us
> has not worked, so we will do it ourselves.

For the record, I offered William and you that I will work on this, if
we manage to find funding for it over the summer. So just to be clear
that I am very interested in this, and not just talking. But
unfortunately I myself didn't find funding for it (and I haven't heard
from you either that you found some possibility), so I had to find an
internship somewhere else related more to my research (finite
elements).

We managed to get one gsoc project that does the assumptions right, so
it may happen anyways over the summer, in fact I very much hope so.

>> > I'm saying this, because I can see this new symbolic stuff exciting,
>> > but without the right amount of interest on it, its development will
>> > always be a little slow
>>
>> Given that we ship SymPy, so anything it can handle we should be able
>> to handle. I imagine when you want to integrate something, it will
>> try to do it natively first, and that failing ask SymPy and/or maxima
>> before returning a failure. Over time we'll depend less and less on
>> the other two.
>
> Yep, that is the only way in my opinion. We control the Sage library
> and can coordinate releases with improvements in the Sage library -
> Sympy has not released as frequently as it used to be and unless
> someone steps up and helps Ondrej out I doubt that is going to change
> any time soon.

We'll see. We managed to get 5 gsoc student this summer, so we are not
at all dead, if this is what you mean. :)

Besides, it will take almost a year to Sage too to release the new
symbolics (started August 2008), so I don't think we are doing that
bad.
Also in terms of developers working just on symbolics, I think sympy
has much more developers. I don't know if it's easy to extract just
patches to Sage symbolics, to compare speeds of developments, but I
would not be surprised if it turns out it's actually very comparable,
or even less people are working on Sage symbolics than on sympy.

That being said, I like that you pursue the pynac way, because first I
also wanted to use ginac but using swig was just not the way, so
William showed me with cython how to actually do it. And second, I
welcome competition, because that's the only way to actually move
forward, but for Sage and sympy. For example thanks to sympy, you
completely abandoned your previous approach for symbolics, because
after many months of development, it was still slower than pure python
sympy and less tested and robust. So now it's our turn to speed up
sympy. :)

Ondrej

Ondrej Certik

unread,
Apr 22, 2009, 12:09:35 PM4/22/09
to sage-...@googlegroups.com
On Wed, Apr 22, 2009 at 9:08 AM, Ondrej Certik <ond...@certik.cz> wrote:
> welcome competition, because that's the only way to actually move
> forward, but for Sage and sympy. For example thanks to sympy, you

but -> both

O.

root

unread,
Apr 22, 2009, 1:29:04 PM4/22/09
to sage-...@googlegroups.com, sage-...@googlegroups.com, mha...@gmail.com, shu...@gmail.com
> >>> Kudos to SymPy!
> >>>
> >>> I'm wondering why the python integration algorithms implemented there
> >>> aren't in the short term adopted by SAGE.
> >>
> >> They are --- you can use them from sympy inside Sage. It's my goal
> >> that all sympy features are nicely integrated in Sage. I work on this
> >> as time permits.
> >
> > Also, in the pynac-based symbolics that Mike Hansen has been polishing
> > up for full inclusion in Sage (to replace the maxima based symbolics),
> > one can just do
> >
> >   f.integrate(algorithm="sympy")
> >
> > and sage will compute the integral using sympy.   He's already
> > implemented this and I've seen it work well when I tried it out.
>
> Yes, that way Sage allows to call any of those libraries very easily
> in Sage and at the same time uses just one library by default, that
> currently is the best at integrating, I guess still Maxima.

There is an integration test suite at:
http://axiom-developer.org/axiom-website/CATS
which has the Schaums integral series along with examples.
Each integral result is subtracted from the Schaums answer
and then simplified, hopefully to a constant.

How does Sympy compare?

Tim Daly

Tim Lahey

unread,
Apr 22, 2009, 3:43:52 PM4/22/09
to sage-...@googlegroups.com
On Wed, Apr 22, 2009 at 1:29 PM, root <da...@axiom-developer.org> wrote:
>
>
> There is an integration test suite at:
> http://axiom-developer.org/axiom-website/CATS
> which has the Schaums integral series along with examples.
> Each integral result is subtracted from the Schaums answer
> and then simplified, hopefully to a constant.
>
> How does Sympy compare?

As far as I know, it hasn't been tested with it yet. I'm planning support for
it, but I need to do some work in order to get the test suite working with it.
It's on my list.

Cheers,

Tim Lahey.

Ondrej Certik

unread,
Apr 22, 2009, 4:10:12 PM4/22/09
to sage-...@googlegroups.com

It's on our list too, so it will happen eventually. We definitely
still need to improve our algorithms a lot, see e.g.:

http://groups.google.com/group/sympy/browse_thread/thread/58916fb31e1ff1ea

but a nice thing is that it's in Python, so it's easy to work with.

Ondrej

Tim Lahey

unread,
Apr 22, 2009, 4:12:54 PM4/22/09
to sage-...@googlegroups.com
On Wed, Apr 22, 2009 at 4:10 PM, Ondrej Certik <ond...@certik.cz> wrote:
>
>
> It's on our list too, so it will happen eventually. We definitely
> still need to improve our algorithms a lot, see e.g.:
>
> http://groups.google.com/group/sympy/browse_thread/thread/58916fb31e1ff1ea
>
> but a nice thing is that it's in Python, so it's easy to work with.
>
> Ondrej

At some point I probably should get around to joining that group.

If you have any suggestions for the testing architecture, I'd be happy to hear
them.

I've now started using GitHub's issues to do bug/issue tracking for the code.

Cheers,

Tim Lahey.

Maurizio

unread,
Apr 22, 2009, 5:53:26 PM4/22/09
to sage-devel
> We managed to get one gsoc project that does the assumptions right, so
> it may happen anyways over the summer, in fact I very much hope so.
>

How does assumptions affect this? If that's so important, you should
probably get a lot of focus on that! But consider also PDE
important ;)
>
> We'll see. We managed to get 5 gsoc student this summer, so we are not
> at all dead, if this is what you mean. :)
>

That's definitely a very good thing, especially in getting people
involved with development... Probably the best result is not just
their short term, but also the chance of some long term commitment of
good people.

> Besides, it will take almost a year to Sage too to release the new
> symbolics (started August 2008), so I don't think we are doing that
> bad.
> Also in terms of developers working just on symbolics, I think sympy
> has much more developers. I don't know if it's easy to extract just
> patches to Sage symbolics, to compare speeds of developments, but I
> would not be surprised if it turns out it's actually very comparable,
> or even less people are working on Sage symbolics than on sympy.
>

That's a very good point. In my opinion, SAGE community is very good,
and vast, and very well driven, but its focus it's actually quite far
from symbolics (there's nothing bad about it), so I can see that for
the time being, SymPy has some more inertia, thanks to the higher
number of developers committed to symbolics (of course) rather than
SAGE. At the same time, I think that SAGE can make a big step forward,
once its community focuses on this task!

Thanks a lot

Maurizio

Ondrej Certik

unread,
Apr 22, 2009, 6:11:25 PM4/22/09
to sage-...@googlegroups.com
On Wed, Apr 22, 2009 at 1:12 PM, Tim Lahey <tim....@gmail.com> wrote:
>
> On Wed, Apr 22, 2009 at 4:10 PM, Ondrej Certik <ond...@certik.cz> wrote:
>>
>>
>> It's on our list too, so it will happen eventually. We definitely
>> still need to improve our algorithms a lot, see e.g.:
>>
>> http://groups.google.com/group/sympy/browse_thread/thread/58916fb31e1ff1ea
>>
>> but a nice thing is that it's in Python, so it's easy to work with.
>>
>> Ondrej
>
> At some point I probably should get around to joining that group.
>
> If you have any suggestions for the testing architecture, I'd be happy to hear
> them.

We use py.test/nosetest compatible tests, but if you prefer Sage like
doctests, that's fine too.

> I've now started using GitHub's issues to do bug/issue tracking for the code.

I only use github for storing the git repository, haven't used their
bugs/issues yet.

Ondrej

Ondrej Certik

unread,
Apr 22, 2009, 6:20:36 PM4/22/09
to sage-...@googlegroups.com
On Wed, Apr 22, 2009 at 2:53 PM, Maurizio <maurizio...@gmail.com> wrote:
>
>> We managed to get one gsoc project that does the assumptions right, so
>> it may happen anyways over the summer, in fact I very much hope so.
>>
>
> How does assumptions affect this? If that's so important, you should
> probably get a lot of focus on that! But consider also PDE
> important ;)

PDE's are of course important, also it's my main research thing, but
for sympy the assumptions are essential, because you cannot build a
more advanced symbolics without a good assumption system. I am
curious which approach Sage is going to use for this, since ginac
doesn't have any assumptions.

>>
>> We'll see. We managed to get 5 gsoc student this summer, so we are not
>> at all dead, if this is what you mean. :)
>>
>
> That's definitely a very good thing, especially in getting people
> involved with development... Probably the best result is not just
> their short term, but also the chance of some long term commitment of
> good people.

Exactly. Getting people to work with you is absolutely essential,
that's my priority number one. If you work alone, you cannot do
anything in the long term.

Ondrej

Tim Lahey

unread,
Apr 22, 2009, 6:27:23 PM4/22/09
to sage-...@googlegroups.com
On Wed, Apr 22, 2009 at 6:11 PM, Ondrej Certik <ond...@certik.cz> wrote:
>
>
> We use py.test/nosetest compatible tests, but if you prefer Sage like
> doctests, that's fine too.

Good to know. I'm not actually using any unit test classes at the moment
since I have the difficulty that each test has different steps to properly
compare the output to the Schaum's result. The multiple solutions also
pose an additional difficulty.

Cheers,

Tim Lahey.

Ondrej Certik

unread,
Apr 22, 2009, 6:31:05 PM4/22/09
to sage-...@googlegroups.com

I think sympy will do very poorly if the assumptions are needed, we
are still working on the assumptions.

So if it turns out too difficult, just skip sympy for the time being,
we'll get back to it later.

Ondrej

Tim Lahey

unread,
Apr 22, 2009, 6:35:25 PM4/22/09
to sage-...@googlegroups.com
On Wed, Apr 22, 2009 at 6:31 PM, Ondrej Certik <ond...@certik.cz> wrote:
>
>
> I think sympy will do very poorly if the assumptions are needed, we
> are still working on the assumptions.
>
> So if it turns out too difficult, just skip sympy for the time being,
> we'll get back to it later.
>

The problem arises with all the different integration systems. Usually some
kind of simplification is needed on the integral returned, even if there aren't
multiple solutions. This complicates the testing procedure since the steps to
perform the simplification are often specific to the returned result.

I'm currently aiming to finish the test suite just for Sage/Maxima and
I'll go back
and address the various issues (like testing SymPy) once that's complete.

Cheers,

Tim Lahey.

Ondrej Certik

unread,
Apr 22, 2009, 6:51:33 PM4/22/09
to sage-...@googlegroups.com

Ok, awesome. I am interested in all bugs that you find.

Thanks for all the work you are doing,
Ondrej

Carl Witty

unread,
Apr 22, 2009, 6:54:52 PM4/22/09
to sage-...@googlegroups.com
On Wed, Apr 22, 2009 at 3:35 PM, Tim Lahey <tim....@gmail.com> wrote:
> The problem arises with all the different integration systems. Usually some
> kind of simplification is needed on the integral returned, even if there aren't
> multiple solutions. This complicates the testing procedure since the steps to
> perform the simplification are often specific to the returned result.
>
> I'm currently aiming to finish the test suite just for Sage/Maxima and
> I'll go back
> and address the various issues (like testing SymPy) once that's complete.

Would it be better to test the results numerically? (For instance,
evaluate the integral returned and the desired result at 100 random
points to high precision, and ensure that the relative error between
the answers at each point is small.)

Of course, this wouldn't count as a proof that the result was correct,
but IMHO it would be good enough (it seems unlikely that integration
bugs would result in wrong answers that were numerically almost
equivalent to the right answer). (Actually, I might actually trust a
numeric result more than a symbolic simplification-based result, given
the theoretical possibility that a simplification bug might cancel out
an integration bug, leading to a false pass in the test suite;
especially if simplification and integration are done in the same
system.)

Carl

Tim Lahey

unread,
Apr 22, 2009, 6:59:32 PM4/22/09
to sage-...@googlegroups.com
On Wed, Apr 22, 2009 at 6:54 PM, Carl Witty <carl....@gmail.com> wrote:
>
>
> Would it be better to test the results numerically?  (For instance,
> evaluate the integral returned and the desired result at 100 random
> points to high precision, and ensure that the relative error between
> the answers at each point is small.)


How can one do that with symbolic variables? Most of the test integrals
have symbolic constants (w.r.t the integration) so it isn't just the
integration
variable. I thought about numerical testing, but it isn't generally feasible.

>
> Of course, this wouldn't count as a proof that the result was correct,
> but IMHO it would be good enough (it seems unlikely that integration
> bugs would result in wrong answers that were numerically almost
> equivalent to the right answer).  (Actually, I might actually trust a
> numeric result more than a symbolic simplification-based result, given
> the theoretical possibility that a simplification bug might cancel out
> an integration bug, leading to a false pass in the test suite;
> especially if simplification and integration are done in the same
> system.)

It's possible to have simplification bugs, but I'll have to rely upon separate
tests of the simplification system.

Cheers,

Tim Lahey.

Carl Witty

unread,
Apr 22, 2009, 7:03:12 PM4/22/09
to sage-...@googlegroups.com
On Wed, Apr 22, 2009 at 3:59 PM, Tim Lahey <tim....@gmail.com> wrote:
>
> On Wed, Apr 22, 2009 at 6:54 PM, Carl Witty <carl....@gmail.com> wrote:
>>
>>
>> Would it be better to test the results numerically?  (For instance,
>> evaluate the integral returned and the desired result at 100 random
>> points to high precision, and ensure that the relative error between
>> the answers at each point is small.)
>
>
> How can one do that with symbolic variables? Most of the test integrals
> have symbolic constants (w.r.t the integration) so it isn't just the
> integration
> variable. I thought about numerical testing, but it isn't generally feasible.

Couldn't you just pick random values for all of the symbolic constants, as well?

Carl

Tim Lahey

unread,
Apr 22, 2009, 7:11:22 PM4/22/09
to sage-...@googlegroups.com
On Wed, Apr 22, 2009 at 7:03 PM, Carl Witty <carl....@gmail.com> wrote:
>
>
> Couldn't you just pick random values for all of the symbolic constants, as well?

Yes, but over what range? If you do that, you've just ensured that it
is correct for
those points. It also could get expensive if you have multiple
constants combined
with the integration variable. It's an option to consider, but I'd
probably only consider
for those that don't simplify to the same result. It could also be
done as an optional
flag. I'll add it to the list of things to consider.

If someone wants to add issues/bugs for me, it can be done at:

http://github.com/tjl/sage_int_testing/issues

It requires a GitHub account (which is free). Or one can e-mail me
directly and I'll
add it.

Cheers,

Tim.

root

unread,
Apr 22, 2009, 5:13:51 PM4/22/09
to sage-...@googlegroups.com, sage-...@googlegroups.com
> > It's on our list too, so it will happen eventually. We definitely
> > still need to improve our algorithms a lot, see e.g.:
> >
> > http://groups.google.com/group/sympy/browse_thread/thread/58916fb31e1ff1ea
> >
> > but a nice thing is that it's in Python, so it's easy to work with.
> >
> > Ondrej
>
> At some point I probably should get around to joining that group.
>
> If you have any suggestions for the testing architecture, I'd be happy to hear
> them.
>
> I've now started using GitHub's issues to do bug/issue tracking for the code.

What's the github name?

Tim Daly

Tim Lahey

unread,
Apr 22, 2009, 9:40:29 PM4/22/09
to sage-...@googlegroups.com, root, Tim Lahey

The issue tracker:
http://github.com/tjl/sage_int_testing/issues

and the main repository:
http://github.com/tjl/sage_int_testing/tree/master

Cheers,

Tim.

---
Tim Lahey
PhD Candidate, Systems Design Engineering
University of Waterloo
http://www.linkedin.com/in/timlahey

Jason Grout

unread,
Apr 23, 2009, 3:51:02 PM4/23/09
to sage-...@googlegroups.com
Ondrej Certik wrote:
> On Wed, Apr 22, 2009 at 2:53 PM, Maurizio <maurizio...@gmail.com> wrote:
>>> We managed to get one gsoc project that does the assumptions right, so
>>> it may happen anyways over the summer, in fact I very much hope so.
>>>
>> How does assumptions affect this? If that's so important, you should
>> probably get a lot of focus on that! But consider also PDE
>> important ;)
>
> PDE's are of course important, also it's my main research thing, but
> for sympy the assumptions are essential, because you cannot build a
> more advanced symbolics without a good assumption system. I am
> curious which approach Sage is going to use for this, since ginac
> doesn't have any assumptions.


Is there anyone working on an assumptions framework in Sage? On the one
hand, you can work in certain domains in Sage (i.e., polynomials over
QQ, etc.), so some of this need may be taken care of there.

But for a more general framework (like declaring that x>0?)
Hmmm...maybe we'll use sympy? :)

Jason

--
Jason Grout

William Stein

unread,
Apr 23, 2009, 4:07:41 PM4/23/09
to sage-...@googlegroups.com
On Thu, Apr 23, 2009 at 12:51 PM, Jason Grout
<jason...@creativetrax.com> wrote:
>
> Ondrej Certik wrote:
>> On Wed, Apr 22, 2009 at 2:53 PM, Maurizio <maurizio...@gmail.com> wrote:
>>>> We managed to get one gsoc project that does the assumptions right, so
>>>> it may happen anyways over the summer, in fact I very much hope so.
>>>>
>>> How does assumptions affect this? If that's so important, you should
>>> probably get a lot of focus on that! But consider also PDE
>>> important ;)
>>
>> PDE's are of course important, also it's my main research thing, but
>> for sympy the assumptions are essential, because you cannot build a
>> more advanced symbolics without a good assumption system.  I am
>> curious which approach Sage is going to use for this, since ginac
>> doesn't have any assumptions.

Could you explain how assumptions are so important? Could you
particularly address how they can (1) be so critically important, and
yet (2) ginac doesn't have them. Incidentaly, to me they are
particularly important in symbolic integration, which ginac doesn't
do. Also, could you explain why the assumption system in Sympy
needs to be rewritten -- in particular, what design decisions were
suboptimal?

-- William


>
>
> Is there anyone working on an assumptions framework in Sage?  On the one
> hand, you can work in certain domains in Sage (i.e., polynomials over
> QQ, etc.), so some of this need may be taken care of there.
>
> But for a more general framework (like declaring that x>0?)
> Hmmm...maybe we'll use sympy? :)
>
> Jason
>
> --
> Jason Grout
>
>
> >
>



Tim Lahey

unread,
Apr 23, 2009, 4:17:50 PM4/23/09
to sage-...@googlegroups.com

On Apr 23, 2009, at 4:07 PM, William Stein wrote:

>
> Could you explain how assumptions are so important? Could you
> particularly address how they can (1) be so critically important, and
> yet (2) ginac doesn't have them. Incidentaly, to me they are
> particularly important in symbolic integration, which ginac doesn't
> do. Also, could you explain why the assumption system in Sympy
> needs to be rewritten -- in particular, what design decisions were
> suboptimal?
>
> -- William
>

One place where I use assumptions regularly is with trig functions. If
you have sin(n*pi) or cos(n*pi), stating the assumption that n is
integer
reduces the first to 0 and the second to alternating +1,-1. These crop
up
in modal analysis of physical systems regularly. I don't know
if sympy's assumptions are suboptimal, since I haven't used sympy much,
but I guess it is since they're planning on improvements.

Ondrej Certik

unread,
Apr 23, 2009, 4:43:54 PM4/23/09
to sage-...@googlegroups.com
On Thu, Apr 23, 2009 at 1:17 PM, Tim Lahey <tim....@gmail.com> wrote:
>
>
> On Apr 23, 2009, at 4:07 PM, William Stein wrote:
>
>>
>> Could you explain how assumptions are so important?  Could you

We already discussed this many times on this list, just search the
archives. Without good assumptions, you cannot implement good
integration, good solvers, good simplification, nothing. Of course,
things like x+x is easy, but once you have for example
Integral(<complex expression involving many constants>), sometimes you
can simplify it, sometimes not and this very much depends on the
assumptions for the constants. Or things like integrate(x**n, x).

>> particularly address how they can (1) be so critically important, and
>> yet (2) ginac doesn't have them.  Incidentaly, to me they are

ginac doesn't have any assumptions and so ginac doesn't have any
advanced symbolic features.
Of course, if the only thing that you are going to use ginac for is
the core engine, then it should not matter much. But what I want from
a CAS is to be able to do advanced symbolic calculations with symbols,
so I need some way to tell it my assumptions about the symbols. Taking
everything as complex numbers will not simplify things in many cases.

>> particularly important in symbolic integration, which ginac doesn't
>> do.    Also, could you explain why the assumption system in Sympy
>> needs to be rewritten -- in particular, what design decisions were
>> suboptimal?

Because we attach assumptions to symbols, e.g. you define

In [2]: x = Symbol("x", positive=True)

and then you work with that:

In [3]: sqrt(x**2)
Out[3]: x

In [4]: x = Symbol("x", negative=True)

In [5]: sqrt(x**2)
Out[5]: -x

In [6]: x = Symbol("x")

In [7]: sqrt(x**2)
Out[7]:
⎽⎽⎽⎽
╱ 2
╲╱ x


But this approach is the wrong one, because then the core engine has
to take this into account and it slows things down. Our new system
will keep the assumptions external, and define methods to work with
them, like in mathematica, e.g. refine(). This should simplify the
core and then we should be able to speed it up.

>>
>> -- William
>>
>
> One place where I use assumptions regularly is with trig functions. If
> you have sin(n*pi) or cos(n*pi), stating the assumption that n is
> integer
> reduces the first to 0 and the second to alternating +1,-1. These crop
> up
> in modal analysis of physical systems regularly. I don't know
> if sympy's assumptions are suboptimal, since I haven't used sympy much,
> but I guess it is since they're planning on improvements.

Exactly, all kinds of these assumptions, that are necessary for
integration, simplifications, e.g. sometimes sqrt(x**2) reduces to x,
sometimes to -x, sometimes to neither.

Ondrej

Burcin Erocal

unread,
Apr 24, 2009, 6:55:11 AM4/24/09
to sage-...@googlegroups.com

Note that ginac has a similar mechanism, which it uses for the
automatic simplifications it does. I haven't exposed this in the Sage
interface. I don't know if Mike has either.

Here is a list of the info flags ginac supports (Scroll down for the
table):

http://www.ginac.de/tutorial/Information-about-expressions.html

Pynac has an "infinity" flag in addition to the ones listed there.


I admit that I have no experience with general simplification
routines, but I have the feeling that it all boils down to answering
is_positive(), is_integer() queries about expressions. The trick is to
be able to deduce the answers to these from a few bits of information
gathered from the user.

I don't know of any open source package which can provide this
capability for Sage. I've heard that Reduce with its Redlog package is
very powerful in this area. I doubt if it can be used from Sage though.


Cheers,
Burcin

Reply all
Reply to author
Forward
0 new messages