From 4ecde037207b0c97fbad6001a75d2cb3f9fff2a7 Mon Sep 17 00:00:00 2001 From: Vincent Macri Date: Wed, 10 Jun 2026 09:04:21 -0600 Subject: [PATCH 1/4] Begin work on new function field integral basis code --- .../function_field/function_field_polymod.py | 188 +++++++++--------- 1 file changed, 97 insertions(+), 91 deletions(-) diff --git a/src/sage/rings/function_field/function_field_polymod.py b/src/sage/rings/function_field/function_field_polymod.py index 82f39fb0b7a..54ad3ba1886 100644 --- a/src/sage/rings/function_field/function_field_polymod.py +++ b/src/sage/rings/function_field/function_field_polymod.py @@ -37,6 +37,7 @@ from sage.categories.number_fields import NumberFields from sage.misc.cachefunc import cached_method from sage.misc.lazy_import import LazyImport +from sage.misc.misc_c import prod from sage.rings.function_field.element import FunctionFieldElement from sage.rings.function_field.element_polymod import FunctionFieldElement_polymod from sage.rings.integer import Integer @@ -2404,97 +2405,102 @@ def _maximal_order_basis(self) -> tuple[FunctionFieldElement]: sage: F.maximal_order_infinite().basis() (1, 1/x*y, 1/x^2*y^2, 1/x^3*y^3, 1/x^4*y^4) """ - from sage.libs.singular.function import lib, singular_function - - k = self.constant_base_field() - K = self.base_field() # rational function field - - # Construct the defining polynomial of the function field as a - # two-variate polynomial g in the ring k[y,x] where k is the constant - # base field. - S, (y, x) = PolynomialRing(k, names='y, x', order='degrevlex').objgens() - v = self.polynomial().list() - g = sum([v[i].numerator().subs(x) * y**i for i in range(len(v))]) - - if self.is_global(): - # If the constant base field is a prime field then we can use - # Singular's integralBasis function, which is much faster. - if k.is_prime_field(): - lib('integralbasis.lib') - integral_basis = singular_function('integralBasis') - - singular_basis, s = integral_basis(g, 1, 'normal') - hom = S.hom([self.gen(), self(K.gen())]) - return tuple(self(hom(b) / hom(s)) for b in singular_basis) - - from sage.env import SAGE_EXTCODE - lib(SAGE_EXTCODE + '/singular/function_field/core.lib') - normalize = singular_function('core_normalize') - - # Singular "normalP" algorithm assumes affine domain over - # a prime field. So we construct gflat lifting g as in - # k_prime[yy,xx,zz]/(k_poly) where k = k_prime[zz]/(k_poly) - R = PolynomialRing(k.prime_subfield(), names='yy,xx,zz') - gflat = R.zero() - for m in g.monomials(): - c = g.monomial_coefficient(m).polynomial('zz') - gflat += R(c) * R(m) # R(m) is a monomial in yy and xx - - k_poly = R(k.polynomial('zz')) - - # invoke Singular - pols_in_R = normalize(R.ideal([k_poly, gflat])) - - # reconstruct polynomials in S - h = R.hom([y, x, k.gen()], S) - pols_in_S = [h(f) for f in pols_in_R] - else: - # Call Singular. Singular's "normal" function returns a basis - # of the integral closure of k(x,y)/(g) as a k[x,y]-module. - pols_in_S = _singular_normal(S.ideal(g))[0] - - from sage.matrix.constructor import matrix - - from .hermite_form_polynomial import reversed_hermite_form - - # reconstruct the polynomials in the function field - x = K.gen() - y = self.gen() - pols = [] - for f in pols_in_S: - p = f.polynomial(S.gen(0)) - s = 0 - for i in range(p.degree() + 1): - s += p[i].subs(x) * y**i - pols.append(s) - - # Now if pols = [g1,g2,...gn,g0], then the g1/g0,g2/g0,...,gn/g0, - # and g0/g0=1 are the module generators of the integral closure - # of the equation order Sb = k[xb,yb] in its fraction field, - # that is, the function field. The integral closure of k[x] - # is then obtained by multiplying these generators with powers of y - # as the equation order itself is an integral extension of k[x]. - d = ~ pols[-1] - _basis = [] - for f in pols: - b = d * f - for i in range(self.degree()): - _basis.append(b) - b *= y - - # Finally we reduce _basis to get a basis over k[x]. This is done of - # course by Hermite normal form computation. Here we apply a trick to - # get a basis that starts with 1 and is ordered in increasing - # y-degrees. The trick is to use the reversed Hermite normal form. - # Note that it is important that the overall denominator l lies in k[x]. - V, fr_V, to_V = self.free_module() - basis_V = [to_V(bvec) for bvec in _basis] - l = lcm([vvec.denominator() for vvec in basis_V]) - - _mat = matrix([[coeff.numerator() for coeff in l * v] for v in basis_V]) - reversed_hermite_form(_mat) - - return tuple(fr_V(v) / l for v in _mat if not v.is_zero()) + Fq = self.constant_base_field() + Fqx = Fq.polynomial_ring() + T = self._polynomial.change_ring(Fqx) + n = T.degree() + D = T.discriminant() + theta = self.gen() + + # Step 1: Factor discriminant of polynomial + one = D.parent().one() + D0 = one + F = one + F_factors = [] + for p, m in D.factor(): + if m % 2 == 0: + pm = p ** (m // 2) + F *= pm + F_factors.append((p, pm)) + else: + D0 *= p + if m > 1: + pm = p ** ((m - 1) // 2) + F *= pm + F_factors.append((p, pm)) + #assert D == D0 * F**2 + assert (D0 * F**2).divides(D) and (D // (D0 * F**2)).degree() == 0 + + # Step 2: Initialize + omega = [theta**i for i in range(n)] + + print(F_factors) + + ## Step 3: Loop on factors of F + #if F.is_one(): + while F_factors: + p, pm = F_factors.pop() + print('Step 3:', F_factors, p, pm) + + # Step 4: Factor modulo p + Fqx_quotient_p = Fqx.quotient(p) + T_bar = T.change_ring(Fqx_quotient_p) + + t_bar = [] + e = [] + for ti_bar, ei in T_bar.factor(): + t_bar.append(ti_bar) + e.append(ei) + assert ti_bar.base_ring() == Fqx_quotient_p + + lift = Fqx_quotient_p.lifting_map() + g_bar = prod(t_bar) + g = g_bar.change_ring(lift) + h_bar = T_bar // g_bar + h = h_bar.change_ring(lift) + f = (g * h - T) // p + f_bar = f.change_ring(Fqx_quotient_p) + Z_bar = f_bar.gcd(g_bar).gcd(h_bar) + U_bar = T_bar // Z_bar + U = U_bar.change_ring(lift) + Z = Z_bar.change_ring(lift) + m = Z.degree() + print(f'Step 4: {Z=}') + + # Step 5: Apply Dedekind + if m.is_zero(): + # 𝒪 is p-maximal + F = F // pm + continue + + v = [(omega[i] * U(theta)).list() for i in range(m)] + v.extend((p * omega[j]).list() for j in range(n)) + from sage.matrix.constructor import matrix + M = matrix(v) + assert M.nrows() == n + m + assert M.ncols() == n + print(M) + H = M.hermite_form(include_zero_rows=False) + print() + print(H) + + for i in range(n): + omega[i] = H.row(i) / p + + # Step 6: Is the new order p-maximal? + if not (p**(m + 1)).divides(F): + # 𝒪 is p-maximal + F = F // pm + continue + + # Step 7: Compute radical + q = p.degree() + return omega, Fqx_quotient_p + while q < n: + q *= p.degree() + print(q, p.degree()) + + return tuple(omega) @cached_method def equation_order(self): From 2af36f7e531714dc1825cc6e7a980b59d730fb44 Mon Sep 17 00:00:00 2001 From: Vincent Macri Date: Thu, 11 Jun 2026 10:22:49 -0600 Subject: [PATCH 2/4] Up to step 6 in Cohen --- .../rings/function_field/element_polymod.pyi | 7 ++- .../rings/function_field/element_polymod.pyx | 3 ++ .../function_field/function_field_polymod.py | 45 ++++++++++++++----- 3 files changed, 41 insertions(+), 14 deletions(-) diff --git a/src/sage/rings/function_field/element_polymod.pyi b/src/sage/rings/function_field/element_polymod.pyi index fad1851c92a..ecee2aa203f 100644 --- a/src/sage/rings/function_field/element_polymod.pyi +++ b/src/sage/rings/function_field/element_polymod.pyi @@ -1,6 +1,6 @@ +from collections.abc import Iterable from typing import Any -from sage.structure.richcmp import richcmp -from sage.structure.element import FieldElement + from sage.rings.function_field.element import FunctionFieldElement class FunctionFieldElement_polymod(FunctionFieldElement): @@ -37,6 +37,9 @@ class FunctionFieldElement_polymod(FunctionFieldElement): def __invert__(self) -> Any: ... + def __iter__(self) -> Iterable: + ... + def list(self) -> list: ... diff --git a/src/sage/rings/function_field/element_polymod.pyx b/src/sage/rings/function_field/element_polymod.pyx index 53d3e13d287..2cb404e7013 100644 --- a/src/sage/rings/function_field/element_polymod.pyx +++ b/src/sage/rings/function_field/element_polymod.pyx @@ -231,6 +231,9 @@ cdef class FunctionFieldElement_polymod(FunctionFieldElement): P = self._parent return P(self._x.xgcd(P._polynomial)[1]) + def __iter__(self): + yield from self.list() + cpdef list list(self): """ Return the list of the coefficients representing the element. diff --git a/src/sage/rings/function_field/function_field_polymod.py b/src/sage/rings/function_field/function_field_polymod.py index 54ad3ba1886..502efa109a6 100644 --- a/src/sage/rings/function_field/function_field_polymod.py +++ b/src/sage/rings/function_field/function_field_polymod.py @@ -2405,6 +2405,8 @@ def _maximal_order_basis(self) -> tuple[FunctionFieldElement]: sage: F.maximal_order_infinite().basis() (1, 1/x*y, 1/x^2*y^2, 1/x^3*y^3, 1/x^4*y^4) """ + from sage.modules.free_module_element import vector + from sage.matrix.constructor import matrix Fq = self.constant_base_field() Fqx = Fq.polynomial_ring() T = self._polynomial.change_ring(Fqx) @@ -2433,6 +2435,7 @@ def _maximal_order_basis(self) -> tuple[FunctionFieldElement]: # Step 2: Initialize omega = [theta**i for i in range(n)] + power_basis = vector(omega) print(F_factors) @@ -2455,7 +2458,10 @@ def _maximal_order_basis(self) -> tuple[FunctionFieldElement]: lift = Fqx_quotient_p.lifting_map() g_bar = prod(t_bar) - g = g_bar.change_ring(lift) + #g = g_bar.change_ring(lift) + #print('g1', g) + #print('g2', prod(ti_bar.change_ring(lift) for ti_bar in t_bar)) + g = prod(ti_bar.change_ring(lift) for ti_bar in t_bar) h_bar = T_bar // g_bar h = h_bar.change_ring(lift) f = (g * h - T) // p @@ -2474,18 +2480,18 @@ def _maximal_order_basis(self) -> tuple[FunctionFieldElement]: continue v = [(omega[i] * U(theta)).list() for i in range(m)] + print(f'{U=}') + print(f'{v=}') v.extend((p * omega[j]).list() for j in range(n)) - from sage.matrix.constructor import matrix - M = matrix(v) + M = matrix(Fqx, v) assert M.nrows() == n + m assert M.ncols() == n - print(M) H = M.hermite_form(include_zero_rows=False) - print() - print(H) + + assert H.det().divides(D) for i in range(n): - omega[i] = H.row(i) / p + omega[i] = (H.row(i) / p).dot_product(power_basis) # Step 6: Is the new order p-maximal? if not (p**(m + 1)).divides(F): @@ -2494,11 +2500,26 @@ def _maximal_order_basis(self) -> tuple[FunctionFieldElement]: continue # Step 7: Compute radical - q = p.degree() - return omega, Fqx_quotient_p - while q < n: - q *= p.degree() - print(q, p.degree()) + SP = matrix(Fqx, n, n, lambda i, j: (omega[i] * omega[j]).trace()) + return SP, Fqx_quotient_p + + #q = p + #while q.degree() < n: + # q *= p + q = self.characteristic() + + A = matrix(Fqx_quotient_p, n) + suborder = self.order_with_basis(omega) + print() + print() + for j in range(n): + coord_vec = suborder.coordinate_vector(omega[j]**q) + print(f'{coord_vec=}') + A.set_row(j, coord_vec) + print(j, A) + print() + return omega, q, Fqx, Fqx_quotient_p, A, p + return tuple(omega) From fab647fbf89e3efe58eae8a3e6bbffbb2963f4fd Mon Sep 17 00:00:00 2001 From: Vincent Macri Date: Thu, 11 Jun 2026 15:16:39 -0600 Subject: [PATCH 3/4] Naming --- .../function_field/function_field_polymod.py | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/src/sage/rings/function_field/function_field_polymod.py b/src/sage/rings/function_field/function_field_polymod.py index 502efa109a6..4f2eec1afdc 100644 --- a/src/sage/rings/function_field/function_field_polymod.py +++ b/src/sage/rings/function_field/function_field_polymod.py @@ -2446,17 +2446,17 @@ def _maximal_order_basis(self) -> tuple[FunctionFieldElement]: print('Step 3:', F_factors, p, pm) # Step 4: Factor modulo p - Fqx_quotient_p = Fqx.quotient(p) - T_bar = T.change_ring(Fqx_quotient_p) + Fqx_mod_p = Fqx.quotient(p) + T_bar = T.change_ring(Fqx_mod_p) t_bar = [] e = [] for ti_bar, ei in T_bar.factor(): t_bar.append(ti_bar) e.append(ei) - assert ti_bar.base_ring() == Fqx_quotient_p + assert ti_bar.base_ring() == Fqx_mod_p - lift = Fqx_quotient_p.lifting_map() + lift = Fqx_mod_p.lifting_map() g_bar = prod(t_bar) #g = g_bar.change_ring(lift) #print('g1', g) @@ -2465,7 +2465,7 @@ def _maximal_order_basis(self) -> tuple[FunctionFieldElement]: h_bar = T_bar // g_bar h = h_bar.change_ring(lift) f = (g * h - T) // p - f_bar = f.change_ring(Fqx_quotient_p) + f_bar = f.change_ring(Fqx_mod_p) Z_bar = f_bar.gcd(g_bar).gcd(h_bar) U_bar = T_bar // Z_bar U = U_bar.change_ring(lift) @@ -2501,14 +2501,15 @@ def _maximal_order_basis(self) -> tuple[FunctionFieldElement]: # Step 7: Compute radical SP = matrix(Fqx, n, n, lambda i, j: (omega[i] * omega[j]).trace()) - return SP, Fqx_quotient_p + lkm = SP.left_kernel_matrix() + return SP, p, Fqx_mod_p, omega, g, theta #q = p #while q.degree() < n: # q *= p q = self.characteristic() - A = matrix(Fqx_quotient_p, n) + A = matrix(Fqx_mod_p, n) suborder = self.order_with_basis(omega) print() print() @@ -2518,7 +2519,7 @@ def _maximal_order_basis(self) -> tuple[FunctionFieldElement]: A.set_row(j, coord_vec) print(j, A) print() - return omega, q, Fqx, Fqx_quotient_p, A, p + return omega, q, Fqx, Fqx_mod_p, A, p return tuple(omega) From 3393e1833919fda0eef3f79f9d8567204be6b5a9 Mon Sep 17 00:00:00 2001 From: Vincent Macri Date: Thu, 11 Jun 2026 17:50:13 -0600 Subject: [PATCH 4/4] Work --- .../function_field/function_field_polymod.py | 87 ++++++++++++++----- 1 file changed, 66 insertions(+), 21 deletions(-) diff --git a/src/sage/rings/function_field/function_field_polymod.py b/src/sage/rings/function_field/function_field_polymod.py index 4f2eec1afdc..ef63b9c991c 100644 --- a/src/sage/rings/function_field/function_field_polymod.py +++ b/src/sage/rings/function_field/function_field_polymod.py @@ -2405,8 +2405,45 @@ def _maximal_order_basis(self) -> tuple[FunctionFieldElement]: sage: F.maximal_order_infinite().basis() (1, 1/x*y, 1/x^2*y^2, 1/x^3*y^3, 1/x^4*y^4) """ - from sage.modules.free_module_element import vector + + def expand_basis(M): + # Cohen Algorithm 2.3.6 (Supplement a basis) + n = M.nrows() + k = M.ncols() + assert k <= n + + # Step 1: Initialize + s = 0 + B = matrix(M.base_ring().fraction_field(), n, n, 1) + + while True: + # Step 2: Finished? + if s == k - 1: + return B + + # Step 3: Search for non-zero + s += 1 + t = None + for j in range(s, n): + if M[j, s] != 0: + t = j + d = ~M[t, s] + if t is None: + raise ValueError('M is of rank less than k') + + # Step 4: Modify basis and eliminate + B.set_column(t, B.column(s)) + B.set_column(s, M.column(s)) + for j in range(s + 1, k): + M[s, j] = d * M[t, j] + M[t, j] = M[s, j] + for i in range(n): + if i == s or i == t: + continue + M[i, j] = M[i, j] - M[i, s] * M[t, j] + from sage.matrix.constructor import matrix + from sage.modules.free_module_element import vector Fq = self.constant_base_field() Fqx = Fq.polynomial_ring() T = self._polynomial.change_ring(Fqx) @@ -2500,26 +2537,34 @@ def _maximal_order_basis(self) -> tuple[FunctionFieldElement]: continue # Step 7: Compute radical - SP = matrix(Fqx, n, n, lambda i, j: (omega[i] * omega[j]).trace()) - lkm = SP.left_kernel_matrix() - return SP, p, Fqx_mod_p, omega, g, theta - - #q = p - #while q.degree() < n: - # q *= p - q = self.characteristic() - - A = matrix(Fqx_mod_p, n) - suborder = self.order_with_basis(omega) - print() - print() - for j in range(n): - coord_vec = suborder.coordinate_vector(omega[j]**q) - print(f'{coord_vec=}') - A.set_row(j, coord_vec) - print(j, A) - print() - return omega, q, Fqx, Fqx_mod_p, A, p + Mp = matrix(Fqx_mod_p, n, n, lambda i, j: (omega[i] * omega[j]).trace()) + ker = Mp.right_kernel_matrix().transpose() + ker = matrix([lift(a) for a in r] for r in ker) + + # Step 8: Compute new basis mod p + l = ker.ncols() + beta_matrix = expand_basis(ker) + beta_matrix_lift = matrix([lift(a) for a in r] for r in beta_matrix) + #beta = [c.dot_product(vector(omega)) for c in beta_matrix_lift.columns()] # TODO: Should this be power basis? + beta = [vector(c.dot_product(vector(omega)).list()) for c in beta_matrix_lift.columns()] # TODO: Should this be power basis? + + # Step 9: Compute big matrix + alpha = [beta[i] for i in range(l)] + alpha.extend(p * beta[i] for i in range(l, n)) + + print(f'{alpha=}') + print(f'{beta=}') + + B = matrix([o * a for a in alpha for o in omega]).transpose() + print(f'{B=}') + A = matrix(alpha).transpose() + print(f'{A=}') + + return A, B, omega, alpha, Fqx_mod_p + #C = A.solve_left(B) + + return omega, alpha, Fqx_mod_p, C + return tuple(omega)