Re: [sage-support] Problem doing symbolic computations (bug in Pynac ?)

9 views
Skip to first unread message

Burcin Erocal

unread,
Sep 7, 2010, 5:26:39 PM9/7/10
to sage-s...@googlegroups.com, pynac...@googlegroups.com
Hi,

Here is a short example to replicate the first error mentioned below:

b = [var('b_%s'%i) for i in range(4)]

precomp = (2^b_2 + 2)*(2^b_1 + 2^(-b_1) + 2^b_1*2^b_0 - 2^b_1*2^(-b_0)
- 2^(-b_1)*2^b_0 - 2^(-b_1)*2^(-b_0) + 2^b_0 + 2^(-b_0) - 9) + (2^b_1 +
2^(-b_1) + 2^b_1*2^b_0 - 2^b_1*2^(-b_0) - 2^(-b_1)*2^b_0 -
2^(-b_1)*2^(-b_0) + 2^b_0 + 2^(-b_0) - 9)/2^b_2

repl_dict = {b_0: b_0, b_3: b_1, b_2: b_3, b_1: b_2}
P = precomp.substitute(repl_dict)
P.expand()


I will take a break now. I'd appreciate any help tracking down the
problem (in pynac) if anybody has the time.


Cheers,
Burcin


On Wed, 1 Sep 2010 08:25:07 -0700 (PDT)
Jean-Pierre Flori <jpf...@gmail.com> wrote:

> Hello,
>
> I'm using Sage to do some symbolic computations.
> What I am basically trying to do is to compute sums made of terms of
> the form :
> cst_{i_0,...,i_d)*b_0^{i_0}*...*b_d^{i_d}*2^{-b_0-...-b_d}
> where the b_i are variables
> (i.e. a product of an exponential with b_i in the exponent by a
> multivariate polynomial in b_i)
> using functions (especially symbolic_sum) from sage.calculus.calculus.
> I am regularly calling expand() and simplify() functions as well as
> some custom functions using substitute() to get a more compact
> expression and avoid redundant terms.
> Everything works fine as long as the sums are not too "complicated",
> but Sage segfaults when they become too big.
> Running Sage with the -gdb flag, gives different errors according to
> which of my functions I am calling, but both errors seem to be related
> to a "compare" function in Pynac (SIGSEGV or SIGBUS):
>
> First error message :
> Program received signal SIGSEGV, Segmentation fault.
> GiNaC::power::compare (this=0x3281c20, other=...) at power.cpp:899
> 899 power.cpp: Aucun fichier ou dossier de ce type.
> in power.cpp
>
> Second error message :
> Program received signal SIGBUS, Bus error.
> GiNaC::mul::compare (this=0x83490e0, other=...) at mul.cpp:1174
> 1174 mul.cpp: Aucun fichier ou dossier de ce type.
> in mul.cpp
>
>
> My box is running a Debian unstable with packages from the
> experimental repository.
> My kernel version is 2.6.35, compiled by myself.
> My cpu is a Intel(R) Core(TM)2 Quad CPU Q6600.
> The error is happening with Sage 4.5.2 and 4.4.4.
> I compiled both of them myself with gcc (Debian 4.4.4-11) 4.4.5
> 20100824 (prerelease).
> I also tried the precompiled sage-4.5.2-linux-64bit-ubuntu_10.04_lts-
> x86_64-Linux which gives me the same result.
>
> I didn't manage to reproduce the bug with a short and simple piece of
> code.
> However I can mail the source code and instructions to reproduce the
> bug to anyone who can help.
> I will also provide any information needed to investigate this
> problem.
> I am currently trying to compile Sage on another box to do another
> test.
> Any help will be appreciated, and thanks for a such a great piece of
> software.
>
> Best regards,
>
> Jean-Pierre Flori
>

Burcin Erocal

unread,
Sep 12, 2010, 9:40:29 AM9/12/10
to sage-s...@googlegroups.com, pynac...@googlegroups.com
Hi Jean-Pierre,

On Thu, 9 Sep 2010 02:13:27 -0700 (PDT)
Jean-Pierre Flori <jpf...@gmail.com> wrote:

