Comparison of sparse polynomials

[bruno]

Currently, the default _cmp_ is written with dense polynomials in mind since it creates an xrange of size the degree of the polynomial. For sparse polynomials which can have very high degree, this is a problem:

sage: R.<x> = PolynomialRing(ZZ, sparse=True)
sage: x^2^100 == x^2^100
Traceback (most recent call last)
...
OverflowError: Python int too large to convert to C long

The patch introduces a _cmp_ method in the class Polynomial_generic_sparse to get rid of this problem:

sage: R.<x> = PolynomialRing(ZZ, sparse=True)
sage: x^2^100 == x^2^100
True

Note, as a side comment, that the comparison is faster for polynomials that are indeed sparse with the new code, while the timings are equivalent for quite dense polynomials.

Without the patch:

sage: %timeit x^100 == x^100
1000 loops, best of 3: 224 µs per loop

With the patch:

sage: %timeit x^100 == x^100
The slowest run took 7.61 times longer than the fastest. This could mean that an intermediate result is being cached 
10000 loops, best of 3: 24.6 µs per loop

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​421bd04

Add _cmp_ in class Polynomial_generic_sparse

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​7c72ce2

Rewritten using lists to correct a bug

[bruno]

There was a bug in my first commit, for some case I didn't think to test. I hope there is no more bug, and I've tested all possibilities.

[vdelecroix]

Hello,

1. You can replace

   keys1 = list(reversed(sorted(d1.keys())))
   

   with

   keys1 = d1.keys()
   keys1.reverse()
   

2. Instead of looking at the length, you can do

   if not keys1 and not keys2: return 0
   if not keys1: return -1
   if not keys2: return 1
   

3. Why the first key is treated differently?

Vincent

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​1726058

list(reversed(sorted(...))) replaced by .sort() and .reverse()

[bruno]

Replying to vdelecroix:

1. You can replace

   keys1 = list(reversed(sorted(d1.keys())))
   

   with

   keys1 = d1.keys()
   keys1.reverse()
   

Done, with the added keys1.sort() before the reverse().

2. Instead of looking at the length, you can do

   if not keys1 and not keys2: return 0
   if not keys1: return -1
   if not keys2: return 1
   

Also done.

3. Why the first key is treated differently?

This is to comply with the default _cmp_ method for polynomials: for instance -2x^2 > x but x^3 - 2x^2 < x^3 + x. A larger degree polynomial is always greater than a lower degree polynomial, but then when the coefficients are compared, 0 is greater than a negative coefficient.

I add a question: I noticed after your remark than the method exponents() could be improved. This is now implemented as

        return [c[0] for c in sorted(self.__coeffs.iteritems())] 

while it would be faster (from my testing) to replace it by

        keys = self.__coeffs.keys()
        keys.sort()
        return keys

Is that right to make the change in this ticket (since it is a rather minor change) or is it more appropriate to open a new one?

[vdelecroix]

Replying to bruno:

Replying to vdelecroix: I add a question: I noticed after your remark than the method exponents() could be improved. This is now implemented as

        return [c[0] for c in sorted(self.__coeffs.iteritems())] 

while it would be faster (from my testing) to replace it by

        keys = self.__coeffs.keys()
        keys.sort()
        return keys

Is that right to make the change in this ticket (since it is a rather minor change) or is it more appropriate to open a new one?

That's fine with me.

And could you add a doctest like

sage: Rd = PolynomialRing(ZZ, 'x', sparse=False)
sage: Rs = PolynomialRing(ZZ, 'x', sparse=True)
sage: for _ in range(100):
....:     pd = Rd.random_element()
....:     qd = Rd.random_element()
....:     assert cmp(pd,qd) == cmp(Rs(pd), Rs(qd))

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​efa67dd	Fasten Polynomial_generic_sparse.exponents()
​c9b2091	Add doctest

[vdelecroix]

keys1.sort(); keys1.reverse() -> keys1.sort(reverse=True)

[vdelecroix]

In

            if c > 0:
                c1 = cmp(d1[k1], 0)
                if c1: return c1

Isn't c1 automatically nonzero? Could there be some zero coefficients in the dictionary?

[vdelecroix]

Instead of comparing with 0 (the Python int zero) it would be better to compare with self.base_ring().zero().

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​1bb00ff	sort + reverse in a single instruction
​3121374	Compare with self.base_ring().zero() instead of Python 0

[bruno]

Replying to vdelecroix:

In

            if c > 0:
                c1 = cmp(d1[k1], 0)
                if c1: return c1

Isn't c1 automatically nonzero? Could there be some zero coefficients in the dictionary?

