Opened 6 years ago

Closed 6 years ago

#16276 closed defect (fixed)

fixes PolynomialElement.mod()

Reported by: Bouillaguet Owned by:
Priority: major Milestone: sage-6.3
Component: basic arithmetic Keywords:
Cc: Merged in:
Authors: Charles Bouillaguet Reviewers: John Cremona, Jean-Pierre Flori, Peter Bruin
Report Upstream: N/A Work issues:
Branch: dca4cc5 (Commits) Commit: dca4cc51ea0d0a6b9f321de8e7815ec5602b0593
Dependencies: Stopgaps:

Description

Before :

           sage: R.<x> = PolynomialRing(ZZ)
            sage: f = x^3 + x + 1
            sage: f.mod(x + 1)
            x^3 + x + 1

After :

           sage: R.<x> = PolynomialRing(ZZ)
            sage: f = x^3 + x + 1
            sage: f.mod(x + 1)
            -1

Change History (29)

comment:1 Changed 6 years ago by Bouillaguet

  • Branch set to u/Bouillaguet/ticket/16276
  • Created changed from 05/02/14 06:17:22 to 05/02/14 06:17:22
  • Modified changed from 05/02/14 06:17:22 to 05/02/14 06:17:22

comment:2 Changed 6 years ago by Bouillaguet

  • Commit set to e48396c79625bb06a06e27552e2cce9ab8bb60e1
  • Status changed from new to needs_review

New commits:

e48396creplace broken default PolynomialElement.mod() with already-present reasonable code

comment:3 Changed 6 years ago by git

  • Commit changed from e48396c79625bb06a06e27552e2cce9ab8bb60e1 to 07221cbd32a822e7d4ba40762bcf843a10a6c715

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

07221cbfix corresponding doctest

comment:4 follow-up: Changed 6 years ago by cremona

Question: are we 100% certain that every instance of the class which has this method does have a valid quo_rem method? If not (and may even if so) it would be good to wrap the call to quo_rem in a try...except block with a failsafe default to just return self.

comment:5 Changed 6 years ago by Bouillaguet

No we are not. As you pointed out, in ZZ[X], quo_rem cannot work reliably in all cases. By design, it does "something" when b is not monic and does not divides a exactly (it tries to reduce "as much as possible").

Do you think we should also patch quo_rem to return an ArithmeticError exception in this case where the result is not well-defined?

comment:6 in reply to: ↑ 4 Changed 6 years ago by Bouillaguet

Replying to cremona:

it would be good to wrap the call to quo_rem in a try...except block with a failsafe default to just return self.

My opinion is that it is better to raise an exception that to return mathematical nonsense (if only for "pedagogical" reasons)

comment:7 Changed 6 years ago by cremona

No, I don't think that would be helpful. We agreed (I think) that for the reduction of X^100 mod 2*X+1 in ZZ[X] there was no better result than leaving it as X^100, so would we really want to raise an error when someone types

sage: x = polygen(ZZ)
sage: (x^100) % (2*x+1)

instead of just returning

x^100

as happens now? Even this

sage: x = polygen(QQ)
sage: (x^100) % (2*x+1)
1/1267650600228229401496703205376

while being quite valid mathematically, is slightly surprising...

comment:8 Changed 6 years ago by Bouillaguet

It turns out that quo_rem does not always have the same specification (which seems bad in my opinion). Examples :

sage: R.<X> = ZZ[]
sage: (2*X^2 + 3*X + 1).quo_rem(2*X + 1)
(X + 1, 0)

So, here, Sage "silently" performs a partial reduction.

sage: RR.<Y> = R[]
RR.random_element().quo_rem((3*X+1)*Y + X + 2)
---------------------------------------------------------------------------
Traceback (most recent call last)
...
ArithmeticError: Nonunit leading coefficient

And indeed, the "generic polynomial" quo_rem (as in src/sage/ring/polynomial/polynomial_element.pyx) advertises that an arithmetic exception shall be raised when a division by a non-monic polynomial is attempted).

comment:9 Changed 6 years ago by jpflori

As far as the current commit is concerned, I'd say a simple addition of

    mod = __mod__

would do the trick rather than copying the __mod__ function. Otherwise, as there is a quo_remfunction defined, I would say we can rely on that, and let it raise an exception if it wants. Note that you can already issue:

f % g

and trigger the call to quo_rem.

comment:10 Changed 6 years ago by git

  • Commit changed from 07221cbd32a822e7d4ba40762bcf843a10a6c715 to 7e902de0445ec47ce835f3cb202b361793c621f5

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

