Force tail reduction in polynomial quotient ring

[rachelplayer]

I'd like to "remove squares" in some polynomials living in a polynomial ring over QQ, in 2 variables: x,y. I tried to implement this by modding out by the ideal (x^2 - x, y^2 - y). Depending on the ordering, the result of .mod() does not always output the polynomial I am looking for.

Without specifying an ordering, everything seems fine:

sage: R1.<x,y> = PolynomialRing(QQ, 2)
sage: I1 = R1.ideal(["x^2 - x", "y^2 - y"])
sage: R1("x^2 + y").mod(I1)
x + y
sage: R1("x + y^2").mod(I1)
x + y

However, when specifying the order lex the reduction of x + y^2 is not as expected:

sage: R2.<x,y> = PolynomialRing(QQ, 2, order="lex")
sage: I2 = R2.ideal(["x^2 - x", "y^2 - y"])
sage: R2("x^2 + y").mod(I2)
x + y
sage: R2("x + y^2").mod(I2)
x + y^2

This issue was reported in sage-support where it was pointed out that it is likely a bug in Singular, or in the Singular interface to Sage.

In particular, using the order lex works when implementation="generic" is also specified:

sage: R3.<x,y> = PolynomialRing(QQ, 2, order="lex", implementation="generic")
sage: I3 = R3.ideal(["x^2 - x", "y^2 - y"])
sage: R3("x^2 + y").mod(I3)
x + y
sage: R3("x + y^2").mod(I3)
x + y

For reference, I am using Sage version 8.6 on macOS Mojave 10.14.3.

PS. see also https://groups.google.com/d/msg/sage-devel/K49-V3BbWbg/pxuoehPvAAAJ

CC: @simon-king-jena @malb @mwageringel

Component: commutative algebra

Keywords: multivariate polynomial, quotient ring, singular

Author: Dima Pasechnik

Branch: 21e4f9a

Reviewer: Markus Wageringel

Issue created by migration from https://trac.sagemath.org/ticket/27508

[dimpase]

comment:1

Thanks, one would put "author" in only if intended to provide a patch soon.

[rachelplayer]

Changed author from Rachel Player to none

[rachelplayer]

comment:3

Thanks, in that case I have removed my name from "author".

[dimpase]

comment:4

Singular understand option "lp" for the lexicographic order.
See https://www.singular.uni-kl.de/Manual/4-0-3/sing_31.htm#SEC43

And if I define R2.<x,y> = PolynomialRing(QQ, 2, order="lp")
then it all works as expect. So the bug is in translating the Sage's option "lex" into Singular's one, should be easy to fix.

[dimpase]

comment:5

The following appears to fix this, although I don't really now why...