sage: R.<x> = PolynomialRing(ZZ, sparse=True)
sage: p = R({0:1, 1:0, 2:1}, check=False)
sage: p
x^2 + 1
sage: p.dict()
{0: 1, 1: 0, 2: 1}

I doubt this happens often in practice, but since this is possible, I am quite reluctant to remove the test.

[vdelecroix]

Replying to bruno:

Replying to vdelecroix:

In

            if c > 0:
                c1 = cmp(d1[k1], 0)
                if c1: return c1

Isn't c1 automatically nonzero? Could there be some zero coefficients in the dictionary?

sage: R.<x> = PolynomialRing(ZZ, sparse=True)
sage: p = R({0:1, 1:0, 2:1}, check=False)
sage: p
x^2 + 1
sage: p.dict()
{0: 1, 1: 0, 2: 1}

I doubt this happens often in practice, but since this is possible, I am quite reluctant to remove the test.

Dough! I was asking because for sparse matrices it is explicitely mentioned in the specification that no value can be zero. And this is very useful!

[vdelecroix]

Though, multiplication is nicer

sage: (p*p).dict()
{0: 1, 2: 2, 4: 1}

What do you think about opening a ticket about enforcing non-zero values in the dictionary?

Vincent

[bruno]

Replying to vdelecroix:

Dough! I was asking because for sparse matrices it is explicitely mentioned in the specification that no value can be zero. And this is very useful!

There is no such specification as for matrices. There is a method __normalize() to remove the zero coefficients, used by default in __init__, and used for instance in the multiplication (since there may be some cancellations).

Replying to vdelecroix:

Though, multiplication is nicer

sage: (p*p).dict()
{0: 1, 2: 2, 4: 1}

What do you think about opening a ticket about enforcing non-zero values in the dictionary?

I am in favor of opening a ticket to add this specification, though I am not in favor of removing the check=False option that can be useful at some places to get (slightly) faster code.

[vdelecroix]

Replying to bruno:

Replying to vdelecroix:

Dough! I was asking because for sparse matrices it is explicitely mentioned in the specification that no value can be zero. And this is very useful!

There is no such specification as for matrices. There is a method __normalize() to remove the zero coefficients, used by default in __init__, and used for instance in the multiplication (since there may be some cancellations).

matrix_generic_sparse.pyx, main docstring line 69-71

    A generic sparse matrix is
    represented using a dictionary whose keys are pairs of integers `(i,j)` and
    values in the base ring. The values of the dictionary must never be zero.

Replying to vdelecroix:

Though, multiplication is nicer

sage: (p*p).dict()
{0: 1, 2: 2, 4: 1}

What do you think about opening a ticket about enforcing non-zero values in the dictionary?

I am in favor of opening a ticket to add this specification, though I am not in favor of removing the check=False option that can be useful at some places to get (slightly) faster code.

I now understand that it only due to check=False. In the case check=True (which is the default) zeros are removed. So I would assume that there never are zero value in the dictionary. One only needs to specify it in the documentation.

You could have a look at the implementation of degree. It does not check if the values are zero or not. It just assumes that values are nonzero.

Vincent

[bruno]

Replying to vdelecroix:

Replying to bruno:

Replying to vdelecroix:

Dough! I was asking because for sparse matrices it is explicitely mentioned in the specification that no value can be zero. And this is very useful!

There is no such specification as for matrices. There is a method __normalize() to remove the zero coefficients, used by default in __init__, and used for instance in the multiplication (since there may be some cancellations).

matrix_generic_sparse.pyx, main docstring line 69-71

    A generic sparse matrix is
    represented using a dictionary whose keys are pairs of integers `(i,j)` and
    values in the base ring. The values of the dictionary must never be zero.

Replying to vdelecroix:

Though, multiplication is nicer

sage: (p*p).dict()
{0: 1, 2: 2, 4: 1}

What do you think about opening a ticket about enforcing non-zero values in the dictionary?

I am in favor of opening a ticket to add this specification, though I am not in favor of removing the check=False option that can be useful at some places to get (slightly) faster code.

I now understand that it only due to check=False. In the case check=True (which is the default) zeros are removed. So I would assume that there never are zero value in the dictionary. One only needs to specify it in the documentation.

Right. I guess we actually agree (though my previous comment was maybe unclear).

You could have a look at the implementation of degree. It does not check if the values are zero or not. It just assumes that values are nonzero.

True. And this is also the case of _repr:

sage: R.<x> = PolynomialRing(ZZ, sparse=True)
sage: R({0:1,1:0}, check=False)
 + 1

So we actually already assume that coefficients in the dictionary are nonzero. I'll add it in the specifications.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​536a7e5	Remove useless test
​c88f1b0	Add (tentative) specification

[vdelecroix]

Looks good to me.

[bruno]

Thanks for your comments!