7e902deavoid code duplication

comment:11 Changed 6 years ago by Bouillaguet

@jpflori : I implemented your suggestion.

comment:12 Changed 6 years ago by jpflori

Great. I would be happy to close the ticket as is cos it's a first step toward something better. Completely fixing pow, powmod, mod, quo_rem and so one looks like it needs more work and thoughts. John: would you mind taking this route?

comment:13 Changed 6 years ago by vbraun_spam

  • Milestone changed from sage-6.2 to sage-6.3

comment:14 Changed 6 years ago by pbruin

Is mod = __mod__ the correct way to say that these two functions should behave identically? If C is a subclass of Polynomial implementing __mod__ but not mod, then C.mod will still be equal to Polynomial.__mod__, not to C.__mod__, if I'm not mistaken.

This is currently the case for Polynomial_dense_modn_ntl_zz, Polynomial_dense_modn_ntl_ZZ and Polynomial_rational_flint. To automatically make mod == __mod__ for these classes, I think Polynomial.mod would have to be defined as

def mod(self, other):
    return self.__mod__(other)

comment:15 Changed 6 years ago by jpflori

Thanks Peter, you might be right, my comment was just a quick remark about code duplication.

comment:16 Changed 6 years ago by jpflori

And I can confirm that just putting mod = __mod__ gives the behavior you mention.

comment:17 Changed 6 years ago by Bouillaguet

Duly noted. I have made the change. I can't see how to avoid doctest duplication though.

comment:18 Changed 6 years ago by git

  • Commit changed from 7e902de0445ec47ce835f3cb202b361793c621f5 to 0061d17ca5e04465c1189d62aa032f98a0c120cf

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

ff93a0aModify FLINT implementation of ZZ[X] so that quo_rem() match specification
d64e7a3Correcting corresponding doctests
0061d17too much code deduplication kills code semantics

comment:19 Changed 6 years ago by pbruin

I don't think the approach in commit ff93a0a (Modify FLINT implementation...) is the right one. While Z[x] is not a Euclidean domain, it is undoubtedly useful for it to have a quo_rem() method that does as much as it can. It is certainly more useful than raising an exception, in my opinion.

If I understand it correctly (I haven't read the code, only the doctests), the existing algorithm for quo_rem(f, g) effectively computes the Hermite normal form of the matrix whose columns (or rows, if you use the Sage convention) are the coefficient vectors of the polynomials f and x^i * g for 0 <= i <= deg f - deg g. One can do this more generally for A[x] with A a Euclidean domain. This is quite well-defined mathematically, and reduces to the Euclidean algorithm in two interesting cases: either (1) the polynomials are constant, or (2) the base ring is a field.

So I think it would be better to change the documentation of quo_rem() to clarify what it actually does, and to keep the code unchanged. Hopefully this can be done while also keeping mod and __mod__ identical to each other and to the rem part of quo_rem().

comment:20 follow-up: Changed 6 years ago by Bouillaguet

Hi pbruin,

I made a case for quo_rem() implementing euclidean division in https://groups.google.com/d/msg/sage-devel/UVrx5CG4qA0/WUz2RItJtF8J Nobody complained (and John approved).

What you ask for can be obtained (as you suggested) by other means.

Note that depending on the ring, and on the implementation, quo_rem does not do the same thing. For instance, the NTL implementation of ZZ[X] does what I ask for, but the FLINT implementation of ZZ[X] does what you requested. This situation calls for clarification.

Note that src/sage/structure/element.pyx asks that quo_rem() and mod() implement euclidean division, or fail if the result is not mathematically defined... I think that this is fine (if anything, it may help students realize that the result is not defined !).

comment:21 in reply to: ↑ 20 Changed 6 years ago by pbruin

Replying to Bouillaguet:

I made a case for quo_rem() implementing euclidean division in https://groups.google.com/d/msg/sage-devel/UVrx5CG4qA0/WUz2RItJtF8J Nobody complained (and John approved).

Sorry, I did not closely follow the sage-devel discussion and only realised I didn't like this approach after seeing your patch.

Actually I interpreted John's comment:7 as saying that he was against raising an error given that there is also a reasonable answer that we can return without raising an error. Maybe he can clarify his opinion?

What you ask for can be obtained (as you suggested) by other means.

What I'm most concerned about is that examples which used to work now raise an error. For example, I think the following examples (of the FLINT implementation) are mathematically perfectly valid:

sage: R.<x> = ZZ[]
sage: (x^100).quo_rem(2*x+1)
(0, x^100)
sage: (5*x+7).quo_rem(3)
(x + 2, 2*x + 1)

Is there a simple way to obtain these results after your patch?

Note that depending on the ring, and on the implementation, quo_rem does not do the same thing. For instance, the NTL implementation of ZZ[X] does what I ask for, but the FLINT implementation of ZZ[X] does what you requested. This situation calls for clarification.

Yes, I agree that we should come up with a consistent approach. I can see that Z[x] (in the FLINT implementation) and K[t][x] (for a field K) behave inconsistently and that this should be fixed.

However, I also think one should think twice before choosing a solution that fails in more cases, rather than a solution that returns a result in more cases.

Note that src/sage/structure/element.pyx asks that quo_rem() and mod() implement euclidean division, or fail if the result is not mathematically defined... I think that this is fine (if anything, it may help students realize that the result is not defined !).

I don't agree that the result is necessarily undefined; that depends on what definition you use. FLINT in fact uses a definition of quo_rem(f, g) which is valid for all f and g != 0 (the one I described), and with your patch we would simply refuse to return this answer in many cases. (Note that the current version of the patch does not seem to conform to the documentation either: you check for the leading coefficient being 1 instead of any unit, and you don't raise an error if the remainder is 0.)