> I created a Ticket on Sage's Trac:
> http://trac.sagemath.org/sage_trac/ticket/9880
>
> It could be a problem with the comparison function
> "expair_rest_is_less" used by "std::sort" function like in the
> following link:
> http://gcc.gnu.org/ml/gcc-bugs/2010-08/msg01163.html

This is indeed the problem. Thanks for tracking it down.

In the process of adapting GiNaC for use in Sage, I changed the
comparison functions to give a consistent order between different runs.
Apparently, the new comparison functions I wrote are not consistent.

For a long time, I wanted to rip out these new functions, and use them
only for printing. The following patch is a first step towards this:

http://sage.math.washington.edu/home/burcin/pynac/pynac_order.patch

After applying this patch to pynac (see here [1] for instructions on
how to set up a pynac development environment), the segfault you
reported, and the bug at #9046 goes away.

Note that this is very much a work in progress. I didn't take care to
make sure these comparison functions are consistent either. In the
process of copying things over, I even flipped a few signs in the
results.

This needs much work to be finished, we need to make sure the order
is consistent and correct, and it matches the one currently used by
Sage. The issue of accessing the operands of a symbolic expression also
needs to be sorted out, but that's merely a user interface issue.

I'd appreciate help with any of this, since I don't have much time to
spend on Sage/pynac these days. I'd be glad to answer any questions
though.


Cheers,
Burcin

Burcin Erocal

unread,
Sep 12, 2010, 12:36:26 PM9/12/10
to sage-s...@googlegroups.com, pynac...@googlegroups.com
Hi,

On Sun, 12 Sep 2010 09:21:12 -0700 (PDT)
Jean-Pierre Flori <jpf...@gmail.com> wrote:

> Thanks for the patch, I'll have a look at the comparison functions
> when I have some time, surely next week.
> Where can I find the one defined by Sage ?

They are spread out through the classes in pynac. Look for compare()
and compare_same_type() functions in mul, power, function, etc.

> I have another question: are the source codes of GiNaC and pynac still
> "synchronized" in some way ?

Not really. I pull some patches from GiNaC when I think they are
relevant.

> I saw that some recent changes are similar in both git repositories,
> but as I was going through the source code, I saw that
> numeric::do_print_csrc is defined in GiNaC but not in pynac.

The numeric class in pynac is completely rewritten to use python
objects instead of CLN. It's normal that those methods are missing.

Cheers,
Burcin

Burcin Erocal

unread,
Sep 29, 2010, 2:28:14 PM9/29/10
to sage-s...@googlegroups.com, pynac...@googlegroups.com
On Wed, 29 Sep 2010 09:48:25 -0700 (PDT)
Jean-Pierre Flori <jpf...@gmail.com> wrote:

> Maybe it is a good thing to keep the same order as ginac internally
> and your more usual ordering for printing.

It is good to keep the ginac ordering internally. The user friendly
ordering is more expensive so it slows down computations, even when
you don't print intermediate expressions.

> However if you'd better not duplicate code, I can look at the "-
> x^2+x^2" part of bug #9046.
> Now I may understand a big enough part of pynac code to do that.

That's right. You're one of the few people who spent so much time
staring at pynac code. :)

> But if you'd better use the above patch, that won't be necessary.

I would like to replace the orders completely. I wanted to attack
this problem for a long time. Now that we started, we might as
well finish. :)

Can you help put the patch in usable shape? At the moment the ordering
it defines is completely different from the one we use in Sage. This
must be because I changed some signs while I was copying code.

If you can make the ordering consistent, I can fix the other relevant
places (op(), etc.).


Thanks.

Burcin

Burcin Erocal

unread,
Oct 1, 2010, 2:59:46 PM10/1/10
to Jean-Pierre Flori, pynac...@googlegroups.com
Hi,

This discussion should really be on pynac-devel. I set the Reply-To
address accordingly. I hope it works. ;)

On Fri, 1 Oct 2010 06:52:11 -0700 (PDT)
Jean-Pierre Flori <jpf...@gmail.com> wrote:

> Some remarks and questions...
>
> I guess the order you want is 'degrevlex' as Sage's default order for
> multivariate polynomial ring and as the name of your functions
> suggest.
> Am I wrong ?

You are right, some approximation of 'degrevlex'.

> At present, I did not modify a lot of things.
> A few signs, replaced seq[0] by the smallest item of the sequence
> sorted with the printing order.

