modulo is broken in some (fairly simple) rings

[Charles Bouillaguet]

Hi all,

I just noticed that we carefully check that Sage computes something completely wrong. From sage/structure/element.pyx, line 2052 :

--------------------------
When little is implemented about a given ring, then mod may
return simply return `f`. For example, reduction is not
implemented for `\ZZ[x]` yet. (TODO!)

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

The problem is that mod() does the following :

from sage.rings.ideal import is_Ideal
if not is_Ideal(I) or not I.ring() is self._parent:
I = self._parent.ideal(I)
#raise TypeError, "I = %s must be an ideal in %s"%(I, self.parent())
return I.reduce(self)

So, mod relies on Ideal.reduce, which by default does nothing and returns its argument. As such, mod is badly broken.

What seems reasonable is to make Ideal.reduce() raise a NotImplementedError by default instead of doing nothing, by this breaks a lot of other stuff (because Ideal.reduce() is called a lot).

What do you think is right ?

*) Go with NotImplementedError and fix 30+ doctests?
*) create a StopGap that warn the user that mod may be complete rubbish?

Cheers,
---
Charles Bouillaguet
http://www.lifl.fr/~bouillaguet/

[Charles Bouillaguet]

On 30 Apr 2014, at 08:28, Charles Bouillaguet <charles.b...@gmail.com> wrote:

> Hi all,
>
> I just noticed that we carefully check that Sage computes something completely wrong. From sage/structure/element.pyx, line 2052 :
>
> --------------------------
> When little is implemented about a given ring, then mod may
> return simply return `f`. For example, reduction is not
> implemented for `\ZZ[x]` yet. (TODO!)
>
> sage: R.<x> = PolynomialRing(ZZ)
> sage: f = x^3 + x + 1
> sage: f.mod(x + 1)
> x^3 + x + 1
> ----------------------

Well, the good thing is that

sage: f % (x+1)
-1

DOES work... (but still, f advertises a broken "mod" method).

[Jean-Pierre Flori]

Wouldn't modify the mod method for polynomial rings over ZZ be enough?

In a similar way, the pow (python?) function accepts an optional "modulus" third args, but Sage just ignores it over a bunch of rings.
For example, over QQ, pow(x,100,x+1) yields x^100.
For non euclidean rings that might be a sensible default choice but for euclidean rings or even just univariate polynomial rings over fields we might want to modify that (maybe something for the category framework as Luca De Feo suggested).

[John Cremona]

Of course ZZ[X] is not Euclidean and it is not clear (to me) what a
good definition of reduction modulo (say) 2*X+1 should be. But
reduction modulo a monic polynomial is easy, and that special case
would be worth having even if one did nothing in the general case.

Quiz question: what do YOU think that the reduction of X^100 modulo
2X+1 should be?

John

>
> --
> You received this message because you are subscribed to the Google Groups
> "sage-devel" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sage-devel+...@googlegroups.com.
> To post to this group, send email to sage-...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sage-devel.
> For more options, visit https://groups.google.com/d/optout.

[Jean-Pierre Flori]

Good question...
We had a look around last week to gather info about what other libs do and IIRC, at least for flint, what happens is that start at top degree monomial you substract the largest multiple of the modulus so that the leading coefficient is positive (did not look at negative leading coeffs...).
In fact, that's also what Sage does when you do (x^100) % (x^2+1) (which does not call the same function as pow(x,100,x^2+1)).

Anyway in my case, where I wanted fast modular exponentiation, so reducing at each step was necessary to keep degrees low, such a case does not make sense as with such a definition the resulting polynomial has high degree anyway (and I don't really care about such a case...).

[John Cremona]

Sure. I think the underlying mathematical issue is that Z[X]/(2*X+1)
is not finitely-generated as an abelian group, (equivalently, -1/2 is
not an algebraic integer), so it is hopless to look for a finite set
of "residues" with the property that everything reduces to a Z-linear
combination of these. As a ring this is Z[1/2] so one choice is to
return the rational number obtained by substituting -1/2 for X. But I
think that if Sage were to return 2^-50 for X^100 % (2*X+1) we would
get complaints....

Sorry for the overkill, I have just been teaching a course in
commutative algebra!

[Charles Bouillaguet]

Anyway, this is now :

http://trac.sagemath.org/ticket/16276

waiting for a ---very simple--- review. (All test pass for me)
---
Charles

[Nils Bruin]

On Wednesday, April 30, 2014 2:13:20 AM UTC-7, Jean-Pierre Flori wrote:

Good question...
We had a look around last week to gather info about what other libs do and IIRC, at least for flint, what happens is that start at top degree monomial you substract the largest multiple of the modulus so that the leading coefficient is positive (did not look at negative leading coeffs...).
In fact, that's also what Sage does when you do (x^100) % (x^2+1) (which does not call the same function as pow(x,100,x^2+1)).

Anyway in my case, where I wanted fast modular exponentiation, so reducing at each step was necessary to keep degrees low, such a case does not make sense as with such a definition the resulting polynomial has high degree anyway (and I don't really care about such a case...).

There's a good theory for this and it's the same as for multivariate polynomial rings: you have to specify a "monomial order" (which should include the integers somehow!)  and then techniques for groebner bases apply . For univariate polynomials this should work out to what you describe: you look at the leading monomials, multiply by a power of x to make the degrees match, and do a euclidean division on the leading coefficients to decide which integer multiple you should describe. Repeat until you cannot reduce the degree any further.

(possibly material for John's follow-up course?)

[Charles Bouillaguet]

Dear all,

I would like to make a case for the quo_rem() method to actually implement euclidean division when it makes sense, or to raise an exception when the result is not defined. After all, this is what it is **supposed** to do, at least according to the docstring in sage/ring/polynomial/polynomial_element.pyx :

"""
@coerce_binop
def quo_rem(self, other):
"""
Returns the quotient and remainder of the Euclidean division of
``self`` and ``other``.

Raises ZerodivisionError if ``other`` is zero. Raises ArithmeticError if ``other`` has
a nonunit leading coefficient.
"""

I personally think that this specification is correct and decent.

Presently, depending on the domain, quo_rem does not always behave the same, and not according to the same specification. Examples :

In ZZ[X] or ZZ[X,Y] :
-------

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

---> Result is not defined, but quo_rem() does not complain. It tries to reduce "as much as possible".

in QQ[x][y] :
----------------

sage: P.<x> = QQ[]
sage: R.<y> = P[]
sage: f = R.random_element(10)
sage: g = y^5+R.random_element(5)
sage: q,r = f.quo_rem(g)
sage: f == q*g + r
True
sage: g = x*y^5
sage: f.quo_rem(g)
Traceback (most recent call last):
...
ArithmeticError: Nonunit leading coefficient

---> quo_rem() complains when result is not defined.

So, I suggest that :
a) we stick to the spec (what happens in QQ[x][y]).
b) we fix the cases like ZZ[X] that do not stick to the spec
c) we make sure, eventually by adding a "reduce" method or something like that, that the previous behavior is still available somehow.

I am willing to implement these changes. Do we agree on this?

[John Cremona]

+1 from me. Will this replace he current patch/branch at #16276?

John

>
> ---
> Charles Bouillaguet
> http://www.lifl.fr/~bouillaguet/
>
>
>