Improve performance of function fields

src/doc/en/reference/function_fields/index.rst

@@ -4,9 +4,8 @@ Algebraic Function Fields
 Sage allows basic computations with elements and ideals in orders of algebraic
 function fields over arbitrary constant fields. Advanced computations, like
 computing the genus or a basis of the Riemann-Roch space of a divisor, are
-available for function fields over finite fields, number fields, and `\QQbar`.
-
-A reference for the basic theory of algebraic function fields is [Stich2009]_.
+available for function fields over finite fields, number fields, and the
+algebraic closure of `\QQ`.
 
 .. toctree::
    :maxdepth: 1
@@ -22,4 +21,14 @@ A reference for the basic theory of algebraic function fields is [Stich2009]_.
    sage/rings/function_field/maps
    sage/rings/function_field/constructor
 
+A basic reference for the theory of algebraic function fields is [Stich2009]_.
+
+A Support Module
+----------------
+
+.. toctree::
+   :maxdepth: 1
+
+   sage/rings/function_field/hermite_form_polynomial
+
 .. include:: ../footer.txt

src/sage/matrix/matrix2.pyx

@@ -2386,14 +2386,14 @@ cdef class Matrix(Matrix1):
         r"""
         Computes the Pfaffian of ``self`` using the Baer-Faddeev-LeVerrier
         algorithm.
-        
+
         .. WARNING::
-        
+
             This method assumes that the base ring is an `\QQ`-algebra.
 
         OUTPUT:
 
-        - an element (possibly coerced) originated from the base ring of 
+        - an element (possibly coerced) originated from the base ring of
           ``self`` representing the Pfaffian
 
         EXAMPLES:
@@ -2409,7 +2409,7 @@ cdef class Matrix(Matrix1):
             ....:             (1/2, -3/2, -1, -5/2, -1/2, 0)])
             sage: A.pfaffian(algorithm='bfl')
             -1/2
-            
+
         TESTS::
 
             sage: A = random_matrix(ZZ[x], 6)
@@ -15313,7 +15313,7 @@ cdef class Matrix(Matrix1):
         cdef Py_ssize_t i = 0
         cdef Py_ssize_t j = 0
 
-        cdef Py_ssize_t k, l
+        cdef Py_ssize_t k, l, c
 
         if transformation:
             from sage.matrix.constructor import identity_matrix
@@ -15349,9 +15349,9 @@ cdef class Matrix(Matrix1):
                             U.set_unsafe(k, c, p * Ukc + q * Ulc)
                             U.set_unsafe(l, c, (-f) * Ukc + e * Ulc)
                 if i != k:
-                    A.swap_rows(i,k)
+                    A.swap_rows_c(i,k)
                     if transformation:
-                        U.swap_rows(i,k)
+                        U.swap_rows_c(i,k)
                 pivot_cols.append(j)
                 i += 1
             j += 1
@@ -15366,9 +15366,9 @@ cdef class Matrix(Matrix1):
                 coeff = normalization(pivot)
                 for c in range(j,n):
                     A.set_unsafe(i, c, A.get_unsafe(i,c) * coeff)
-                    if transformation:
-                        for c in range(m):
-                            U.set_unsafe(i, c, U.get_unsafe(i,c) * coeff)
+                if transformation:
+                    for c in range(m):
+                        U.set_unsafe(i, c, U.get_unsafe(i,c) * coeff)
 
             pivot = A.get_unsafe(i,j)
             for k in range(i):

src/sage/rings/function_field/function_field.py

@@ -3534,6 +3534,7 @@ def _maximal_order_basis(self):
         in some algorithms.
         """
         from sage.matrix.constructor import matrix
+        from .hermite_form_polynomial import reversed_hermite_form
 
         k = self.constant_base_field()
         K = self.base_field() # rational function field
@@ -3608,14 +3609,8 @@ def _maximal_order_basis(self):
         basis_V = [to_V(bvec) for bvec in _basis]
         l = lcm([vvec.denominator() for vvec in basis_V])
 
-        # Why do we have 'reversed' here? I don't know. But without it, the
-        # time to get hermite_form_reversed dramatically increases.
-        _mat = matrix([[coeff.numerator() for coeff in l*v] for v in reversed(basis_V)])
-
-        # compute the reversed hermite form
-        _mat.reverse_rows_and_columns()
-        _mat._hermite_form_euclidean(normalization=lambda p: ~p.lc())
-        _mat.reverse_rows_and_columns()
+        _mat = matrix([[coeff.numerator() for coeff in l*v] for v in basis_V])
+        reversed_hermite_form(_mat)
 
         basis = [fr_V(v) / l for v in _mat if not v.is_zero()]
         return basis

src/sage/rings/function_field/hermite_form_polynomial.pyx

