Simplify __getitem__ for polynomials and deprecate slicing

[pbruin]

For polynomials, we implement a generic __getitem__() method, which depends on a new Cython method get_unsafe() to get single coefficients. Overall, this greatly simplifies the code. Except for truncation pol[:n], slicing for polynomials is deprecated, since there is no clear use case and it is not mathematically well-defined what the answer should be (see also comment 39 and further below).

This fix allows us to simplify the method PowerSeries_poly.__getitem__(). Also, the patch corrects some (previously wrong) slice bounds in LaurentPolynomial_univariate.__getitem__().

[pbruin]

Fixing this for all polynomials is somewhat involved because every polynomial class defines its own __getitem__() method. The following 14 methods are relevant (line numbers for SageMath 6.8.rc0):

laurent_polynomial.pyx:321:    def __getitem__(self, i):
padics/polynomial_padic_capped_relative_dense.py:381:    def __getitem__(self, n):
polynomial_element_generic.py:300:    def __getitem__(self,n):
polynomial_element.pyx:7608:    def __getitem__(self, n):
polynomial_gf2x.pyx:64:    def __getitem__(self, i):
polynomial_integer_dense_flint.pyx:415:    def __getitem__(self, n):
polynomial_integer_dense_ntl.pyx:284:    def __getitem__(self, n):
polynomial_modn_dense_ntl.pyx:180:    def __getitem__(self, n):
polynomial_modn_dense_ntl.pyx:658:    def __getitem__(self, n):
polynomial_modn_dense_ntl.pyx:1223:    def __getitem__(self, n):
polynomial_rational_flint.pyx:379:    def __getitem__(self, n):
polynomial_real_mpfr_dense.pyx:183:    def __getitem__(self, ix):
polynomial_zmod_flint.pyx:231:    def __getitem__(self, i):
polynomial_zz_pex.pyx:151:    def __getitem__(self,i):

[jdemeyer]

For matrices and vectors, this is done using get_unsafe() which is the implementation of a[i] for one index i. Then the slicing is done by a generic __getitem__ function.

[pbruin]

Replying to jdemeyer:

For matrices and vectors, this is done using get_unsafe() which is the implementation of a[i] for one index i. Then the slicing is done by a generic __getitem__ function.

I implemented this suggestion, now testing.

[jdemeyer]

I might review this once the dependency is merged.

[jdemeyer]

Just a remark: did you have a look at ​PySlice_GetIndicesEx? That might be used to avoid slice arithmetic like

+            start, stop, step = n.start, n.stop, n.step
+            if step is None:
+                step = 1
+            if start is None:
+                start = 0
+            elif start < 0:
+                start %= step
+            if stop is None or stop > d:
+                stop = d

[pbruin]

Replying to jdemeyer:

Just a remark: did you have a look at ​PySlice_GetIndicesEx?

I tried this, but unfortunately it does not return the correct indices for polynomials. For example, it wraps around negative start values (as in the Python convention that L[-1] with L a list returns the last item in L); this is incompatible with Sage's convention that P[n] returns 0 if P is a polynomial and n < 0.

[bruno]

Replying to pbruin:

Replying to jdemeyer:

Just a remark: did you have a look at ​PySlice_GetIndicesEx?

I tried this, but unfortunately it does not return the correct indices for polynomials. For example, it wraps around negative start values (as in the Python convention that L[-1] with L a list returns the last item in L); this is incompatible with Sage's convention that P[n] returns 0 if P is a polynomial and n < 0.

I think this behavior should be documented: I'd add a Warning block to the documentation to specify that it doesn't (and cannot!) respect Python convention. This may be obvious for Laurent polynomials, but it is not clear for polynomials.

[git]

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

​b12ea0b	Trac 18940: document behaviour of __getitem__ for polynomials vs. lists
​e70430a	Trac 18940: forbid step < 0 in Polynomial.__getitem__()