--- a/src/sage/libs/singular/ring.pyx
+++ b/src/sage/libs/singular/ring.pyx
@@ -50,7 +50,7 @@ from collections import defaultdict
 order_dict = {
     "dp": ringorder_dp,
     "Dp": ringorder_Dp,
-    "lp": ringorder_lp,
+    "lex": ringorder_lp,
     "rp": ringorder_rp,
     "ds": ringorder_ds,
     "Ds": ringorder_Ds,

"lp" order still works, no idea why.

[kwankyu]

comment:7

The file "src/sage/rings/polynomial/term_order.py" might be related.

[dimpase]

comment:8

This bug is already in Sage 8.1, that's as far back as I went.
Perhaps it was always there, I don't know.

Broken lex order should be a blocker, I think...

[jdemeyer]

comment:9

Replying to @dimpase:

This bug is already in Sage 8.1, that's as far back as I went.

If a bug is that old and nobody noticed before, then surely it's not serious enough to warrant blocker status? But I'll let the release manager judge that.

[nbruin]

comment:10

Replying to @dimpase:

"lp" order still works, no idea why.

What do you mean with "still works"? At the places where order_dict is used, a default value of ringorder_dp is used, so if lp is not a key in the dict anymore, I would expect that you end up with dp instead.

[dimpase]

comment:11

Replying to @nbruin:

Replying to @dimpase:

"lp" order still works, no idea why.

What do you mean with "still works"? At the places where order_dict is used, a default value of ringorder_dp is used, so if lp is not a key in the dict anymore, I would expect that you end up with dp instead.

No, an unknown value of order leads to error:

sage: R2.<x,y> = PolynomialRing(QQ, 2, order="baz")
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-5-4b6111e40869> in <module>()
----> 1 R2 = PolynomialRing(QQ, Integer(2), order="baz", names=('x', 'y',)); (x, y,) = R2._first_ngens(2)

/mnt/opt/Sage/sage-dev/local/lib/python2.7/site-packages/sage/rings/polynomial/polynomial_ring_constructor.pyc in PolynomialRing(base_ring, *args, **kwds)
    646 
    647     if multivariate or len(names) != 1:
--> 648         return _multi_variate(base_ring, names, **kwds)
    649     else:
    650         return _single_variate(base_ring, names, **kwds)

/mnt/opt/Sage/sage-dev/local/lib/python2.7/site-packages/sage/rings/polynomial/polynomial_ring_constructor.pyc in _multi_variate(base_ring, names, sparse, order, implementation)
    761     from sage.rings.polynomial.term_order import TermOrder
    762     n = len(names)
--> 763     order = TermOrder(order, n)
    764 
    765     # "implementation" must be last

/mnt/opt/Sage/sage-dev/local/lib/python2.7/site-packages/sage/rings/polynomial/term_order.pyc in __init__(self, name, n, force)
    729                     else: # simple order
    730                         if name not in print_name_mapping.keys() and name not in singular_name_mapping.values():
--> 731                             raise ValueError("unknown term order {!r}".format(name))
    732                         self._length = n
    733                         self._name = name

ValueError: unknown term order 'baz'

[dimpase]

comment:12

I suspect it is easy to fix, as it's some sort of translation issue between orders of Sage and orders of Singular, and people who wrote it are reading this ticket, I gather.

Please, please, do something, it is very serious, a broken Groebner basis for lex order is really a blocker.

[jdemeyer]

comment:13

If the bug is already in Sage 8.1, Sage 8.2, Sage 8.3, Sage 8.4, Sage 8.5 and Sage 8.6, would it really make such a big difference to have the same bug also in Sage 8.7?

But it doesn't matter what I think, the release manager will decide.

[dimpase]

comment:14

I also checked now that it is already in 7.4. This is a bug of the type "2+2=5", if you ask me, it's a shame to release anything with this kind of thing knowingly broken without either a fix or a fat warning somewhere...

[kwankyu]

comment:15

I chased this issue a little more. It seems that there is no problem in the translation. Perhaps it is a deeper issue.

Ultimately the normal form computation occurs in the code starting from line 283 in src/sage/libs/singular/groebner_strategy.pyx

Someone else who knows libSingular well, not me, is invited to look into it.

[dimpase]

comment:16

Replying to @kwankyu:

I chased this issue a little more. It seems that there is no problem in the translation. Perhaps it is a deeper issue.

This does not explain why after the change in comment:5
one can do

sage: R2.<x,y> = PolynomialRing(QQ, 2, order="lex")
sage: I2 = R2.ideal(["x^2 - x", "y^2 - y"])
sage: R2("x + y^2").mod(I2)
x + y
sage: R2._singular_()
polynomial ring, over a field, global ordering
// coefficients: QQ
// number of vars : 2
//        block   1 : ordering lp
//                  : names    x y
//        block   2 : ordering C
sage: R2.<x,y> = PolynomialRing(QQ, 2, order="lp")
sage: I2 = R2.ideal(["x^2 - x", "y^2 - y"])
sage: R2("x + y^2").mod(I2)
x + y
sage: R2._singular_()
polynomial ring, over a field, global ordering
// coefficients: QQ
// number of vars : 2
//        block   1 : ordering lp
//                  : names    x y
//        block   2 : ordering C

even though if order_dict is a (mathematical) function then "lp" should not have worked. (And indeed replacing lp by baz gives an error at once.

So you see above the order is lp, a.k.a. lex, all good. Well, without the comment:5 change:

sage: R2.<x,y> = PolynomialRing(QQ, 2, order="lex")
sage: R2._singular_() # all good, no?
polynomial ring, over a field, global ordering
// coefficients: QQ
// number of vars : 2
//        block   1 : ordering lp
//                  : names    x y
//        block   2 : ordering C
sage: I2 = R2.ideal(["x^2 - x", "y^2 - y"])
sage: R2("x + y^2").mod(I2) # OOPS!!!!
x + y^2

Ultimately the normal form computation occurs in the code starting from line 283 in src/sage/libs/singular/groebner_strategy.pyx

at this point the order is already set. From the conversation about it appears that nobody knows the real purpose order_dict.

Someone else who knows libSingular well, not me, is invited to look into it.

[dimpase]

comment:17

Replying to @nbruin:

Replying to @dimpase:

"lp" order still works, no idea why.

What do you mean with "still works"? At the places where order_dict is used, a default value of ringorder_dp is used, so if lp is not a key in the dict anymore, I would expect that you end up with dp instead.

please also see comment:16 about this - no, at no point here the "internal" Singular's order becomes dp.

[dimpase]

comment:18

More about the bug: it appears that only the 1st term is reduced:

sage: R2.<y,z,t,x> = PolynomialRing(QQ, 4, order="lex")
....: I2 = R2.ideal([x^2 + x, y^2 +2*y, z^2 - z, t^2 - t])
....: p = x^2 + y^2+z^2+t^2; p.mod(I2)
....: 
-2*y + z^2 + t^2 + x^2
sage: I2.reduce(p)
-2*y + z^2 + t^2 + x^2

(similarly with 3 variables...)

Now, reading code tells me that p.mod(I2) does I2.reduce(p)
which in turn tries the strat thing, which "works":

        try:
            strat = self._groebner_strategy()
            return strat.normal_form(f)
        except (TypeError, NotImplementedError, ValueError):
            pass

and gets the wrong answer. If it didn't work then it'd try to reduce using an explicit Groebner basis for I2, something that gives a correct answer:

sage: p.reduce(I2.groebner_basis())
-2*y + z + t - x

So indeed it's that strat thing. Note that if I specify order="lp" then calling _groebner_strategy() fails, and the result is (correctly) computed reducing with the Groebner basis.

[vbraun]

comment:19

We don't really promise anywhere to do tail reduction I think; Singular does the right thing by default but also not necessarily:

> ring r=0,(x,y),lp;
// ** redefining r (ring r=0,(x,y),lp;)
> ideal i=x2-x,y2-y;
> ideal gb=groebner(i);
> reduce(x+y2, gb, 0);     // Tail reduction (default)
x+y
> reduce(x+y2, gb, 1);     // No tail reduction
x+y2

While we should do tail reduction there is nothing mathematically wrong, so imho no blocker.

[dimpase]

comment:21

Replying to @vbraun:

We don't really promise anywhere to do tail reduction I think;

yes we do:

age: I2.reduce?
Signature:      I2.reduce(f)
Docstring:     
   Reduce an element modulo the reduced Groebner basis for this ideal.
   This returns 0 if and only if the element is in this ideal. In any
   case, this reduction is unique up to monomial orders.

and here we see different reductions depending upon the chosen implementation for the ring, nothing one can call "unique up to monomial orders".
(the default implementation gives one result, the generic one gives another result).

Perhaps Sage should allow for not doing the tail reduction, but it must not be the default in any case (it's not the default in Singular, the optional 3rd argument in reduce is by default 0 in Singular)

[nbruin]

comment:44

git blame ascribes this to the original check-in of the file groebner_strategy.pyx (by Martin Albrecht), so we've never done differently. Given that we have pretty authoritative advice now, I think it would be reasonable to just ascribe it to incomplete information at the time, and make the change now. If there was a good reason, we'll discover it ...

[dimpase]

comment:45

the advice is from 3rd author of groebner_strategy.pyx, so that's a bit worrying :-)

[malb]

comment:46

To add to the worry: I have no recollection why I did what I did.

[dimpase]

Changed author from Dima Pasechnik to none

[embray]

comment:48

As the Sage-8.8 release milestone is pending, we should delete the sage-8.8 milestone for tickets that are not actively being worked on or that still require significant work to move forward. If you feel that this ticket should be included in the next Sage release at the soonest please set its milestone to the next release milestone (sage-8.9).

[dimpase]

Branch: u/dimpase/singular/singu_nf

[dimpase]

comment:50

by popular demand, a minimal patch.

---

New commits:

21e4f9a

a minimal fix to ensure that 'lex' order always reduces

[dimpase]

Commit: 21e4f9a

[dimpase]

Author: Dima Pasechnik

[mwageringel]

comment:52

Here is another example that seems to be solved by this. This does not involve a lex ordering I think, but only happens over non-prime fields.

sage: R.<x,y> = GF(101^2)['X,Y'].quotient('X*Y - 1')
sage: R('X^3 - X^2*Y')
x^3 - x^2*y

However, I could not find any documentation whatsoever on this noTailReduction flag. I think we should avoid using low-level functionality like this which we do not understand, as this is hard to maintain. Upstream suggested to call kNF() instead, so we should at least try that, unless we know more about why this was not done in the first place. In fact, looking at the code in multi_polynomial_libsingular.pyx:

sage: R2.<x,y> = PolynomialRing(QQ, 2, order="lex")
sage: I2 = R2.ideal(["x^2 - x", "y^2 - y"])
sage: R2("x + y^2").reduce(I2)          # calls _groebner_strategy().normal_form()
x + y^2
sage: R2("x + y^2").reduce(I2.gens())   # calls kNF()
x + y

[dimpase]

comment:53

Frankly, the whole libsingular interface is hard to maintain, as it relies on undocumented functionality.

I went and read source to figure out what that flag does.

[dimpase]

comment:54

I think making use of kNF is a bit too major a rewrite of the libsingular interface for me to attempt it.

[mwageringel]

comment:55

Ok, fair enough. The whole groebner_strategy module directly uses Singular internals, so this patch is certainly in the spirit of that module. The module is only used for normal form computations, so it seems unlikely this change has unexpected effects that are worse than what we currently have.

Therefore, I am setting this to positive. The issue this ticket solves is a huge pain, as equality testing in some quotient rings is completely broken, which defeats the point of using quotient rings in the first place. Thanks for finding a solution.

[mwageringel]

Reviewer: Markus Wageringel

[vbraun]

Changed branch from u/dimpase/singular/singu_nf to 21e4f9a

[slel]

Changed commit from 21e4f9a to none

[slel]

comment:57

Replying to @dimpase:

And by the way there is this: #12529

Doctest added there (needs review).