Nice. I totally forgot about having to sort subexpressions.

> I guess we want symbolic variables to be bigger than numerics.
> So I changed the sign for different typeid's.

OK.

> What is your opinion about symbolic exponents ?
> (If I am not wrong, in mul object, they are encapsulated in power
> object so it counts for 1.
> So total_degree(a^b*x) == 2.)
> In the compare_mul_power, you decided to compare the exponent of the
> power object to the one of the smallest item (which could in fact be
> the same power object) of the mul object if it is symbolic.
> So 3*x^y will be smaller than x^y (considering that y > 1).
> That is not a "real" problem because 3*x^y+x^y will be simplified into
> 4*x^y.
> But we also get that a^b > a^b*x (if b > 1) with the current strategy
> and it prints a^b + a^b*x...

It looks like you just need to return a^b < a^b*x. This is also what we
have now:

sage: var('a,b,x')
(a, b, x)
sage: a^b + a^b*x
a^b*x + a^b

Note that the patch reverses the order used for printing. AFAIR, it
used to be that we said x^3 < x^2 to be able to print x^3 + x^2.

> Shouldn't we consider the degree of the power is 1, so at least the
> above case is solved.
> Or do something more complicated comparing exponents...

If you always take the degree of power objects to be 1, you end up
saying a^b == a^c. Something more complicated is required there...

> I do not understand the comment of William Stein about symbols
> ordering.
> When a polynomial ring is created in Sage, the variables are ordered
> according to the order they were "created", i.e.
> PR.<x,y> = PolynomialRing(QQ) will result in x > y
> and PR.<y,x> = PolynomialRing(QQ) in y > x

He's probably trying to describe this order:

sage: var('z,x,y,a')
(z, x, y, a)
sage: a+x+y+z
a + x + y + z


Cheers,
Burcin

Jean-Pierre Flori

unread,
Oct 2, 2010, 6:30:27 AM10/2/10
to pynac-devel
Hi,

I managed (?) to get mercurial to produce a patch to be applied after
yours at :
http://perso.telecom-paristech.fr/~flori/pynac_order-2.patch
Hope that step did not introduce any problem.

I left the code for comparing symbols as is (i.e. it reverse the
string comparison of C++).
And I did not have much time to test it, but it is a first step.

Cheers,

Jean-Pierre Flori

unread,
Oct 1, 2010, 5:05:10 PM10/1/10
to pynac-devel
Hi,

On 1 oct, 20:59, Burcin Erocal <bur...@erocal.org> wrote:
> Hi,
>
> This discussion should really be on pynac-devel. I set the Reply-To
> address accordingly. I hope it works. ;)
Thanks !
> It looks like you just need to return a^b < a^b*x. This is also what we
> have now:
>
> sage: var('a,b,x')
> (a, b, x)
> sage: a^b + a^b*x
> a^b*x + a^b
As I changed some function it turned out to be the converse because a
was greater than x (don't remember why, I must have changed something
somewhere...), so it compared to compare the exponent in a^b (b) and
the exponent of x (1) which returned b > 1 and a^b > x so that I got
a^b*x < a^b

>
> Note that the patch reverses the order used for printing. AFAIR, it
> used to be that we said x^3 < x^2 to be able to print x^3 + x^2.
>
> > Shouldn't we consider the degree of the power is 1, so at least the
> > above case is solved.
> > Or do something more complicated comparing exponents...
>
> If you always take the degree of power objects to be 1, you end up
> saying a^b == a^c. Something more complicated is required there...
>
For such expression, it is the compare_s_t_power function which is
called and that one will compare the exponents, so that should not be
a problem.

> > I do not understand the comment of William Stein about symbols
> > ordering.
> > When a polynomial ring is created in Sage, the variables are ordered
> > according to the order they were "created", i.e.
> > PR.<x,y> = PolynomialRing(QQ) will result in x > y
> > and PR.<y,x> = PolynomialRing(QQ) in y > x
>
> He's probably trying to describe this order:
>
> sage: var('z,x,y,a')
> (z, x, y, a)
> sage: a+x+y+z
> a + x + y + z
>
I think you are right.
My concern is that it is not the kind of lexicographic order that is
used by PolynomialRing, but that is just a detail.
Both should be ok, but one has to be chosen.