I think that we should interpret the documentation as what it is, namely documentation (descriptive), not a specification (prescriptive). In particular, it may be clarified/amended if necessary, and it should not force methods in derived classes to raise errors in all cases where the original method would raise an error. In this particular case it should allow (not necessarily force) quo_rem() to use a more general algorithm if the base ring is a Euclidean domain that is not a field, e.g. if the polynomial ring is Z[x] (base ring Z) or K[t][x] (base ring K[t]).

comment:22 Changed 6 years ago by cremona

I'm not sure which John is which,a nd I cannot remember exactly what I said before. However: please let us not raise errors in quo_rem unless there is something really invalid in the input such as quotient by 0! If there is nothing else available having q=0 and r = dividend is surely better. And as Peter says, we swrite the documentation afterwards.

comment:23 Changed 6 years ago by Bouillaguet

Hi all,

In fact, most of my rant comes of out of the following: the name "quo_rem" strongly suggests euclidean division (at least for me). In particular, the QUOtient bit. Because of how the system works, the quo_rem() method is also available on non-euclidean rings, where euclidean division does not make sense (but where exact division sometimes make sense, e.g. over integral domains).

I thought that the "remainder" (which makes sense mathematically as everyone agrees) could be obtained by :

sage: R.<x> = ZZ[]
sage: Ideal(2*x+1).reduce(x^100)
x^100

(and given the remainder, the "quotient" part could be obtained by exact division).

But it turns out that this is not always the case:

sage: Ideal(3).reduce(5*x+7)
Traceback (most recent call last)
...
TypeError: not a constant polynomial

sage: Ideal(R(3)).reduce(5*x+7)
5*x + 7  # <--- not what you want, and probably incorrect?

This is yet another inconsistency.

Lastly, let me mention that the NTL implementation of ZZ[X] already does what I thought was right (git blame tells dmharvey wrote it):

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

This is actually where I copied-pasted the error message (and some of the code I think).

Anyway, I don't want to put up a fight about this. I don't care. I will thus retract the modifications to FLINT's polynomial (which will remain inconsistent with NTL). The fix to the "mod" method should still be accepted by everyone though.

comment:24 Changed 6 years ago by jpflori

I also suggest to close this one to get mod fixed and open a follow up ticket to unify quo_rem behavior.

comment:25 Changed 6 years ago by git

  • Commit changed from 0061d17ca5e04465c1189d62aa032f98a0c120cf to dca4cc51ea0d0a6b9f321de8e7815ec5602b0593

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

dca4cc5revert non-consensual modification to quo_rem() for FLINT's implementation of ZZ[X]

comment:26 Changed 6 years ago by Bouillaguet

That should do it

comment:27 Changed 6 years ago by pbruin

  • Reviewers set to John Cremona, Jean-Pierre Flori, Peter Bruin

Thanks! I think we all agree that this neatly solves the original problem. The patchbot doesn't seem to have tested the current state of affairs yet; now running doctests.

comment:28 Changed 6 years ago by pbruin

  • Status changed from needs_review to positive_review

comment:29 Changed 6 years ago by vbraun

  • Branch changed from u/Bouillaguet/ticket/16276 to dca4cc51ea0d0a6b9f321de8e7815ec5602b0593
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.