Add a method `divides` to Polynomial

[bruno]

The generic method divides that can be found in src/sage/structure/element.pyx uses quo_rem (via %) to test if an element divides another one: return (x % self) == 0.

For polynomials, depending on the implementations, the method quo_rem may raise an error if the divisor is not monic (see #16649 for more on this).

This ticket aims at implementing a method divides for the class Polynomial, so that it catches the errors in quo_rem to return False when it is needed.

Example of a problematic behavior:

sage: R.<x> = PolynomialRing(ZZ, implementation="FLINT")
sage: p = 2*x + 1
sage: q = R.random_element(10)
sage: p.divides(q)
False
sage: R.<x> = PolynomialRing(ZZ, implementation="NTL")
sage: p = R(p)
sage: q = R(q)
sage: p.divides(q)
Traceback (most recent call last):
...
ArithmeticError: division not exact in Z[x] (consider coercing to Q[x] first)

[bruno]

New commits:

​98b7ca1	Rebased on 6.9beta1
​6d0da24	One last zero_element removed + correct one doctest
​dbfb518	Remove useless colons
​331c713	#19171: New method divides

[git]

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

​27040e0

#19171: New method divides

[git]

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

​9a189c7

#19171: New method divides

[bruno]

I added some noise with two forced pushes (due to some mistake I made), but there is really only one commit! Now the ticket passes the tests (as far as I can tell) and is very ready for review!

[vdelecroix]

You would better use coerce_binop from sage.structure.element rather than a manually handled coercion.

Why is this code specific to polynomials?

[git]

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

​d03a75e	#19171: New method divides
​ac43c02	19171: Use coerce_binop

[bruno]

Replying to vdelecroix:

You would better use coerce_binop from sage.structure.element rather than a manually handled coercion.

Done.

Why is this code specific to polynomials?

I tried to explain the reasons in the description: For polynomials over (say) ZZ, the euclidean division is not well defined. The problem occurs when the leading coefficient is not invertible. There are several ways to deal with this problem, and different implementations in Sage use different conventions (see #16649 for pointers on this). As a result, testing divisibility by invoking self % p == 0 can lead to ArithmeticError for polynomials: But when this exception is raised, we know that self does not divide p.¹ That's why I decided to add this new method which only differs from the generic one by the fact it catches the exception.

I could have changed the generic method instead to catch an ArithmeticError and return False in such case, but though I am pretty confident that it makes sense for polynomials, I am much less confident concerning other rings.

¹ "We know" is maybe a too strong assumption. In #16649, I did changes that ensure that the assumption holds for all current implementations of polynomials in Sage. (This was not the case before.) Yet, future implementations may break this rule. This may imply that we shouldn't simply rely on quo_rem to test divisibility.

[vdelecroix]

Since one is always a unit the if self.is_one(): return True would better be inside the block except NotImplementedError:.

Please add a more exotic example like GF(4)['x']['y'].

[vdelecroix]

And why not using quo_rem instead of catching a potential ArithmeticError?

[bruno]

I am not sure I fully understand your latest comments:

Since one is always a unit the if self.is_one(): return True would better be inside the block except NotImplementedError:.

Do you mean adding a try: ... except NotImplementedError: construction? Because if I add this test say right below the test if self.is_unit(): return True, an exception will be raised in the cases where is_unit is not implemented but is_one is.

And why not using quo_rem instead of catching a potential ArithmeticError?

Using p % self is really the same as p.quo_rem(self)[1] so I do not see the difference.

[vdelecroix]

Replying to bruno:

I am not sure I fully understand your latest comments:

Since one is always a unit the if self.is_one(): return True would better be inside the block except NotImplementedError:.

Do you mean adding a try: ... except NotImplementedError: construction? Because if I add this test say right below the test if self.is_unit(): return True, an exception will be raised in the cases where is_unit is not implemented but is_one is.

I meant

try:
    if self.is_unit(): return True   # units divide everything
except NotImplementedError:
    if self.is_one(): return True    # if is_unit is not implemented

And why not using quo_rem instead of catching a potential ArithmeticError?

Using p % self is really the same as p.quo_rem(self)[1] so I do not see the difference.

Not exactly. NTL makes a distinction... and you should actually use pseudo_quo_rem when it is there:

sage: R.<x> = PolynomialRing(ZZ, implementation="NTL")
sage: p = 2*x + 1
sage: q = R.random_element(10)
sage: p.quo_rem(q)
Traceback (most recent call last):
...
ArithmeticError: division not exact in Z[x] (consider coercing to Q[x] first)
sage: p.pseudo_quo_rem(q)
(0, 2*x + 1)

[git]

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

​0906dfd

19171: Add more involved example + move a line

[bruno]

Replying to vdelecroix:

...

I meant

...

Done.

And why not using quo_rem instead of catching a potential ArithmeticError?

Using p % self is really the same as p.quo_rem(self)[1] so I do not see the difference.

Not exactly. NTL makes a distinction... and you should actually use pseudo_quo_rem when it is there:

sage: R.<x> = PolynomialRing(ZZ, implementation="NTL")
sage: p = 2*x + 1
sage: q = R.random_element(10)
sage: p.quo_rem(q)
Traceback (most recent call last):
...
ArithmeticError: division not exact in Z[x] (consider coercing to Q[x] first)
sage: p.pseudo_quo_rem(q)
(0, 2*x + 1)

Not done: I do not think pseudo_quo_rem is applicable, because of the following:

sage: R.<x> = ZZ[] # works the same with NTL
sage: p = x^2 - 1
sage: q = 2*x + 2
sage: p.pseudo_quo_rem(q)
(2*x - 2, 0)

So, the remainder in p.pseudo_quo_rem(q) can be zero even if q does not divide p.

P.S.: I won't be able to work on this ticket until mid-august.

[git]

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

​f85a7e0

19171: conflict resolution

[tscrim]

Needs rebasing.

[git]

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

​04b0032	Merge branch 'develop' into 19171_divides_poly
​a0a5bad	Merge branch 'u/bruno/t19171_divides_poly' of git://trac.sagemath.org/sage into 19171_divides_poly

[mmezzarobba]

Hi,

This looks reasonable to me; I would be willing to set the ticket to positive review. Vincent, do you agree?

[vdelecroix]

why return (p % self) == 0 instead of return (p % self).is_zero()?

[vdelecroix]

How can you be sure that the following

        except ArithmeticError:
            return False                     # if division is not exact

only catches non exact divisions?

[vdelecroix]

Could you add check on more exotic rings like Zmod(6)['x'], ZZ['x']['y'], GF(2)['x,y']['z']?

[vdelecroix]

Replying to vdelecroix:

How can you be sure that the following

        except ArithmeticError:
            return False                     # if division is not exact

only catches non exact divisions?

BTW, each except will cost a non trivial amount of time. Isn't there a fastest way to see if the division is exact? (eg calling quo_rem)

[mmezzarobba]

FWIW, I agree that using is_zero() would be better, but I don't think it matters enough for blocking the ticket. Considering that an ArithmeticError implies that there is no divisibility looks very reasonable to me, though testing divisibility of the leading coefficients would likely be faster.

[git]

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

​f350d87	Merge branch 'u/bruno/t19171_divides_poly' of git://trac.sagemath.org/sage into t/19171/divides
​bcc4534	19171: Add leading coeff test + change '== 0' to is_zero
​5bac984	19171: Add tests

[bruno]

I've pushed my branch but some doctest fails. Actually I need some advice.

• The try... except ArithmeticError is needed for PolynomialRing(ZZ, 'x', implementation = 'NTL') since in this implementation quo_rem raises an exception for non-divisible polynomials, even when the leading coefficients divide each other. See 16649 for details.

• Contrary to what my (broken) commit message might indicate, I did not add a leading coefficient test that would fasten the computation. The problem comes from the Zmod(6)['x'] example: In Zmod(6), computing the remainder raise exceptions. (Test for instance Zmod(6)(1) % Zmod(6)(4).)

• As mentioned before, there is still a failing doctest, also due to Zmod(6)['x']. We face the same kind problem as before with the 'NTL' implementation of Z[x]: Computing quo_rem in Zmod(6)['x'] is only possible if the leading coefficient of the divisor is a unit (an exception is raised otherwise). The exception is not caught since contrary to the case of ZZ['x'], the raised exception is a ValueError.

I am not sure to see what the best solutions are to make it work. Some possibilities:

• We make this ticket depend on another one that unifies all the quo_rem implementations for polynomials over non-field, so that we can base this function on the common behavior of these methods.¹
• We make this ticket depend on another one with more limited goals, such as unifying the errors to ArithmeticError.
• We add some specific code, to catch the ValueError.

What do you think?

The first solution seem to me unreachable given the endless discussions about the behavior of quo_rem (links in 16649). The second one is maybe more doable. The third one is probably the faster to close this ticket but I am not sure whether there is a good way to implement it.

¹ More generally, there should exist a ticket for unifying the diverse behaviors of polynomial rings, but 1. it is clearly off-topic this ticket and 2. given the difficulty with unifying "content" in another ticket, I am not at all ready to go for this!

[vdelecroix]

An ArithmeticError in % does not imply not divisible as shown with

sage: R.<x> = Zmod(6)[]
sage: S.<z> = R[]
sage: p = x*z + 1
sage: q = z + 1
sage: (p*q) % p
Traceback (most recent call last):
...
ArithmeticError: Division non exact (consider coercing to polynomials over the fraction field)

The reason being that

sage: Z6 = Zmod(6)
sage: a = Z6(4)
sage: b = Z6(2)
sage: b.divides(a)
True
sage: a / b
Traceback (most recent call last):
...
ZeroDivisionError: Inverse does not exist.

From there three possibilities:

• "fix" divisions in integer mod ring Zmod so that a / b above does work. Though in this situation there are several possible answers as 4 = 2 x 5 = 2 x 2 mod 6.
• restrict divides on polynomial whose base rings are integral domains (so that division is well defined: quotient is unique when it exists).

[vdelecroix]

A fix is provided in #25277 for Zmod ring.

[bruno]

The difficulties with rings such as Zmod(6) let me think that we should indeed restrict to polynomials over integral domains. The problem being that I actually do not know for instance how to test divisibility over Zmod(6).

If we decide this, the problem with PolynomialRing(ZZ, 'x', implementation="NTL") remains since quo_rem cannot be used there without catching an exception, which is not very appealing... Should we convert first the polynomial to polynomials over the fraction field of the base ring?

[vdelecroix]

Replying to bruno:

The difficulties with rings such as Zmod(6) let me think that we should indeed restrict to polynomials over integral domains. The problem being that I actually do not know for instance how to test divisibility over Zmod(6).

There is #25277 for that (it is merged in beta3 already). So you can test for division... but you can not divide yet.

If we decide this, the problem with PolynomialRing(ZZ, 'x', implementation="NTL") remains since quo_rem cannot be used there without catching an exception, which is not very appealing... Should we convert first the polynomial to polynomials over the fraction field of the base ring?

It is fine to catch an exception, but it should be clear what it means.

[bruno]

Replying to vdelecroix:

Replying to bruno:

The difficulties with rings such as Zmod(6) let me think that we should indeed restrict to polynomials over integral domains. The problem being that I actually do not know for instance how to test divisibility over Zmod(6).

There is #25277 for that (it is merged in beta3 already). So you can test for division... but you can not divide yet.

I meant testing divisibility of polynomials over Zmod(6). I am not sure that it can be reduced to simply testing divisibility of elements of Zmod(6), I think you need to perform some divisions.

[vdelecroix]

Replying to bruno:

Replying to vdelecroix:

Replying to bruno:

The difficulties with rings such as Zmod(6) let me think that we should indeed restrict to polynomials over integral domains. The problem being that I actually do not know for instance how to test divisibility over Zmod(6).

There is #25277 for that (it is merged in beta3 already). So you can test for division... but you can not divide yet.

I meant testing divisibility of polynomials over Zmod(6). I am not sure that it can be reduced to simply testing divisibility of elements of Zmod(6), I think you need to perform some divisions.

A far from efficient method consist in

1. testing divisibility of leading coefficients b_m | a_n
2. if yes, test all possible quotients q so that a_n = q b_m

But for that, you need a method that return all quotients. Hence my remark.

[bruno]

So do you agree with me that we should better abandon non-integral rings (at least for now)?

[vdelecroix]

Replying to bruno:

So do you agree with me that we should better abandon non-integral rings (at least for now)?

yes

[git]

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

​fbc129c	Merge branch 'u/bruno/t19171_divides_poly' of git://trac.sagemath.org/sage into ticket/19171
​213ccd8	Add degree test and lc test

[git]

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

​235735f

19171: only for integral domains

[bruno]

Is the proposed solution ok for you?

[vdelecroix]

You are not doing the right thing on non-integral domains. You should test is_integral_domain much earlier. With the ticket applied

sage: R.<x> = Zmod(6)[]
sage: (2*x +1).divides(3)
False
sage: (2*x + 1) * 3 == 3
True

(testing divisibility of leading coefficient is not even valid... I was wrong)

You need a line break after TESTS::

It should be documented that this method only works for integral domains. And the examples with Zmod(6) (raising error) should be in the EXAMPLES block rather than TESTS.

[git]

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

​94c0390

19171: Divides implemented only for integral domains

[bruno]

I've followed your advice, and changed my mind: My idea was to give an answer (for non-integral domains) as long as I am able to do so. The problem with the divisibility of leading coefficients made me adopt the opposite viewpoint: Divisibility is not implemented for non-integral domains, that's it. I am not sure that an implementation which only works in some uninteresting cases would be helpful to anyone!

[vdelecroix]

update milestone 8.3 -> 8.4

[vdelecroix]

ok