[git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

​6842d86

Trac 18940: forbid step < 0 in Polynomial.__getitem__()

[pbruin]

Replying to bruno:

Replying to pbruin:

Replying to jdemeyer:

Just a remark: did you have a look at ​PySlice_GetIndicesEx?

I tried this, but unfortunately it does not return the correct indices for polynomials. For example, it wraps around negative start values (as in the Python convention that L[-1] with L a list returns the last item in L); this is incompatible with Sage's convention that P[n] returns 0 if P is a polynomial and n < 0.

I think this behavior should be documented: I'd add a Warning block to the documentation to specify that it doesn't (and cannot!) respect Python convention. This may be obvious for Laurent polynomials, but it is not clear for polynomials.

OK, I added a warning and an example for this.

Also, in the meantime I realised that step <= 0 is (1) not necessarily well-defined, and (2) may cause Sage to crash with the current implementation. I added a few lines to detect this case and raise an IndexError.

[git]

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

​1c56779

Merge branch 'develop' into ticket/18940-getitem_step

[jdemeyer]

For some classes, you still reimplement __getitem__ for slices. Why is that? Why cannot you always use the generic implementation for slicing?

[jdemeyer]

Why not allow step < 0? Actually, I think that step = -1 might be the most mathematically useful value of step.

[jdemeyer]

Are you sure that the current behaviour is really what we want?

sage: pol = PolynomialRing(QQ, 'x')(range(10))
sage: pol
9*x^9 + 8*x^8 + 7*x^7 + 6*x^6 + 5*x^5 + 4*x^4 + 3*x^3 + 2*x^2 + x
sage: pol[2:]
9*x^9 + 8*x^8 + 7*x^7 + 6*x^6 + 5*x^5 + 4*x^4 + 3*x^3 + 2*x^2

I would have expected

sage: pol[2:]
9*x^7 + 8*x^6 + 7*x^5 + 6*x^4 + 5*x^3 + 4*x^2 + 3*x + 2

I know that this branch just extends the current behaviour. However, I think we need to step back and define mathematically what should be returned when asking for pol[a:b:c]. The current definition is not very useful.

I think that list(pol[a:b:c]) should be the same as list(pol)[a:b:c], except for different rules when the range is outside of 0 ... pol.degree().

[jdemeyer]

Peter, what's your use-case for pol[a:b:c], i.e. why did you create this ticket?

[jdemeyer]

For the record, the current definition of pol[a:b] was introduced with

commit ff63a1c2b5798a777dcef6a4366a43dc60a6bebd
Author: William Stein <wstein@gmail.com>
Date:   Fri Jan 26 02:42:10 2007 -0800

    A huge amount -- starting with writing tons of doctests for power series, pyrex integer_ring, etc.
    
        * changed before of polynomial and power series slicing to return objects of the same type,
          which is much more consistent.
        * optimized integer constructor
        * improved random number generation (e.g., ZZ.random_element()) ten times faster and more
          flexible.
        * separated out pmem_malloc

[tscrim]

We currently have

sage: R.<x> = QQ[]
sage: p = R.random_element((-1, 10)); p
-2/3*x^10 + 3/2*x^9 - 1/2*x^8 + 10*x^7 + x^6 + 3/2*x^5 + 1/3*x^4 + 5/2*x^3 - x^2 + 29*x - 1
sage: list(p)
[-1, 29, -1, 5/2, 1/3, 3/2, 1, 10, -1/2, 3/2, -2/3]

So list does go low to high and the current behavior is consistent with converting to a list.

+1 for supporting negative steps.

[jdemeyer]

Replying to tscrim:

the current behavior is consistent with converting to a list.

What do you mean? I think that these two outputs should be equal:

sage: pol = PolynomialRing(QQ, 'x')(range(10))
sage: list(pol[5:])
[0, 0, 0, 0, 0, 5, 6, 7, 8, 9]
sage: list(pol)[5:]
[5, 6, 7, 8, 9]

[tscrim]

Ah, I see what you mean. I'm sorry, I didn't read your example closely enough and though you were truncating off the 2 leading terms.

However, list and slice already don't commute because of the python wrap around, so I'm not so sure on how much we should enforce the exact same behavior (granted, an error gets raised in the polynomial case). I personally would find that the degree was lowered more surprising if we still wanted slices to return polynomials as opposed to a list of the coefficients. Yet I don't really have much invested in this.

[jdemeyer]

Replying to tscrim:

I personally would find that the degree was lowered more surprising

So what's in your opinion the correct answer for

sage: pol = PolynomialRing(QQ, 'x')(range(10))
sage: list(pol[2::2])

If we cannot agree on a "natural" definition for pol[a:b:c], perhaps the best compromise is simply

if isinstance(i, slice):
    raise NotImplementedError("polynomial slicing is not defined")

[tscrim]

Replying to jdemeyer:

Replying to tscrim:

I personally would find that the degree was lowered more surprising

So what's in your opinion the correct answer for

sage: pol = PolynomialRing(QQ, 'x')(range(10))
sage: list(pol[2::2])

That depends on if we want a slice of a polynomial returns a polynomial or a simple list of the coefficients. I am (very) slightly in favor of returning a polynomial, so with that [0, 2, 4, 6, 8].

If we cannot agree on a "natural" definition for pol[a:b:c], perhaps the best compromise is simply

if isinstance(i, slice):
    raise NotImplementedError("polynomial slicing is not defined")

This would be okay with me as well, but I don't think my opinion here should weigh very much here because I am mostly ambivalent about this behavior since I don't have a use for it.

[jdemeyer]

Replying to tscrim:

That depends on if we want a slice of a polynomial returns a polynomial or a simple list of the coefficients. I am (very) slightly in favor of returning a polynomial, so with that [0, 2, 4, 6, 8].

Sorry, I asked a bad question so I cannot interpret your answer. Please try again with this question:

sage: pol = PolynomialRing(QQ, 'x')(range(1,10))
sage: list(pol[3::2])

[tscrim]

Replying to jdemeyer:

Replying to tscrim:

That depends on if we want a slice of a polynomial returns a polynomial or a simple list of the coefficients. I am (very) slightly in favor of returning a polynomial, so with that [0, 2, 4, 6, 8].

Sorry, I asked a bad question so I cannot interpret your answer. Please try again with this question:

sage: pol = PolynomialRing(QQ, 'x')(range(1,10))
sage: list(pol[3::2])

I would want this:

sage: list(pol[3::2])
[0, 0, 0, 4, 0, 6, 0, 8]

Actually, my answer above was wrong and shouw be [0, 0, 2, 0, 4, 0, 6, 0, 8]. [Edit - Some changes due to me not thinking things completely through]

However, when trying to think about how this would work for Laurent polynomials, we currently have this behavior for list:

sage: R.<x> = LaurentPolynomialRing(QQ)
sage: p = ~x + 2 + 3*x + 4*x^2
sage: list(p)
[1, 2, 3, 4]
sage: p = 3*x + 4*x^2
sage: list(p)
[3, 4]

In some ways I would say the iterator for polynomials should be changed to agree with that of Laurent polynomials, in that it starts from the lowest degree. Although I understand there are many (good) reasons why we do not want that behavior. I feel that at some point we just have to choose the least bad of all options and then put up warning messages about our choices.

[jdemeyer]

Replying to tscrim:

I would want this:

sage: list(pol[3::2])
[0, 0, 0, 4, 0, 6, 0, 8]

Interesting. This is indeed what this branch does, but it doesn't correspond with my intuition about what it should do (I would expect [4, 6, 8]).

[pbruin]

Replying to jdemeyer:

Peter, what's your use-case for pol[a:b:c], i.e. why did you create this ticket?

Because of comments 22, 25 and 26 on ticket #15601.

[pbruin]

Replying to jdemeyer:

Are you sure that the current behaviour is really what we want?

sage: pol = PolynomialRing(QQ, 'x')(range(10))
sage: pol
9*x^9 + 8*x^8 + 7*x^7 + 6*x^6 + 5*x^5 + 4*x^4 + 3*x^3 + 2*x^2 + x
sage: pol[2:]
9*x^9 + 8*x^8 + 7*x^7 + 6*x^6 + 5*x^5 + 4*x^4 + 3*x^3 + 2*x^2

Yes. The fact that polynomials are usually written with decreasing exponents is just a historical convention. I think it is absolutely reasonable that pol[2:] should return the monomials of degree 2 and higher, whether as a polynomial or as a list. And in fact

sage: list(pol)
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]

