Force tail reduction in polynomial quotient ring

[gh-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

[dimpase]

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

[gh-rachelplayer]

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

[dimpase]

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]

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

src/sage/libs/singular/ring.pyx

@@ -50 +50 @@
 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.

[klee]

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

[dimpase]

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]

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]

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]

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]

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]

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]

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...

[klee]

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]

Replying to klee:

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]

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]

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]

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]

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]

Replying to dimpase:

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

I looked at comment:16 and I don't see a conclusive argument there that checks that the internal libsingular ring does *not* have dp as its ordering. In ​singular_ring_new we see that the singular termorder string is looked up from the term_order.singular_str() attribute. That's what's used as key into the order_dict dictionary, with ​order_dict.get(s, ringorder_dp), so if a key is not found then you end up with dp as ordering. We have:

sage: R1.<x,y> = PolynomialRing(QQ, 2)
sage: R1.term_order().singular_str()
'dp'
sage: R2.<x,y> = PolynomialRing(QQ, 2, order="lex")
sage: R2.term_order().singular_str()
'lp'
sage: R3.<x,y> = PolynomialRing(QQ, 2, order="lp")
sage: R3.term_order().singular_str()
'lp'

so changing "lp" as a key in order_dict to "lex" would result in all three of these resulting into a libsingular ring with dp as ordering; where we are not observing failure of tail reduction on the examples we're trying.

It's a little strange that lp is also accepted as ordering, since that is not a value that Sage really knows how to work with otherwise. That explicitly coded, however: ​if name not in print_name_mapping.keys() and name not in singular_name_mapping.values() (that's why "baz" is not accepted as a term order and "lp" is)

So the way I read the code is that the currently proposed modification maps both degrevlex and lex (and also "lp") to ringorder_dp in libsingular.

We should probably be using examples where "lex" and "degrevlex" give different normal forms (and hopefully also different tailreduction results). I'm pretty sure the real problem is *not* in order_dict.

[dimpase]

Replying to nbruin:

[...]

We should probably be using examples where "lex" and "degrevlex" give different normal forms (and hopefully also different tailreduction results). I'm pretty sure the real problem is *not* in order_dict.

It seems that Volker might be right about the real problem: the tail is not rewritten fully simply because due to Singular's options in place when Singular's redtailBba function is called in normal_form(). It's undocumented code in Singular, and Sage's interface to it dates back to 2009. Perhaps malb can comment on it...

[dimpase]

And by the way there is this: #12529

[nbruin]

Replying to dimpase:

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_()
[...]

Beware, the ._singular_() code has to do with the Singular *interface*. Whatever it does is quite independent of libsingular. If you want to know what libsingular is doing, you have to query *it*. The design of libsingular doesn't make it particularly easy to dig up stuff about it interactively.

[dimpase]

Replying to nbruin:

Replying to dimpase:

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_()
[...]

Beware, the ._singular_() code has to do with the Singular *interface*. Whatever it does is quite independent of libsingular. If you want to know what libsingular is doing, you have to query *it*. The design of libsingular doesn't make it particularly easy to dig up stuff about it interactively.

Yes, I became aware of this. At some point on this ticket I started testing things after removing SAGE_LOCAL/bin/Singular :-)

[dimpase]

It's not obvious which of several prototypes of redtailBba is being called from normal_form(). So this is the question for the original authors: how does one change the default---which is, apparently, assuming that this is indeed not a bug, that the tail reduction is not carried out?

[malb]

Replying to dimpase:

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

src/sage/libs/singular/ring.pyx

@@ -50 +50 @@
 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.

Doesn't that simply revert back to "dp"? _order[idx] = order_dict.get(s, ringorder_dp)

[malb]

I'm guessing this is a question for the Singular team, they might not like tail reductions for lex?

[embray]

Moving all blocker/critical issues from 8.7 to 8.8.

[dimpase]

Replying to malb:

Replying to dimpase:

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

src/sage/libs/singular/ring.pyx

@@ -50 +50 @@
 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.

Doesn't that simply revert back to "dp"? _order[idx] = order_dict.get(s, ringorder_dp)

no, sorry, this actually just breaks the creation of the strategy thing completely, and resorts to reducing using an explicitly computed Groebner basis, as I sort of said later in comment:18.

[dimpase]

Replying to malb:

I'm guessing this is a question for the Singular team, they might not like tail reductions for lex?

Which version of redtailBba() do you call? In kutil.h one sees

KINLINE poly redtailBba (poly p,int end_pos,kStrategy strat,
       BOOLEAN normalize=FALSE);
poly redtailBba (LObject *L, int end_pos,kStrategy strat,
               BOOLEAN withT = FALSE,BOOLEAN normalize=FALSE);
poly redtailBba (TObject *T, int end_pos,kStrategy strat);

so that normalize=FALSE (and mysteriously named withT = FALSE) suggests you might be calling with these defaults.

The call is _p = redtailBba(_p, max_ind, self._strat) with cdef poly *_p, so it might be it picks up the 2nd prototype, with these defaults...

[malb]

Indeed, that's the one. However, this doesn't seem to make a difference:

src/sage/libs/singular/decl.pxd

@@ -890 +890 @@
     # head reduction
     poly *redNF(poly *p, int index, int nonorm, skStrategy *strat)
     # tail reduction
-    poly *redtailBba(poly *p, int index, skStrategy *strat)
+    poly *redtailBba(poly *p, int index, skStrategy *strat, bint normalize)
 
     cdef int CMD_1
     cdef int CMD_2

src/sage/libs/singular/groebner_strategy.pyx