@@ -0,0 +1,187 @@
+r"""
+Hermite form computation for function fields
+
+This module provides an optimized implementation of the algorithm computing
+Hermite forms of matrices over polynomials. This is the workhorse of the
+function field machinery of Sage.
+
+EXAMPLES::
+
+    sage: P.<x> = PolynomialRing(QQ)
+    sage: A = matrix(P,3,[-(x-1)^((i-j+1) % 3) for i in range(3) for j in range(3)])
+    sage: A
+    [        -x + 1             -1 -x^2 + 2*x - 1]
+    [-x^2 + 2*x - 1         -x + 1             -1]
+    [            -1 -x^2 + 2*x - 1         -x + 1]
+    sage: from sage.rings.function_field.hermite_form_polynomial import reversed_hermite_form
+    sage: B = copy(A)
+    sage: U = reversed_hermite_form(B, transformation=True)
+    sage: U * A == B
+    True
+    sage: B
+    [x^3 - 3*x^2 + 3*x - 2                     0                     0]
+    [                    0 x^3 - 3*x^2 + 3*x - 2                     0]
+    [        x^2 - 2*x + 1                 x - 1                     1]
+
+The function :func:`reversed_hermite_form` computes the reversed hermite form,
+which is reversed both row-wise and column-wise from the usual hermite form.
+Let us check it::
+
+    sage: A.reverse_rows_and_columns()
+    sage: C = copy(A.hermite_form())
+    sage: C.reverse_rows_and_columns()
+    sage: C
+    [x^3 - 3*x^2 + 3*x - 2                    0                     0]
+    [                    0 x^3 - 3*x^2 + 3*x - 2                     0]
+    [        x^2 - 2*x + 1                 x - 1                     1]
+    sage: C == B
+    True
+
+AUTHORS:
+
+- Kwankyu Lee (2021-05-21): initial version
+
+"""
+#*****************************************************************************
+#       Copyright (C) 2021 Kwankyu Lee <ekwankyu@gmail.com>
+#
+#  Distributed under the terms of the GNU General Public License (GPL)
+#  as published by the Free Software Foundation; either version 2 of
+#  the License, or (at your option) any later version.
+#                  http://www.gnu.org/licenses/
+#*****************************************************************************
+
+from sage.matrix.matrix cimport Matrix
+from sage.rings.polynomial.polynomial_element cimport Polynomial
+from sage.matrix.constructor import identity_matrix
+
+def reversed_hermite_form(Matrix mat, bint transformation=False):
+    """
+    Transform the matrix in place to reversed hermite normal form and
+    optionally return the transformation matrix.
+
+    INPUT:
+
+    - ``transformation`` -- boolean (default: ``False``); if ``True``,
+      return the transformation matrix
+
+    EXAMPLES::
+
+        sage: from sage.rings.function_field.hermite_form_polynomial import reversed_hermite_form
+        sage: P.<x> = PolynomialRing(QQ)
+        sage: A = matrix(P,3,[-(x-1)^((i-2*j) % 4) for i in range(3) for j in range(3)])
+        sage: A
+        [                    -1         -x^2 + 2*x - 1                     -1]
+        [                -x + 1 -x^3 + 3*x^2 - 3*x + 1                 -x + 1]
+        [        -x^2 + 2*x - 1                     -1         -x^2 + 2*x - 1]
+        sage: B = copy(A)
+        sage: U = reversed_hermite_form(B, transformation=True)
+        sage: U * A == B
+        True
+        sage: B
+        [                        0                         0                         0]
+        [                        0 x^4 - 4*x^3 + 6*x^2 - 4*x                         0]
+        [                        1             x^2 - 2*x + 1                         1]
+    """
+    cdef Matrix A = mat
+    cdef Matrix U
+
+    cdef Py_ssize_t m = A.nrows()
+    cdef Py_ssize_t n = A.ncols()
+    cdef Py_ssize_t i = m - 1
+    cdef Py_ssize_t j = n - 1
+    cdef Py_ssize_t k, l, c, ip
+
+    cdef Polynomial a, b, d, p, q, e, f, Aic, Alc, Uic, Ulc
+    cdef Polynomial zero = A.base_ring().zero()
+
+    cdef int di, dk
+
+    cdef list pivot_cols = []
+
+    if transformation:
+        U = <Matrix> identity_matrix(A.base_ring(), m)
+
+    while j >= 0:
+        # find the first row with nonzero entry in jth column
+        k = i
+        while k >= 0 and not A.get_unsafe(k, j):
+            k -= 1
+        if k >= 0:
+            # swap the kth row with the ith row
+            if k < i:
+                A.swap_rows_c(i, k)
+                if transformation:
+                    U.swap_rows_c(i, k)
+            k -= 1
+            # put the row with the smalllest degree to the ith row
+            di = A.get_unsafe(i, j).degree()
+            while k >= 0 and A.get_unsafe(k, j):
+                dk = A.get_unsafe(k, j).degree()
+                if dk < di:
+                    A.swap_rows_c(i, k)
+                    di = dk
+                    if transformation:
+                        U.swap_rows_c(i, k)
+                k -= 1
+            l = i - 1
+            while True:
+                # find a row with nonzero entry in the jth column
+                while l >= 0 and not A.get_unsafe(l, j):
+                    l -= 1
+                if l < 0:
+                    break
+
+                a = <Polynomial> A.get_unsafe(i, j)
+                b = <Polynomial> A.get_unsafe(l, j)
+                d, p, q = a.xgcd(b) # p * a + q * b = d = gcd(a,b)
+                e = a // d
+                f = -(b // d)
+
+                A.set_unsafe(i, j, d)
+                A.set_unsafe(l, j, zero)
+
+                for c in range(j):
+                    Aic = <Polynomial> A.get_unsafe(i, c)
+                    Alc = <Polynomial> A.get_unsafe(l, c)
+                    A.set_unsafe(i, c, p * Aic + q * Alc)
+                    A.set_unsafe(l, c, f * Aic + e * Alc)
+                if transformation:
+                    for c in range(m):
+                        Uic = <Polynomial> U.get_unsafe(i, c)
+                        Ulc = <Polynomial> U.get_unsafe(l, c)
+                        U.set_unsafe(i, c, p * Uic + q * Ulc)
+                        U.set_unsafe(l, c, f * Uic + e * Ulc)
+            pivot_cols.append(j)
+            i -= 1
+        j -= 1
+
+    # reduce entries below pivots
+    for i in range(len(pivot_cols)):
+        j = pivot_cols[i]
+        ip = m - 1 - i
+        pivot = <Polynomial> A.get_unsafe(ip, j)
+
+        # normalize the leading coefficient to one
+        coeff = ~pivot.lc()
+        for c in range(j + 1):
+            A.set_unsafe(ip, c, <Polynomial> A.get_unsafe(ip, c) * coeff)
+        if transformation:
+            for c in range(m):
+                U.set_unsafe(ip, c, <Polynomial> U.get_unsafe(ip, c) * coeff)
+
+        pivot = <Polynomial> A.get_unsafe(ip,j)
+        for k in range(ip + 1, m):
+            q = -(<Polynomial> A.get_unsafe(k, j) // pivot)
+            if q:
+                for c in range(j + 1):
+                    A.set_unsafe(k, c, <Polynomial> A.get_unsafe(k, c)
+                                 + q * <Polynomial> A.get_unsafe(ip, c))
+                if transformation:
+                    for c in range(m):
+                        U.set_unsafe(k, c, <Polynomial> U.get_unsafe(k, c)
+                                     + q * <Polynomial> U.get_unsafe(ip, c))
+
+    if transformation:
+        return U
+

src/sage/rings/function_field/ideal.py

@@ -110,6 +110,8 @@
 
 from .divisor import divisor
 
+from .hermite_form_polynomial import reversed_hermite_form
+
 
 class FunctionFieldIdeal(Element):
     """
@@ -1591,10 +1593,7 @@ def intersect(self, other):
         M = block_matrix([[I,I],[A,O],[O,B]])
 
         # reversed Hermite form
-        M.reverse_rows_and_columns()
-        U = M._hermite_form_euclidean(transformation=True,
-                                      normalization=lambda p: ~p.lc())
-        U.reverse_rows_and_columns()
+        U = reversed_hermite_form(M, transformation=True)
 
         vecs = [U[i][:n] for i in range(n)]
 

src/sage/rings/function_field/order.py

@@ -137,6 +137,8 @@
     FunctionFieldIdealInfinite_rational,
     FunctionFieldIdealInfinite_polymod)
 
+from .hermite_form_polynomial import reversed_hermite_form
+
 
 class FunctionFieldOrder_base(CachedRepresentation, Parent):
     """
@@ -1475,10 +1477,8 @@ def _ideal_from_vectors_and_denominator(self, vecs, d=1, check=True):
         # so that we get a unique hnf. Here the hermite form
         # algorithm also makes the pivots monic.
 
-        # compute the reverse hermite form with zero rows deleted
-        mat.reverse_rows_and_columns()
-        mat._hermite_form_euclidean(normalization=lambda p: ~p.lc())
-        mat.reverse_rows_and_columns()
+        # compute the reversed hermite form with zero rows deleted
+        reversed_hermite_form(mat)
         i = 0
         while i < mat.nrows() and mat.row(i).is_zero():
             i += 1
@@ -2605,9 +2605,7 @@ def ideal_with_gens_over_base(self, gens):
                 k = x * k
 
                 h2 = block_matrix([[h],[k]])
-                h2.reverse_rows_and_columns()
-                h2._hermite_form_euclidean(normalization=lambda p: ~p.lc())
-                h2.reverse_rows_and_columns()
+                reversed_hermite_form(h2)
                 i = 0
                 while i < h2.nrows() and h2.row(i).is_zero():
                     i += 1