I have working code for printing.
I must now produce a clean patch file (which looks harder to me) so
other people can test it.

Cheers,

Jean-Pierre Flori

unread,
Oct 4, 2010, 1:00:51 PM10/4/10
to pynac-devel
Hi,

I put a different version of the patch at the same adress.
It treats comparison differently according to whether we are in a mul
or an add object.
However I did it in a kind of dirty way :
- virtualized get_sorted_seq, overrided it in mul objects...
- ..using a flag to know if we are inside a mul or add object.
Don't know how you would like to do it.

As far as symbolic exponents are concerned, we could add an exponent
member to mul object and compute it when we want to print as you did
with sorted_seq.
We would use it to do a good comparison of exponents.
Of course that would be slower...
Right now the total_degree function does not make much sense because
of the symbolic exponents within power objects which count for 1 and
we still have :
sage: x*a^b+a^b
a^b + a^b*x
and inconsistent ordering for expressions mixing power objects and mul
objects with different kind of exponents :
sage: x^(a+b)+x^3 + x^4+a^4+a^2*x^2
a^2*x^2 + x^(a + b) + a^4 + x^4 + x^3
because I always compare exponents for power objects, but not for mul
and power objects.
If you agree with that idea, i can do it and we would have something
much more consistent.

Cheers,

Burcin Erocal

unread,
Oct 6, 2010, 3:15:06 AM10/6/10
to pynac...@googlegroups.com
Hi Jean-Pierre,

Thank you for the patch. I don't have time to test it right now, but
at least it applies cleanly on top of mine.

I will be offline for a couple of days at a meeting by a lake side
conference center in Austria, which unfortunately doesn't have wifi. I
should be able to respond to your message properly when I get back in
the weekend.

BTW, does your patch add a comparison operator for functions? We could
sort them alphabetically based on the function name for example. This
would resolve #9632.


Cheers,
Burcin

Jean-Pierre Flori

unread,
Oct 6, 2010, 4:15:27 AM10/6/10
to pynac-devel
Hi,

Answer when you can, no problem, I also have lots of urgent things to
do right now...

- I fear I introduced some new inconsistencies while tweaking
everything, I will run more tests, but I posted a irst version of the
patch to get some feedback.
- I didn't think of sorting the functions because i don't use any,
I'll add that.
- I'll also add functions to get "real" total exponents in my previous
post and (try to) properly compare symbolic (or not with hold ?)
exponents of pow objects wrapped inside mul objects.
- I did not take into account the fact that expressions can now be
held, so maybe some of the assumptions I made are not true anymore (no
mul inside a mul e.g., everything is noted in comments whithin the
source code), I'll try to get what hold exactly does, test and fix it.

Cheers,

Jean-Pierre Flori

unread,
Oct 7, 2010, 6:39:57 AM10/7/10
to pynac-devel
Hi,

> - I didn't think of sorting the functions because i don't use any,
> I'll add that.
A new version is up with code for functions and fderivatives (http://
perso.telecom-paristech.fr.fr/~flori/pynac_order-2.patch).

> - I'll also add functions to get "real" total exponents in my previous
> post and (try to) properly compare symbolic (or not with hold ?)
> exponents of pow objects wrapped inside mul objects.
Not done yet, will do it later, surely not before end of october.

> - I did not take into account the fact that expressions can now be
> held, so maybe some of the assumptions I made are not true anymore (no
> mul inside a mul e.g., everything is noted in comments whithin the
> source code), I'll try to get what hold exactly does, test and fix it.
Not done yet, working on it soon.
I made some comments in a new thread.

Cheers,

Jean-Pierre Flori

unread,
Oct 8, 2010, 6:03:04 AM10/8/10
to pynac-devel
By the way, looking at http://trac.sagemath.org/sage_trac/ticket/10053
, I realized that the restored GiNaC order (so mainly comparing hash
values) is called when one calls cmp(a,b) or a.__cmp__(b).
I guess that should also be changed to use the order used for printing.

Jean-Pierre Flori

unread,
Oct 8, 2010, 8:00:29 PM10/8/10
to pynac-devel
The version up now seems in a good shape.
http://perso.telecom-paristech.fr/pynac_order-2.patch