so doing list(pol)[2:] also extracts the coefficients of the monomials of degree at least 2.

Similarly, in PARI, my experience is that Vecrev (which does the equivalent of list above) is more useful than Vec to extract the list of coefficients of a polynomial.

[pbruin]

Replying to jdemeyer:

Replying to tscrim:

the current behavior is consistent with converting to a list.

What do you mean? I think that these two outputs should be equal:

sage: pol = PolynomialRing(QQ, 'x')(range(10))
sage: list(pol[5:])
[0, 0, 0, 0, 0, 5, 6, 7, 8, 9]
sage: list(pol)[5:]
[5, 6, 7, 8, 9]

What would you want the output of pol[5:] itself to be?

[pbruin]

Replying to jdemeyer:

Why not allow step < 0? Actually, I think that step = -1 might be the most mathematically useful value of step.

We already have the reverse() method for that. I can think of two possible results for pol[a:b:-1], one of which would be the reverse() of the other. I am afraid allowing step = -1 would only create confusion...

[jhpalmieri]

Mathematically, I can certainly imagine wanting the part of the polynomial between degrees 2 and 5: I might want this:

sage: f
-10*x^6 + 1/2*x^5 - x^4 + 7*x^3 - 2*x^2 + x - 3
sage: f[2:6]
1/2*x^5 - x^4 + 7*x^3 - 2*x^2

