Trac #29625: support for weighted term orders in normal_basis

As a follow-up to #29543, this ticket changes the `normal_basis` method
of ideals to handle the case of weighted term orders.

With this change, the degree of the monomials in the normal basis is
taken with respect to the weighted degree (which agrees with Sage's
notion of degree).

{{{
sage: R.<x,y,z> = PolynomialRing(QQ, order=TermOrder('wdegrevlex', (1,
2, 3)))
sage: I = R.ideal(x*y^2 + x^5, z*y + x^3*y)
sage: I.normal_basis(degree=9)
[x^2*y^2*z, x^3*z^2, x*y*z^2, z^3]
sage: all(f.degree() == 9 for f in _)
True
}}}

This also came up in an [https://ask.sagemath.org/question/47623 Ask
SageMath question].

The implementation uses the Singular function [https://www.singular.uni-
kl.de/Manual/4-1-2/sing_377.htm#SEC417 weightKB].

URL: https://trac.sagemath.org/29625
Reported by: gh-mwageringel
Ticket author(s): Markus Wageringel
Reviewer(s): Travis Scrimshaw

Author: Release Manager

src/sage/rings/polynomial/multi_polynomial_ideal.py

@@ -2888,7 +2888,7 @@ def hilbert_numerator(self, grading = None, algorithm = 'sage'):
             raise ValueError("'algorithm' must be one of 'sage' or 'singular'")
 
     @require_field
-    def _normal_basis_libsingular(self, degree):
+    def _normal_basis_libsingular(self, degree, weights=None):
         r"""
         Return the normal basis for a given Groebner basis.
 
@@ -2899,6 +2899,9 @@ def _normal_basis_libsingular(self, degree):
 
         - ``degree`` -- ``None`` or integer
 
+        - ``weights`` -- tuple of positive integers (default: ``None``); if not
+          ``None``, compute the degree with respect to these weights
+
         OUTPUT:
 
         If ``degree`` is an integer, only the monomials of the given degree in
@@ -2913,6 +2916,9 @@ def _normal_basis_libsingular(self, degree):
             sage: J = R.ideal(x^2-2*x*z+5)
             sage: J.normal_basis(3)  # indirect doctest
             [z^3, y*z^2, x*z^2, y^2*z, x*y*z, y^3, x*y^2]
+            sage: [J._normal_basis_libsingular(d, (2, 2, 3)) for d in (0..8)]
+            [[1], [], [x, y], [z], [x*y, y^2], [x*z, y*z], [x*y^2, y^3, z^2],
+             [x*y*z, y^2*z], [x*y^3, y^4, x*z^2, y*z^2]]
 
         TESTS:
 
@@ -2926,7 +2932,13 @@ def _normal_basis_libsingular(self, degree):
         from sage.rings.polynomial.multi_polynomial_ideal_libsingular import kbase_libsingular
         from sage.rings.polynomial.multi_polynomial_sequence import PolynomialSequence
         gb = self._groebner_basis_libsingular()
-        res = kbase_libsingular(self.ring().ideal(gb), degree)
+        J = self.ring().ideal(gb)
+        if weights is None:
+            res = kbase_libsingular(J, degree)
+        else:
+            from sage.libs.singular.function_factory import ff
+            res = ff.weightKB(J, -1 if degree is None else degree,
+                              tuple(weights), attributes={J: {'isSB': 1}})
         if len(res) == 1 and res[0].is_zero():
             res = []
         return PolynomialSequence(self.ring(), res, immutable=True)
@@ -2943,7 +2955,8 @@ def normal_basis(self, degree=None, algorithm='libsingular',
         - ``degree`` -- integer (default: ``None``)
 
         - ``algorithm`` -- string (default: ``"libsingular"``); if not the
-          default, this will use the ``kbase()`` command from Singular
+          default, this will use the ``kbase()`` or ``weightKB()`` command from
+          Singular
 
         - ``singular`` -- the singular interpreter to use when ``algorithm`` is
           not ``"libsingular"`` (default: the default instance)
@@ -2974,18 +2987,35 @@ def normal_basis(self, degree=None, algorithm='libsingular',
             sage: [J.normal_basis(d, algorithm='singular') for d in (0..3)]
             [[1], [z, y, x], [z^2, y*z, x*z, x*y], [z^3, y*z^2, x*z^2, x*y*z]]
 
+        In case of a polynomial ring with a weighted term order, the degree of
+        the monomials is taken with respect to the weights.  ::
+
+            sage: T = TermOrder('wdegrevlex', (1, 2, 3))
+            sage: R.<x,y,z> = PolynomialRing(QQ, order=T)
+            sage: B = R.ideal(x*y^2 + x^5, z*y + x^3*y).normal_basis(9); B
+            [x^2*y^2*z, x^3*z^2, x*y*z^2, z^3]
+            sage: all(f.degree() == 9 for f in B)
+            True
+
         TESTS:
 
         Check that this method works over QQbar (:trac:`25351`)::
 
-            sage: P.<x,y,z> = QQbar[]
+            sage: R.<x,y,z> = QQbar[]
             sage: I = R.ideal(x^2+y^2+z^2-4, x^2+2*y^2-5, x*z-1)
             sage: I.normal_basis()
             [y*z^2, z^2, y*z, z, x*y, y, x, 1]
             sage: J = R.ideal(x^2+y^2+z^2-4, x^2+2*y^2-5)
             sage: [J.normal_basis(d) for d in (0..3)]
             [[1], [z, y, x], [z^2, y*z, x*z, x*y], [z^3, y*z^2, x*z^2, x*y*z]]
 
+        Check the option ``algorithm="singular"`` with a weighted term order::
+
+            sage: T = TermOrder('wdegrevlex', (1, 2, 3))
+            sage: S.<x,y,z> = PolynomialRing(GF(2), order=T)
+            sage: S.ideal(x^6 + y^3 + z^2).normal_basis(6, algorithm='singular')
+            [x^4*y, x^2*y^2, y^3, x^3*z, x*y*z, z^2]
+
         Check the deprecation::
 
             sage: R.<x,y> = PolynomialRing(QQ)
@@ -3001,15 +3031,22 @@ def normal_basis(self, degree=None, algorithm='libsingular',
             algorithm = degree
             degree = None
 
+        weights = tuple(x.degree() for x in self.ring().gens())
+        if all(w == 1 for w in weights):
+            weights = None
+
         if algorithm == 'libsingular':
-            return self._normal_basis_libsingular(degree)
+            return self._normal_basis_libsingular(degree, weights=weights)
         else:
             gb = self.groebner_basis()
             R = self.ring()
             if degree is None:
                 res = singular.kbase(R.ideal(gb))
-            else:
+            elif weights is None:
                 res = singular.kbase(R.ideal(gb), int(degree))
+            else:
+                res = singular.weightKB(R.ideal(gb), int(degree),
+                                        singular(weights, type='intvec'))
             return PolynomialSequence(R, [R(f) for f in res], immutable=True)