The few tests I ran for held expressions were ok.

I forgot the idea of computing total exponents, that would slow
everything down and would only introduce problems.
The comparison is made kind of like before : symbolic exponents are
wrapped inside power objects so that they count for 1 in total
degrees, and those power kind of act like "new" variables.
We get x^2 > x^a > x^b > x for example.

The main thing I'm not completely happy with is the way we compare
complex numbers with norm.
We could maybe use GiNaC function.
However when I tried it with hold, I got strange errors :
HIT STUB : Number_T_to... unable...
happening in numeric.cpp.
The code I used is still within order.cpp but commented out.

Burcin Erocal

unread,
Oct 10, 2010, 6:34:24 PM10/10/10
to pynac...@googlegroups.com
Hi,

On Fri, 8 Oct 2010 17:00:29 -0700 (PDT)
Jean-Pierre Flori <jpf...@gmail.com> wrote:

> The version up now seems in a good shape.
> http://perso.telecom-paristech.fr/pynac_order-2.patch

This address gives a 404 error for me. Can you attach the patch to the
trac ticket?

> The few tests I ran for held expressions were ok.
>
> I forgot the idea of computing total exponents, that would slow
> everything down and would only introduce problems.
> The comparison is made kind of like before : symbolic exponents are
> wrapped inside power objects so that they count for 1 in total
> degrees, and those power kind of act like "new" variables.
> We get x^2 > x^a > x^b > x for example.

At the moment, we have this:

sage: var('a,b')
(a, b)
sage: x^2 + x^a + x^b + x
x^b + x^a + x^2 + x

How come the new order is different? Note that all these are power
objects, so the comparison functions of mul are not relevant here.

> The main thing I'm not completely happy with is the way we compare
> complex numbers with norm.
> We could maybe use GiNaC function.
> However when I tried it with hold, I got strange errors :
> HIT STUB : Number_T_to... unable...
> happening in numeric.cpp.
> The code I used is still within order.cpp but commented out.

I can take a look at this once I download the patch.

Thank you.

Burcin

Burcin Erocal

unread,
Oct 10, 2010, 6:58:28 PM10/10/10
to pynac...@googlegroups.com

I guess we'll have to do this.

Using a consistent order for __cmp__() is enough, so the printing order
is not necessary. Keeping the GiNaC order for __cmp__() would make it
much faster, though this will be

* very confusing to users,
* not consistent between different runs of Sage

There is not much of a choice really.

We will also have to change the operands(), iterator(), etc. functions
to address the operands in the printing order.

Cheers,
Burcin

Jean-Pierre Flori

unread,
Oct 11, 2010, 2:23:55 AM10/11/10
to pynac-devel


On 11 oct, 00:34, Burcin Erocal <bur...@erocal.org> wrote:
> Hi,
>
> On Fri, 8 Oct 2010 17:00:29 -0700 (PDT)
>
> Jean-Pierre Flori <jpfl...@gmail.com> wrote:
> > The version up now seems in a good shape.
> >http://perso.telecom-paristech.fr/pynac_order-2.patch
>
> This address gives a 404 error for me. Can you attach the patch to the
> trac ticket?
>
I'll post both your patch and mine on the trac server asap.

> > I forgot the idea of computing total exponents, that would slow
> > everything down and would only introduce problems.
> > The comparison is made kind of like before : symbolic exponents are
> > wrapped inside power objects so that they count for 1 in total
> > degrees, and those power kind of act like "new" variables.
> > We get x^2 > x^a > x^b > x for example.
>
> At the moment, we have this:
>
> sage: var('a,b')
> (a, b)
> sage: x^2 + x^a + x^b + x
> x^b + x^a + x^2 + x
>
> How come the new order is different? Note that all these are power
> objects, so the comparison functions of mul are not relevant here.
>
I was trying to get something logical (or at least consistent) about
total degrees.
Previously we got x^b*y < y^3 and x^b*x^2 < x^4 but x^b > x^4...
So my idea was to consider "x^b" as a "new symbol", bigger than "x"
and of degree "1".
So we still have x^b*y < y^3 and x^b*x^2 < x^4 but now x^b < x^4.

JP
Reply all
Reply to author
Forward
0 new messages