I don't know if I would ever want f[2:6] to return the polynomial with same coefficients but translated down in degree (i.e., 1/2*x^3 - x^2 + 7*x - 2). Computationally, though, if I have that second result, I can easily get back to the first one by multiplying by x^2. If I have the first, I can get to the second by dividing, but that puts the result in the fraction field, which is probably not what I want. Alternatively, I can also get to the second by creating a polynomial from list(f)[2:6], which is not as intuitive, but not too hard.

Surely different mathematicians will have different points of view on this.

I think this is the issue for me: do we want

sage: list(pol[5:]) == list(pol)[5:]
True

or do we want

sage: pol[:5] + pol[5:] == pol
True

We can't have both, right? I prefer the second. That is, I tend to think of polynomials as graded objects and slicing is returning some of their homogeneous components, and it seems most natural to do this without regrading the result: if a term is in degree 3 in the original polynomial, it should still be in degree 3 after slicing. I would also like this:

sage: pol[0:10:2] + pol[1:10:2] == pol[0:10]
True

Adding the even terms to the odd terms should yield the original polynomial.

A related question: should f[2] return -2 or -2 * x^2?

---

When it comes down to it, I'm not sure I care which option is implemented as __getitem__, but perhaps both should be implemented as methods which are visible via tab-completion, or as a single method with a keyword which determines which option to use.

[jdemeyer]

Replying to pbruin:

Replying to jdemeyer:

Peter, what's your use-case for pol[a:b:c], i.e. why did you create this ticket?

Because of comments 22, 25 and 26 on ticket #15601.

You're just justifying a dubious use-case with a different dubious use-case.

[jdemeyer]

Replying to pbruin:

Similarly, in PARI, my experience is that Vecrev (which does the equivalent of list above) is more useful than Vec to extract the list of coefficients of a polynomial.

I think you are totally misunderstanding me. This is not at all about in which order the monomials are extracted. Obviously, PARI's Vecrev() is the right thing, I'm not arguing against that.

It's about

9*x^9 + 8*x^8 + 7*x^7 + 6*x^6 + 5*x^5 + 4*x^4 + 3*x^3 + 2*x^2

versus

9*x^7 + 8*x^6 + 7*x^5 + 6*x^4 + 5*x^3 + 4*x^2 + 3*x + 2

[jdemeyer]

Replying to pbruin:

What would you want the output of pol[5:] itself to be?

sage: pol = PolynomialRing(QQ, 'x')(range(10))
sage: pol[5:]
9*x^4 + 8*x^3 + 7*x^2 + 6*x + 5