@@ -288 +288 @@
         cdef int max_ind = 0
         cdef poly *_p = redNF(p_Copy(p._poly, self._parent._ring), max_ind, 0, self._strat)
         if likely(_p!=NULL):
-            _p = redtailBba(_p, max_ind, self._strat)
+            _p = redtailBba(_p, max_ind, self._strat, True)
         return new_MP(self._parent, _p)
 
 cdef class NCGroebnerStrategy(SageObject):
@@ -517 +517 @@
         if unlikely(self._parent._ring != currRing):
             rChangeCurrRing(self._parent._ring)
 
-        cdef int max_ind
+        cdef int max_ind = 0
         cdef poly *_p = redNF(p_Copy(p._poly, self._parent._ring), max_ind, 0, self._strat)
         if likely(_p!=NULL):
-            _p = redtailBba(_p, max_ind, self._strat)
+            _p = redtailBba(_p, max_ind, self._strat, True)
         return new_NCP(self._parent, _p)
 
 def unpickle_NCGroebnerStrategy0(I):

And I still get:

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

[dimpase]

Note that the class ​skStrategy has a noTailReduction member. Could it be that it's True by default?

[dimpase]

Replying to dimpase:

Note that the class ​skStrategy has a noTailReduction member. Could it be that it's True by default?

YES!

src/sage/libs/singular/decl.pxd

@@ -317 +317 @@
         void (*initEcartPair)(LObject * h, poly *f, poly *g, int ecartF, int ecartG)
         int (*posInLOld)(LObject * Ls, int Ll, LObject* Lo, skStrategy *strat)
         LObject P
+        bint noTailReduction
         ideal *Shdl
         ideal *D
         ideal *M

src/sage/libs/singular/groebner_strategy.pyx

@@ -127 +127 @@
         self._strat.sl = -1
         #- init local data struct
         initS(i, NULL, self._strat)
+        self._strat.noTailReduction = False
 

restores the sanity...

I'll make a patch.

[malb]

Nice find! It's strange that

sage: from sage.libs.singular.option import LibSingularOptions
sage: LibSingularOptions()["redTail"]
True

doesn't seem to affect it, but perhaps that only affects GB computations not manual normal forms.

[dimpase]

Replying to malb:

Nice find! It's strange that

sage: from sage.libs.singular.option import LibSingularOptions
sage: LibSingularOptions()["redTail"]
True

doesn't seem to affect it, but perhaps that only affects GB computations not manual normal forms.

What I don't understand is why self._strat.noTailReduction = False is needed explicitly, isn't the False value the default result of the initialisation of the class?

Or is it implementation-dependent?

[dimpase]

Replying to malb:

Nice find! It's strange that

sage: from sage.libs.singular.option import LibSingularOptions
sage: LibSingularOptions()["redTail"]
True

doesn't seem to affect it, but perhaps that only affects GB computations not manual normal forms.

Already creating a ring with lex term ordering option flips this option:

sage: from sage.libs.singular.option import LibSingularOptions
sage: LibSingularOptions()["redTail"]
True
sage: R2.<x,y> = PolynomialRing(QQ, 2, order="lex")
sage: LibSingularOptions()["redTail"]
False

in fact, I saw the following comment in kernel/GBEngine/kutil.cc:

/* alway use tailreduction, except:
  * - in local rings, - in lex order case, -in ring over extensions */

Thus, lex is special...

[dimpase]

I opened ​https://github.com/Singular/Sources/issues/920 to ask upstream for confirmation that we understand everything well here.

[nbruin]

Replying to dimpase:

I opened ​https://github.com/Singular/Sources/issues/920 to ask upstream for confirmation that we understand everything well here.

It doesn't look like github issues are a regular discussion forum for Singular issues, so you might want to send it elsewhere.

I can imagine that particularly for lex, tail reduction might be more expensive (e.g., there's generally an infinite number of monomials smaller than the leading one, so it's harder to bound the process). So there might be something to say for amending the documentation of reduce to allow for non-tailreduced outcomes and allowing for an option to set whether tailreduction is required or not.

[dimpase]

Replying to nbruin:

Replying to dimpase:

I opened ​https://github.com/Singular/Sources/issues/920 to ask upstream for confirmation that we understand everything well here.

It doesn't look like github issues are a regular discussion forum for Singular issues, so you might want to send it elsewhere.

in my experience, it's a place where one of the authors of sage/libs/singular/groebner_strategy looks very often, and replies quickly.

I can imagine that particularly for lex, tail reduction might be more expensive (e.g., there's generally an infinite number of monomials smaller than the leading one, so it's harder to bound the process). So there might be something to say for amending the documentation of reduce to allow for non-tailreduced outcomes and allowing for an option to set whether tailreduction is required or not.

yes, with lex Groebner bases are pretty bad, that's no surprise that one does not want full reduction (and there is no harm in postponing it, except one wants it, of course, too).

I'm just surprised that libSingular appears to change an option without asking.

[dimpase]

on ​https://github.com/Singular/Sources/issues/920 we get an advice to call kNF() rather than redNF() and redtailBba(), and indeed kNF() will call redNF() and redtailBba() for you, and do some other "admin" stuff with options.

The question is, why was this not done in the 1st place?

[nbruin]

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]

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

[malb]

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

[embray]

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]

by popular demand, a minimal patch.

---

New commits:

​21e4f9a

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

[gh-mwageringel]

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]

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]

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

[gh-mwageringel]

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.

[slelievre]

Replying to dimpase:

And by the way there is this: #12529

Doctest added there (needs review).