(although, I don't really want it to be this, but this is what I would expect)

[jdemeyer]

What I really want is probably

sage: pol = PolynomialRing(QQ, 'x')(range(10))
sage: pol[5:]
DeprecationWarning: slicing polynomials is deprecated...

[jdemeyer]

In the Sage library, there are use cases for pol[:n] (i.e. truncation) and that's also the only case (I think) where both definitions agree.

[jhpalmieri]

So deprecate slicing unless it's a truncation?

[jdemeyer]

One more thought. Let's call this degree-keeping slicing:

sage: pol = PolynomialRing(QQ, 'x')(range(10))
sage: pol[3::2]
9*x^9 + 7*x^7 + 5*x^5 + 3*x^3

and let's call this degree-reducing slicing:

sage: pol = PolynomialRing(QQ, 'x')(range(10))
sage: pol[3::2]
9*x^3 + 7*x^2 + 5*x + 3

The implementation of certain divide-and-conquer algorithms (e.g. Karatsuba multiplication and FFT) can be expressed very easily using degree-reducing slicing.

I'm not saying that this is a particularly strong argument for degree-reducing slicing, since most likely you'd implement those algorithms using lists instead of polynomials anyway.

[jdemeyer]

A new idea, let's call it double slicing. The argument to __getitem__ can be a pair of slices i, j, where the slice i is used to define the coefficients and the second slice j is used to define the exponents. It would allow for example:

sage: pol = PolynomialRing(QQ, 'x')(range(10))
sage: pol[:, ::2]
9*x^18 + 8*x^16 + 7*x^14 + 6*x^12 + 5*x^10 + 4*x^8 + 3*x^6 + 2*x^4 + x^2
sage: pol[:5, 100:]
4*x^104 + 3*x^103 + 2*x^102 + x^101

We could define shorthands pol[i] equivalent to pol[i, i] for degree-keeping slicing and pol[i,] equivalent to pol[i, :] for degree-reducing slicing.

This would be very general and cover all possible use-cases. The question is: is it over-engineered? Do we really need this?

[jdemeyer]

I can confirm that the Sage library doesn't use any polynomial slicing except for truncation pol[:n]. After adding

if start or step is not None:
    raise NotImplementedError

to the __getitem__ methods, the only doctest failures were with slicing in the actual doctest example.

[pbruin]

Replying to jdemeyer:

Replying to pbruin: I think you are totally misunderstanding me. This is not at all about in which order the monomials are extracted. Obviously, PARI's Vecrev() is the right thing, I'm not arguing against that.

It's about

9*x^9 + 8*x^8 + 7*x^7 + 6*x^6 + 5*x^5 + 4*x^4 + 3*x^3 + 2*x^2

versus

9*x^7 + 8*x^6 + 7*x^5 + 6*x^4 + 5*x^3 + 4*x^2 + 3*x + 2

OK, I see. As you can guess from how I interpreted your proposal, for me this is not the expected answer at all.

I would personally be quite happy with just deprecating the start and stop parameters and allowing polynomial slicing only for pol[:n], i.e. truncation. I think users can easily implement other forms of slicing via lists.

[tscrim]

I think the double-slicing is probably something we don't really need. So I'm also for deprecating the (non-zero?) start and step conditions.

[jdemeyer]

If we deprecate general slicing, shouldn't we deprecate all kinds of slicing, including pol[:n] in favour of the truncate() method?

[jdemeyer]

Truncation is actually implemented using slicing, so we cannot easily deprecate that.

I am working on a patch.

[git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

​9478486

Deprecate slicing with start index

[jdemeyer]

Sorry, forgot src/sage/rings/polynomial/padics/polynomial_padic_capped_relative_dense.py

[git]

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

​b559028	Use __new__ to create an empty RealNumber
​44cca40	Deprecate slicing with start index also in padics

[mmezzarobba]

FWIW, I think I have personal code that uses (“degree-keeping”) slicing of the form pol[n:]. I can port it, but I may not be alone...

[jdemeyer]

Cool, what's your use case? I'm curious because I couldn't think of any.

[mmezzarobba]

Replying to jdemeyer:

Cool, what's your use case? I'm curious because I couldn't think of any.

Splitting polynomials that actually represent power series expansions into an “approximation” and a “remainder“; removing the constant term of a polynomial; removing low-order terms that are know to be zero mathematically but are computed as intervals containing zero from a residual Φ(pol) where Φ is an operator and pol is a polynomial close to a zero of Φ.

[jdemeyer]

@mmezzarobba: what is your opinion on 40?

[mmezzarobba]

Replying to jdemeyer:

@mmezzarobba: what is your opinion on 40?

I don't think I'd need it—any of the two basic forms of slicing (either “degree-keeping” or “degree-reducing”) is enough for me. Regardless whether “double slicing” is implemented, I'd prefer pol[i:] to keep working (with one of the two basic meanings, I don't really care which, though I find the “degree-keeping” version more intuitive with polynomials). It wouldn't be a big effort to change my code to work without slicing in any case.

[pbruin]

Fixed a merge conflict.

[rws]

This is one of those enhancements I like to do myself. It is also thoroughly done. Since patchbot is happy too, I feel confident to set it positive.

[roed]

Replying to jdemeyer:

I can confirm that the Sage library doesn't use any polynomial slicing except for truncation pol[:n]. After adding

if start or step is not None:
    raise NotImplementedError

to the __getitem__ methods, the only doctest failures were with slicing in the actual doctest example.

For a use case where deprecating pol[n:] is annoying, if you have an Eisenstein polynomial pol and want to compute p/x^n modulo pol, it's natural to express it as -(pol[n:] >> n) / (pol[:n] // p), where the division is taken modulo pol. Obviously, this is doable with lists.

[jdemeyer]

Maybe a different solution would be to re-introduce the slicing but with a different interface (not __getitem__ but an explicit method).