From bdb62fdd9fdfa8a81265ae225dae6ae1f841e3e9 Mon Sep 17 00:00:00 2001 From: Montessinos Mickael Gerard Bernard Date: Wed, 13 Mar 2024 13:03:43 +0200 Subject: [PATCH] Added algorithms related to splitting vector bundles --- README.rst | 2 +- makefile | 3 + setup.py | 9 +- vector_bundle/algebras.py | 559 ++++++++++++++++++++ vector_bundle/constructions.py | 33 +- vector_bundle/ext_group.py | 25 +- vector_bundle/function_field_utility.py | 675 ++++++++++++++++++++++-- vector_bundle/hom_bundle.py | 268 +++++++++- vector_bundle/vector_bundle.py | 562 ++++++++++++++++---- 9 files changed, 1968 insertions(+), 168 deletions(-) create mode 100644 vector_bundle/algebras.py diff --git a/README.rst b/README.rst index 7c43dd3..865c698 100644 --- a/README.rst +++ b/README.rst @@ -14,7 +14,7 @@ Local install from source Download the source from the git repository:: - $ git clone https://git.mif.vu.lt/mimo7829/vector-bundles.git + $ git clone https://git.disroot.org/montessiel/vector-bundles-sagemath.git Change to the root directory and run:: diff --git a/makefile b/makefile index 146fd76..8bd202b 100644 --- a/makefile +++ b/makefile @@ -20,6 +20,9 @@ develop: test: $(SAGE) setup.py test +debug: + $(SAGE) setup.py debug + coverage: $(SAGE) -coverage $(PACKAGE)/* diff --git a/setup.py b/setup.py index f20e898..1b2ef17 100644 --- a/setup.py +++ b/setup.py @@ -18,6 +18,13 @@ class SageTest(TestCommand): if errno != 0: sys.exit(1) +# For the tests +class SageDebug(TestCommand): + def run_tests(self): + errno = os.system("sage -t --debug --force-lib vector_bundle") + if errno != 0: + sys.exit(1) + setup( name = "vector_bundle", version = readfile("VERSION").strip(), # the VERSION file is shared with the documentation @@ -44,7 +51,7 @@ setup( ], # classifiers list: https://pypi.python.org/pypi?%3Aaction=list_classifiers keywords = "Algebraic Geometry Number Theory Curves Vector Bundles", packages = ['vector_bundle'], - cmdclass = {'test': SageTest}, # adding a special setup command for tests + cmdclass = {'test': SageTest, 'debug': SageDebug}, # adding a special setup command for tests and debugs setup_requires = ['sage-package'], install_requires = ['sage-package', 'sphinx'], ) diff --git a/vector_bundle/algebras.py b/vector_bundle/algebras.py new file mode 100644 index 0000000..8504a51 --- /dev/null +++ b/vector_bundle/algebras.py @@ -0,0 +1,559 @@ +r""" + Implementations of algorithms for associative algebras. Many of those + appear in the reference of Sage but the class is somewhat buggy for + the moment. + + A lot of it is very hacky but I just want to get something that works + until the Sage implementation can be fixed. + + It's also pretty unoptimized, there is a lot of room for improvement + if it is too slow. +""" +########################################################################### +# Copyright (C) 2024 Mickaël Montessinos (mickael.montessinos@mif.vu.lt),# +# # +# Distributed under the terms of the GNU General Public License (GPL) # +# either version 3, or (at your option) any later version # +# # +# http://www.gnu.org/licenses/ # +########################################################################### +import itertools +from sage.misc.cachefunc import cached_function +from sage.categories.algebras import Algebras +from sage.misc.misc_c import prod +from sage.algebras.all import FiniteDimensionalAlgebra +from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing +from sage.arith.misc import CRT_list +from sage.matrix.matrix_space import MatrixSpace +from sage.matrix.constructor import matrix +from sage.matrix.special import block_matrix +from sage.modules.free_module_element import vector + + +def subalgebra(A, basis, category=None): + r""" + The main issue with the algebra class in Sage is that the submodule method + used to generate subalgebras raises an error. We implement this function + instead which generates structure constants for a subalgebra. + + This is not very efficient since we have to recompute the multiplication + table instead of using multiplication in A, but it works + + An advantage is that the identity of the subalgebra needs not be the + identity of A, so we may represent two-sided ideals as algebras as well. + + INPUT: + + -``A`` -- A finite dimensional algebra with basis + + -``basis`` -- A free family in A generating a subalgebra + + OUTPUT: + + -``S`` -- The subalgebra of ``A`` generated by basis + + -``to_S`` -- A left inverse to ``from_S`` (raises an error if the element + is not in S) + + -``from_S`` -- The inclusion map from S to A + """ + k = A.base_ring() + mat = matrix([e.vector() for e in basis]).transpose() + tables = [matrix([mat.solve_right((f*e).vector()) for f in basis]) + for e in basis] + if category is None: + category = A.category() + S = FiniteDimensionalAlgebra(k, tables, + assume_associative = True, + category = category) + to_S = lambda a : S(mat.solve_right(a.vector())) + from_S = lambda s : A(mat * s.vector()) + return S, to_S, from_S + + +def in_column_space(mat, vec): + r""" + Checks if vec is in the space spanned by the columns of mat. + + Not very pretty but does the job + + EXAMPLES :: + + sage: from vector_bundle.algebras import in_column_space + sage: mat = matrix(QQ,2,1,[1,2]) + sage: in_column_space(mat, vector([2,4])) + True + sage: in_column_space(mat, vector([0,1])) + False + """ + try: + mat.solve_right(vec) + except ValueError: + return False + return True + + +def clear_zero_columns(mat): + return mat.matrix_from_columns([i for i in range(mat.ncols()) + if not mat[:,i].is_zero()]) + + +def subalgebra_from_gens(A, gens, category=None): + mat = matrix([gen.vector() for gen in gens]).echelon_form() + basis = [A(row) for row in mat.rows() if not row.is_zero()] + return subalgebra(A, basis, category=category) + + +@cached_function +def radical(A): + return subalgebra(A, A.radical_basis()) + + +@cached_function +def center(A): + return subalgebra(A, A.center_basis(), category=A.category().Commutative()) + + +def wedderburn_malcev_basis(A): + r""" + Compute a basis of A of the form ``[r1,r2,...,rk,s1,...,sm]`` such that + ``[r1,...,rk]`` is a basis of the radical and ``[s1,...,sm]`` is a basis + of a semisimple subalgebra of ``A``. + + EXAMPLES :: + + sage: from vector_bundle import VectorBundle + sage: from vector_bundle.algebras import wedderburn_malcev_basis + sage: F. = FunctionField(GF(7)) + sage: L1 = VectorBundle(F, 3 * x.zeros()[0].divisor()) + sage: L2 = VectorBundle(F, -2 * x.poles()[0].divisor()) + sage: V = L1.direct_sum(L2) + sage: T = matrix(F,2,2,[5*x^2 + 3*x + 2, x^2 + 4*x + 4, 4*x^2 + x + 1, 6*x^2 + 4*x + 5]) + sage: V = V.apply_isomorphism(T) + sage: End = V.end() + sage: A, to_A, from_A = End.global_algebra() + sage: basis = wedderburn_malcev_basis(A) + sage: [T^-1 * from_A(b) * T for b in basis] + [ + [0 0] [1 0] + [0 1], [0 0] + ] + """ + if not A.radical_basis(): + return A.basis() + #A is not semi-simple + k = A.base_ring() + n = A.dimension() + basis = A.basis() + #Compute bases of the powers of the radical. + gens = [A.radical_basis()] + while gens[-1]: + gens.append([a*b for a in gens[0] + for b in gens[-1] if not (a*b).is_zero()]) + gens = gens[:-1] + gen_vecs = [[g.vector() for g in gen] for gen in gens] + #Organise them into a hierchical basis + mat = clear_zero_columns(matrix(gen_vecs[-1]).echelon_form().transpose()) + rad_power_steps = [n, n - mat.ncols()] + for i in range(len(gens)-2, -1, -1): + new_cols = matrix([c for c in gen_vecs[i] + if not in_column_space(mat, c)]) + new_cols = clear_zero_columns(new_cols.echelon_form().transpose()) + mat = block_matrix([[new_cols, mat]]) + rad_power_steps.append(n - mat.ncols()) + rad_power_steps.append(0) + rad_power_steps.reverse() + #Now, add a complement which we'll turn into an algebra + for i in range(n): + col = matrix(k,n,1,[1 if j == i else 0 for j in range(n)]) + if not in_column_space(mat, col): + mat = block_matrix([[col, mat]]) + if mat.ncols() == n: + break + r = rad_power_steps[1] + imat = mat**-1 + #We may now begin running the algorithm in earnest. + basis = [A(c) for c in mat.columns()] + table = [imat * c.left_matrix().transpose() * mat for c in basis[:r]] + for first, stop in zip(rad_power_steps[:-1], rad_power_steps[1:]): + b_pairs = list(itertools.product(enumerate(basis[:r]), repeat=2)) + d_tuples = [[0]*i + [d] + [0]*(r-1-i) + for i in range(r) for d in basis[first:stop]] + #(first, stop) marks the positions of the vectors generating V_i. + eq = matrix([sum([list(mat.solve_right( + (bi*ds[j] + + ds[i]*bj + - sum([table[i][s,j]*d + for (s, d) in enumerate(ds)])).vector())[first:stop]) + for (i, bi),(j, bj) in b_pairs], []) + for ds in d_tuples]) + target = vector(sum([list(mat.solve_right( + (sum([table[i][s,j]*b for s, b in enumerate(basis[:r])]) + - bi*bj).vector())[first:stop]) + for (i, bi), (j, bj) in b_pairs], [])) + sol = eq.solve_left(target) + sol_chunks = [sol[i: i + stop - first] + for i in range(0, len(sol), stop - first)] + d_sols = [A(mat[:, first: stop] * d) for d in sol_chunks] + for (i, d) in enumerate(d_sols): + basis[i] += d + return basis[:r] + + +@cached_function +def wedderburn_malcev_complement(A): + category = A.category().Semisimple() + return subalgebra(A, wedderburn_malcev_basis(A), category=category) + + +def idems_from_element(e): + factorization = e.minimal_polynomial().factor() + factors = [fa[0] for fa in factorization] + n = len(factors) + hs = [CRT_list([0]*i + [1] + [0]*(n-i-1), factors) for i in range(n)] + return [h(e) for h in hs] + + +def splitting_from_idems(C, idems): + return [subalgebra_from_gens( + C, + [idem*c*idem for c in C.basis()]) + for idem in idems] + + +def wedderburn_splitting_large_field(C): + r""" + Assume that C is semi-simple and commutative, and that the base field of C + is finite but larger that 2*C.dimension() + """ + k = C.base_ring() + d = C.dimension() + search_set = [k(i) for i in range(2*C.dimension())] + while True: + a = C(vector([k.random_element() for _ in range(d)])) + f = a.minimal_polynomial() + if f.degree() == d: + break + splits = splitting_from_idems(C, idems_from_element(a)) + return [(split[0], [split[2](b) for b in split[0].basis()]) + for split in splits] + +def wedderburn_splitting_small_field(C): + r""" + Assume that C is semi-simple and commutative, and that the base field of C + is finite. + """ + k = C.base_ring() + test = matrix([(a**k.cardinality()- a).vector() for a in C.basis()]) + basis = test.left_kernel() + if len(basis) == 1: + return [(C, C.basis())] + one = C.one() + e = basis[0] + mat = matrix([one.vector(), e]) + if mat.rank() == 1: + e = basis[1] + partial_split = splitting_from_idems(C,idems_from_element(C(e))) + splits_of_splits = [wedderburn_splitting_comm(S[0]) for S in partial_split] + return sum([[(split[0], + [S[2](b) for b in split[1]]) + for split in splits] + for S, splits in zip(partial_split, splits_of_splits)], []) + + +def wedderburn_splitting_comm(C): + r""" + Assume that C is semi-simple commutative. Return an isomorphism from A to a direct + sum of simple algebras. + """ + n = C.dimension() + if 2*n < C.base_ring().cardinality(): + return wedderburn_splitting_large_field(C) + return wedderburn_splitting_small_field(C) + + +def wedderburn_splitting(A): + r""" + Assume that A is semi-simple. Return a list of simple algebras with maps + to and from A. A is isomorphic to the direct sum of these algebras. + + EXAMPLES :: + + sage: from vector_bundle import VectorBundle + sage: from vector_bundle.algebras import (wedderburn_malcev_complement, + ....: wedderburn_splitting) + sage: F. = FunctionField(GF(7)) + sage: L1 = VectorBundle(F, 3 * x.zeros()[0].divisor()) + sage: L2 = VectorBundle(F, -2 * x.poles()[0].divisor()) + sage: V = L1.direct_sum(L2) + sage: T = matrix(F,2,2,[5*x^2 + 3*x + 2, + ....: x^2 + 4*x + 4, 4*x^2 + x + 1, + ....: 6*x^2 + 4*x + 5]) + sage: V = V.apply_isomorphism(T) + sage: End = V.end() + sage: A, to_A, from_A = End.global_algebra() + sage: S, to_S, from_S = wedderburn_malcev_complement(A) + sage: simples = wedderburn_splitting(S) + sage: idems = [T^-1 * from_A(from_S(s[2](s[0].one()))) * T for s in simples] + sage: len(idems) + 2 + sage: matrix(GF(7),[[1, 0], [0, 0]]) in idems + True + sage: matrix(GF(7),[[0, 0], [0, 1]]) in idems + True + """ + C, to_C, from_C = center(A) + split = wedderburn_splitting_comm(C) + return [subalgebra_from_gens(A, [a * from_C(b) + for a in A.basis() for b in s[1]]) + for s in split] + + +def subfield_as_field(A, basis): + r""" + Return a bijection from the center of A to the corresponding extension + of the base field of A + """ + n = len(basis) + if n == 1: + to_k = lambda a : matrix([A.one().vector()]).solve_left(a.vector())[0] + from_k = lambda a : a*A.one() + return A.base_ring(), to_k, from_k + k = A.base_ring() + q = k.cardinality() + mat = matrix([(c**(q**(n-1)) - c).vector() for c in basis]) + ker = mat.left_kernel().tranpose() + b = [c for c in basis() if not in_column_space(ker, c)][0] + f = b.minimal_polynomial() + K = k.extension(f) + mat = matrix([(b**n).vector() for n in range(f.degree() - 1)]) + def to_K(a): + sol = mat.solve_right(a.vector()) + return K(sol) + from_K = lambda a: A(mat*vector(a.list())) + return K, to_K, from_K + + +def min_poly_over_subfield(A, basis, e): + K, to_K, from_K = subfield_as_field(A, basis) + mat = matrix([c.vector() for c in basis]).transpose() + n = 0 + acc = e + while(True): + try: + sol = mat.solve_right(acc.vector()) + break + except ValueError: + n +=1 + mat = block_matrix([[ + mat, + matrix([(c*acc).vector() for c in basis]).transpose()]]) + acc *= e + d = len(basis) + coeffs = [-to_K(sum([c*v for c, v in zip(sol[i:i+d], basis)])) + for i in range(0, len(sol), d)] + [1] + return PolynomialRing(K, 't')(coeffs), mat + + +def zero_div_from_min_poly(a, f): + facto = f.factor() + if len(facto) > 1: + return((facto[0][0]**facto[0][1])(a), + prod([(fa[0]**fa[1])(a) for fa in facto[1:]])) + if facto[0][1] > 1: + return((facto[0][0])(a), (facto[0][0]**(facto[0][1]-1))(a)) + return None + + +def solve_norm_equation(K, L, a): + r""" + This is not a general norm equation solver. It works only in the cases + relevant for zero_div! + """ + d, r = L.degree().quorem(K.degree()) + assert r == 0 + if d % 2 == 1: + g = UnivariatePolynomialRing(L)([-a] + [0]*(d-1) + [1]) + facto = g.factor() + assert all([f[0].degree() == 1 for f in facto]) + h = facto[0][0] + return -h.coefficient(0)/h.coefficient(1) + if d == 2: + b = sqrt(L(-a)) + if b not in K: + return b + q = K.cardinality() + f = UnivariatePolynomialRing([1] + [0]*1 + [1]) + a = f.roots()[0][0] + return a*b + raise NotImplementedError + + +def zero_div(A): + r""" Return a pair of zero divisors in A or conclude that A is a field + extension of its base field (which, as always, is assumed finite) + """ + rad = A.radical_basis() + if rad: + a = rad[0] + n = [a**n for n,_ in enumerate(rad)].index(0) + return a, a**(n-1) + #A is semisimple + split = wedderburn_splitting(A) + if len(split) > 1: + s0, to_s0, from_s0 = split[0] + s1, to_s1, from_s1 = split[1] + return (from_s0(s0.one()), from_s1(s1.one())) + #A is simple + center_mat = matrix([c.vector() for c in A.center_basis()]).transpose() + q = A.base_ring().cardinality()**center_mat.ncols() + if A.is_commutative(): + return None, None + #A is not a field + b = [a for a in A.basis() if not in_column_space(center_mat, a.vector())][0] + f, Cb_mat = min_poly_over_subfield(A, A.center_basis(), b) + t = zero_div_from_min_poly(b, f) + if t is not None: + return t + #C(b) is a field + if f.degree() % 2 == 0 and f.degree() > 2: + Cb_basis = [A(c) for c in Cb_mat.columns()] + eq = matrix([(a**(q**2) - a).vector() for a in Cb_basis]) + ker = eq.left_kernel() + b = [sum([v[i]*b for i, b in enumerate(Cb_basis)]) for v in ker + if any(v[q:])] + #C(b) has rank odd or 2 over C + eq = matrix([(b*c - c*b**q).vector() for c in A.basis()]) + c = A(eq.left_kernel()[0]) + fc, _ = min_poly_over_subfield(A, A.center_basis(), c) + t = zero_div_from_min_poly(c, f) + if t is not None: + return t + #c is invertible and c^-1bc = b^q + fb, mat_Ab = min_poly_over_subfield(A, A.center_basis(), b) + basis_Ab = [A(v) for v in mat_Ab.columns()] + _, mat_Abc = min_poly_over_subfield(A, basis_Ab, c) + if mat_Abc.ncols() < A.dimension(): + basis = [A(v) for v in mat_Abc.columns()] + S, to_S, from_S = subalgebra(A, basis) + z1, z2 = zero_divisor(S) + return from_S(z1), from_S(z2) + #The subalgebra generated by the center, b and c is equal to A + n = fc.degree() + assert all([fc.coefficient(i) == 0 for i in range(1,n)]) + alpha = -fc.coefficient(0) + K, to_K, from_K = subfield_as_field(A, A.center_basis()) + L, to_L, from_L = subfield_as_field(A, basis_Ab) + d = from_L(solve_norm_equation(K, L, 1/alpha)) + x = 1 - c*d + y = sum([c**i * d**sum([q**j for j in range(i)]) for i in range(n)]) + return x, y + + +def idempotent_from_zero_div(A, z): + if z is None: + return A.one(), 1 + mat = matrix([(a*z).vector() for a in A.basis()]).echelon_form() + ideal_basis = [A(r) for r in mat.rows() if not r.is_zero()] + eq = matrix([sum([list((a*e).vector()) for a in ideal_basis], []) + for e in ideal_basis]) + target = vector(sum([list(a.vector()) for a in ideal_basis], [])) + e = sum([c*v for c, v in zip(eq.solve_left(target), ideal_basis)]) + n = (A.dimension() // len(A.center_basis())).sqrt() + r = len(ideal_basis) // (len(A.center_basis())*n) + if r > n/2: + return 1 - e, n - r + return e, r + + +def rank_one_idempotent(A): + r = len(A.center_basis()) + n = (A.dimension() // r).sqrt() + z1, _ = zero_div(A) + e, d = idempotent_from_zero_div(A, z1) + if d == 1: + return e + else: + S, to_S, from_S = subalgebra_from_gens( + A, + [a*e*b for a in A.basis() for b in A.basis()]) + e = rank_one_idempotent(S) + return from_S(e) + + +class SplitAlgebra: + r""" + Class for a splitting of an algebra. Contains the isomorphic matrix algebra + and maps to and from it + """ + def __init__(self, A, to_A=None, from_A = None): + z = rank_one_idempotent(A) + self.A = A + self._to_A = to_A + self._from_A = from_A + ideal_gens = [a*z for a in A.basis()] + center_basis = A.center_basis() + self._r = len(center_basis) + basis_over_center = [] + mat = matrix(A.base_ring(), A.dimension(), 0, []) + for g in ideal_gens: + if not in_column_space(mat, g.vector()): + basis_over_center.append(g) + mat = block_matrix([[ + mat, + matrix([(f*g).vector() for f in center_basis]).transpose()]]) + self._K, self._to_K, self._from_K = subfield_as_field(A, center_basis) + self._n = len(basis_over_center) + self.M = MatrixSpace(self._K, self._n) + self._eq = None + self._ideal_matrix = mat + self._basis_over_center = basis_over_center + + def _coords_over_center(self, a): + sol = self._ideal_matrix.solve_right(a.vector()) + coefs = [sum([c*v for c, v in zip(sol[i: i+self._r], self.A.center_basis())]) + for i in range(0, self._n, self._r)] + return vector(self._K, [self._to_K(c) for c in coefs]) + + def to_M(self, a): + if self._to_A is not None: + a = self._to_A(a) + return matrix( self._K, [self._coords_over_center(a*v) + for v in self._basis_over_center]).transpose() + + def from_M(self, mat): + if self._eq is None: + self._eq = matrix([sum([list((self._coords_over_center(a*b))) + for b in self._basis_over_center],[]) + for a in self.A.basis()]) + target = vector(mat.transpose().list()) + sol = self._eq.solve_left(target) + res = sum([self._from_K(s)*a for s, a in zip(sol, self.A.basis())]) + if self._from_A is not None: + res = self._from_A(res) + return res + + +def full_split(A): + r""" + A is any associative algebra over a finite field. + Returns the wedderburn_malcev complement S of A, and then a list of + splittings of all the simple factors of S. + """ + S = wedderburn_malcev_complement(A) + ss_split = wedderburn_splitting(S[0]) + splits = [SplitAlgebra(ss[0], ss[1], ss[2]) for ss in ss_split] + return S, splits + + +def random_central_simple_algebra(k, n): + Mat = MatrixSpace(k, n) + while True: + Isom = matrix(k, n**2, n**2, [k.random_element() for _ in range(n**4)]) + if Isom.is_unit(): + break + basis = [Mat(c) for c in Isom.columns()] + table = [matrix([Isom.solve_right(vector(bi * bj)) for bi in basis]) for bj in basis] + category = Algebras(k).FiniteDimensional().WithBasis().Semisimple() + A = FiniteDimensionalAlgebra(k, table, category=category) + return A, Isom diff --git a/vector_bundle/constructions.py b/vector_bundle/constructions.py index 3279d54..652cf07 100644 --- a/vector_bundle/constructions.py +++ b/vector_bundle/constructions.py @@ -57,7 +57,7 @@ def canonical_bundle(K): sage: L == trivial_bundle(K) True """ - pi,_ = function_field_utility.safe_uniformizer(K) + pi = function_field_utility.safe_uniformizers(K)[0] return VectorBundle(K, pi.differential().divisor()) def _euclid(a,b): @@ -94,7 +94,10 @@ def atiyah_bundle(field, rank, degree, base=None): sage: F. = FunctionField(GF(11)) sage: R. = F[] sage: K. = F.extension(y^2 - x^3 - x) - sage: base = VectorBundle(K, K.places_finite()[0].divisor()) + sage: base = VectorBundle( + ....: K, + ....: K.places_finite()[0].divisor() + ....: - K.places_infinite()[0].divisor()) sage: E = atiyah_bundle(K, 5, 3, base) sage: E.rank() 5 @@ -120,22 +123,28 @@ def atiyah_bundle(field, rank, degree, base=None): if degree < 0: return atiyah_bundle(field, rank, -degree, base).dual() divisor = field.places_infinite()[0].divisor() - gcd = _euclid(rank, degree) - plan = [(i % 2,q) for i,(_, _, q, _) in enumerate(gcd)] - a, b = plan[-1] - plan[-1] = (a, b - 1) - starting_rank = gcd[-1][1] + if degree == 0: + plan = [] + starting_rank = rank + else: + gcd = _euclid(rank, degree) + plan = [(i % 2,q) for i,(_, _, q, _) in enumerate(gcd)] + a, b = plan[-1] + plan[-1] = (a, b - 1) + starting_rank = gcd[-1][1] result = trivial_bundle(field) line_bundle = VectorBundle(field, divisor) for _ in range(starting_rank - 1): result = result.extension_by_global_sections() - result = result.tensor_product(line_bundle) - for move, reps in reversed(plan): + result = result.tensor_product(base) + if degree > 0: + result = result.tensor_product(line_bundle) + for op, reps in reversed(plan): for _ in range(reps): - if move == 0: - result = result.extension_by_global_sections() - else: + if op: result = result.tensor_product(line_bundle) + else: + result = result.extension_by_global_sections() return result diff --git a/vector_bundle/ext_group.py b/vector_bundle/ext_group.py index 0798b7c..0b85f0f 100644 --- a/vector_bundle/ext_group.py +++ b/vector_bundle/ext_group.py @@ -135,11 +135,30 @@ class ExtGroup(SageObject): """ return self._right + def dual_bundle(self): + r""" + Return the dual bundle of the `\mathcal{Ext}^1` bundle. This bundle + is `\omega \otimes \mathrm{right}^\vee \otimes \mathrm{left}`, where + `\omega` is a canonical line bundle of `K`. + + EXAMPLES :: + + sage: from vector_bundle import VectorBundle + sage: F. = FunctionField(GF(11)) + sage: R. = F[] + sage: K. = F.extension(y^2 - x^5 - 1) + sage: ksi = VectorBundle(K, K.places_finite()[0].divisor()) + sage: ext = ksi.extension_group(ksi.dual()); ext.dual_bundle() + Homomorphism bundle from Vector bundle of rank 1 over Function + field in y defined by y^2 + 10*x^5 + 10 to Vector bundle of rank 1 + over Function field in y defined by y^2 + 10*x^5 + 10 + """ + return self._ext_dual_bundle + def dual_basis(self): r""" - Return a basis of the dual of the `Ext^1` group. This is a basis of - `H^0(\omega \otimes \mathrm{right}^\vee \otimes \mathrm{left})`, where - `\omega` is a canonical line bundle. + This is the same as ``self.dual_bundle().h0()`` + The output basis is dual to ``self.basis()`` under Serre duality. EXAMPLES :: diff --git a/vector_bundle/function_field_utility.py b/vector_bundle/function_field_utility.py index 3089556..b4a48e5 100644 --- a/vector_bundle/function_field_utility.py +++ b/vector_bundle/function_field_utility.py @@ -1,3 +1,16 @@ +r""" + Implementations of all sort of algorithms for function field that are not + yet implemented in Sage. + + A lot of code for the FunctionFieldCompletionCustom class comes directly + from the sage source, and was written by Kwankyu Lee. + + REFERENCE: +.. [Coh00] H. Cohen + Advanced topics in computational number theory + Springer + 2000 +""" ########################################################################### # Copyright (C) 2024 Mickaël Montessinos (mickael.montessinos@mif.vu.lt),# # # @@ -8,7 +21,9 @@ ########################################################################### from sage.matrix.constructor import matrix -from sage.rings.infinity import Infinity +from sage.categories.map import Map +from sage.modules.free_module_element import vector +from sage.rings.infinity import infinity from copy import copy from sage.misc.cachefunc import cached_function from sage.misc.misc_c import prod @@ -19,6 +34,8 @@ from sage.rings.function_field.function_field_rational\ import RationalFunctionField from sage.rings.function_field.order_rational\ import FunctionFieldMaximalOrderInfinite_rational +from sage.rings.function_field.order import FunctionFieldOrderInfinite +from sage.rings.function_field.ideal import FunctionFieldIdealInfinite @cached_function def all_infinite_places(K): @@ -49,7 +66,7 @@ def infinite_valuation(a): 1 """ if a == 0: - return Infinity + return infinity return a.denominator().degree() - a.numerator().degree() @@ -71,7 +88,7 @@ def infinite_mod(a,i): def infinite_integral_matrix(mat): r""" - Return an matrix with coefficient in the infinite maximal order and its denominator. + Return a matrix with coefficient in the infinite maximal order and its denominator. INPUT: @@ -116,11 +133,16 @@ def infinite_hermite_form(mat,include_zero_cols=True,transformation=False): sage: K. = FunctionField(GF(3)) sage: R = K.maximal_order_infinite() sage: mat = matrix(R,[[1, x**-1, x**-2, (x**3+1) / x**3], [(2*x+2) / (x**3+2), x**-2, (x**2+2) / (x**4+1), 1]]) - sage: H,T = infinite_hermite_form(mat,transformation=True); H + sage: H, T = infinite_hermite_form(mat,transformation=True); H [0 0 1 0] [0 0 0 1] sage: mat*T == H True + sage: H, T = infinite_hermite_form(mat, False, True); H + [1 0] + [0 1] + sage: mat*T == H + True TESTS: @@ -165,6 +187,7 @@ def infinite_hermite_form(mat,include_zero_cols=True,transformation=False): H *= E if not include_zero_cols: H = H[:,r-n:] + T = T[:,r-n:] if transformation: return H,T return H @@ -225,6 +248,9 @@ def infinite_order_xgcd(ideals): sage: all([a[i] in primes[i] for i in range(2)]) True """ + s = len(ideals) + non_zero_indices = [i for i, ideal in enumerate(ideals) if ideal != 0] + ideals = [ideal for i, ideal in enumerate(ideals) if ideal != 0] order_basis = ideals[0].ring().basis() if order_basis[0] != 1: raise ValueError('The first element of the basis of the order should' @@ -241,8 +267,12 @@ def infinite_order_xgcd(ideals): H,U = infinite_hermite_form(C, include_zero_cols=False, transformation=True) if not (H/den).is_one(): raise ValueError("The ideals should be coprime.") - v = U[:,-n].list() - return [sum([ideals_bases[i][j]*v[n*i+j] for j in range(n)]) for i in range(k)] + v = U[:,0].list() + coefs = [sum([ideals_bases[i][j]*v[n*i+j] for j in range(n)]) for i in range(k)] + res = [0]*s + for i, c in zip(non_zero_indices, coefs): + res[i] = c + return res def infinite_approximation(places,valuations,residues): @@ -272,55 +302,296 @@ def infinite_approximation(places,valuations,residues): @cached_function -def safe_uniformizer(K): +def safe_uniformizers(K): r""" Return a safe uniformizer and an infinite place of self._function_field A uniformizer is safe if its valuation at other infinite places is 0. EXAMPLES: - sage: from vector_bundle.function_field_utility import safe_uniformizer + sage: from vector_bundle.function_field_utility import safe_uniformizers sage: F. = FunctionField(GF(3)) sage: R. = F[] sage: K. = F.extension(y^2 - x**-5 - 1) sage: places = K.places_infinite() - sage: pi, place = safe_uniformizer(K); pi - ((2*x + 1)/x)*y + (2*x + 2)/x - sage: place == places[0] + sage: pis = safe_uniformizers(K) + sage: all([(pi.valuation(place) == 1 and i == j) or (pi.valuation(place) == 0 and i != j) for (i,pi) in enumerate(pis) for (j, place) in enumerate(places)]) True - sage: pi.valuation(place) - 1 - sage: pi.valuation(places[1]) - 0 - TESTS: - - sage: from vector_bundle.function_field_utility import all_infinite_places - sage: F. = FunctionField(GF(3)) - sage: R. = F[] - sage: K. = F.extension(y^2 + x + 2) - sage: places = K.places_infinite() - sage: pi, place = safe_uniformizer(K); pi - 1/x*y - sage: place == places[0] - True - sage: pi.valuation(place) - 1 - sage: R. = F[] - sage: K. = F.extension(y^4 + (2*x^2 + 2)/x^2) - sage: pi, _ = safe_uniformizer(K) - sage: [pi.valuation(place) for place in all_infinite_places(K)] - [1, 0, 0] - sage: safe_uniformizer(F) - (1/x, Place (1/x)) """ places = all_infinite_places(K) n = len(places) - return (infinite_approximation( + return [infinite_approximation( places, - [2] + ([1]*(n-1)), - [places[0].local_uniformizer()] + ([1]*(n-1))), - places[0]) + [2 if p == place else 1 for p in places], + [place.local_uniformizer() if p == place else 1 for p in places]) + for place in places] + + +class FunctionFieldCompletionCustom(Map): + """ + Completions on function fields. + + Allows for choice of uniformizer. + + INPUT: + + - ``field`` -- function field + + - ``place`` -- place of the function field + + - ``pi`` -- a local uniformizer at place + + - ``name`` -- string for the name of the series variable + + - ``prec`` -- positive integer; default precision + + - ``gen_name`` -- string; name of the generator of the residue + field; used only when place is non-rational + + EXAMPLES:: + + sage: from vector_bundle.function_field_utility import FunctionFieldCompletionCustom + sage: K. = FunctionField(GF(2)); _. = K[] + sage: L. = K.extension(Y^2 + Y + x + 1/x) + sage: p = L.places_finite()[0] + sage: m = FunctionFieldCompletionCustom(L,p) + sage: m + Completion map: + From: Function field in y defined by y^2 + y + (x^2 + 1)/x + To: Laurent Series Ring in s over Finite Field of size 2 + sage: m(x) + s^2 + s^3 + s^4 + s^5 + s^7 + s^8 + s^9 + s^10 + s^12 + s^13 + + s^15 + s^16 + s^17 + s^19 + O(s^22) + sage: m(y) + s^-1 + 1 + s^3 + s^5 + s^7 + s^9 + s^13 + s^15 + s^17 + O(s^19) + sage: m(x*y) == m(x) * m(y) + True + sage: m(x+y) == m(x) + m(y) + True + + The variable name of the series can be supplied. If the place is not + rational such that the residue field is a proper extension of the constant + field, you can also specify the generator name of the extension:: + + sage: p2 = L.places_finite(2)[0] + sage: p2 + Place (x^2 + x + 1, x*y + 1) + sage: m2 = FunctionFieldCompletionCustom(L, p2, name='t', gen_name='b') + sage: m2(x) + (b + 1) + t + t^2 + t^4 + t^8 + t^16 + O(t^20) + sage: m2(y) + b + b*t + b*t^3 + b*t^4 + (b + 1)*t^5 + (b + 1)*t^7 + b*t^9 + b*t^11 + + b*t^12 + b*t^13 + b*t^15 + b*t^16 + (b + 1)*t^17 + (b + 1)*t^19 + O(t^20) + + The choice of local uniformizer used for the expansion can be supplied. + + sage: from vector_bundle.function_field_utility import safe_uniformizers + sage: from vector_bundle.function_field_utility import all_infinite_places + sage: F. = FunctionField(GF(3)) + sage: R. = F[] + sage: K. = F.extension(y^2 - x**-5 - 1) + sage: pi = safe_uniformizers(K)[0] + sage: place = all_infinite_places(K)[0] + sage: f = 1 / (1-pi) + sage: m3 = FunctionFieldCompletionCustom(K, place, pi) + sage: m3(f) + 1 + s + s^2 + s^3 + s^4 + s^5 + s^6 + s^7 + s^8 + s^9 + s^10 + s^11 + + s^12 + s^13 + s^14 + s^15 + s^16 + s^17 + s^18 + s^19 + O(s^20) + """ + def __init__(self, field, place, pi=None, name=None, prec=None, gen_name=None): + """ + Initialize. + + EXAMPLES:: + + sage: # needs sage.rings.finite_rings sage.rings.function_field + sage: K. = FunctionField(GF(2)); _. = K[] + sage: L. = K.extension(Y^2 + Y + x + 1/x) + sage: p = L.places_finite()[0] + sage: m = L.completion(p) + sage: m + Completion map: + From: Function field in y defined by y^2 + y + (x^2 + 1)/x + To: Laurent Series Ring in s over Finite Field of size 2 + """ + if name is None: + name = 's' # default + + if gen_name is None: + gen_name = 'a' # default + + k, from_k, to_k = place.residue_field(name=gen_name) + + self._place = place + if pi is None: + self._pi = place.local_uniformizer() + else: + self._pi = pi + + self._gen_name = gen_name + + if prec is infinity: + from sage.rings.lazy_series_ring import LazyLaurentSeriesRing + codomain = LazyLaurentSeriesRing(k, name) + self._precision = infinity + else: # prec < infinity: + # if prec is None, the Laurent series ring provides default precision + from sage.rings.laurent_series_ring import LaurentSeriesRing + codomain = LaurentSeriesRing(k, name=name, default_prec=prec) + self._precision = codomain.default_prec() + + Map.__init__(self, field, codomain) + + def _repr_type(self) -> str: + """ + Return a string containing the type of the map. + + EXAMPLES:: + + sage: # needs sage.rings.finite_rings sage.rings.function_field + sage: K. = FunctionField(GF(2)); _. = K[] + sage: L. = K.extension(Y^2 + Y + x + 1/x) + sage: p = L.places_finite()[0] + sage: m = L.completion(p) + sage: m # indirect doctest + Completion map: + From: Function field in y defined by y^2 + y + (x^2 + 1)/x + To: Laurent Series Ring in s over Finite Field of size 2 + """ + return 'Completion' + + def _call_(self, f): + """ + Call the completion for f + + EXAMPLES:: + + sage: # needs sage.rings.finite_rings sage.rings.function_field + sage: K. = FunctionField(GF(2)); _. = K[] + sage: L. = K.extension(Y^2 + Y + x + 1/x) + sage: p = L.places_finite()[0] + sage: m = L.completion(p) + sage: m(y) + s^-1 + 1 + s^3 + s^5 + s^7 + s^9 + s^13 + s^15 + s^17 + O(s^19) + """ + if f.is_zero(): + return self.codomain().zero() + if self._precision is infinity: + return self._expand_lazy(f) + else: + return self._expand(f, prec=None) + + def _call_with_args(self, f, args, kwds): + """ + Call the completion with ``args`` and ``kwds``. + + EXAMPLES:: + + sage: # needs sage.rings.finite_rings sage.rings.function_field + sage: K. = FunctionField(GF(2)); _. = K[] + sage: L. = K.extension(Y^2 + Y + x + 1/x) + sage: p = L.places_finite()[0] + sage: m = L.completion(p) + sage: m(x+y, 10) # indirect doctest + s^-1 + 1 + s^2 + s^4 + s^8 + O(s^9) + """ + if f.is_zero(): + return self.codomain().zero() + if self._precision is infinity: + return self._expand_lazy(f, *args, **kwds) + else: + return self._expand(f, *args, **kwds) + + def _expand(self, f, prec=None): + """ + Return the Laurent series expansion of f with precision ``prec``. + + INPUT: + + - ``f`` -- element of the function field + + - ``prec`` -- positive integer; relative precision of the series + + EXAMPLES:: + + sage: # needs sage.rings.finite_rings sage.rings.function_field + sage: K. = FunctionField(GF(2)); _. = K[] + sage: L. = K.extension(Y^2 + Y + x + 1/x) + sage: p = L.places_finite()[0] + sage: m = L.completion(p) + sage: m(x, prec=20) # indirect doctest + s^2 + s^3 + s^4 + s^5 + s^7 + s^8 + s^9 + s^10 + s^12 + s^13 + s^15 + + s^16 + s^17 + s^19 + O(s^22) + """ + if prec is None: + prec = self._precision + + place = self._place + F = place.function_field() + der = F.higher_derivation() + + k, from_k, to_k = place.residue_field(name=self._gen_name) + sep = self._pi + + val = f.valuation(place) + e = f * sep**(-val) + + coeffs = [to_k(der._derive(e, i, sep)) for i in range(prec)] + return self.codomain()(coeffs, val).add_bigoh(prec + val) + + def _expand_lazy(self, f): + """ + Return the lazy Laurent series expansion of ``f``. + + INPUT: + + - ``f`` -- element of the function field + + EXAMPLES:: + + sage: # needs sage.rings.finite_rings sage.rings.function_field + sage: K. = FunctionField(GF(2)); _. = K[] + sage: L. = K.extension(Y^2 + Y + x + 1/x) + sage: p = L.places_finite()[0] + sage: m = L.completion(p, prec=infinity) + sage: e = m(x); e + s^2 + s^3 + s^4 + s^5 + s^7 + s^8 + ... + sage: e.coefficient(99) # indirect doctest + 0 + sage: e.coefficient(100) + 1 + """ + place = self._place + F = place.function_field() + der = F.higher_derivation() + + k, from_k, to_k = place.residue_field(name=self._gen_name) + sep = self._pi + + val = f.valuation(place) + e = f * sep**(-val) + + def coeff(s, n): + return to_k(der._derive(e, n - val, sep)) + + return self.codomain().series(coeff, valuation=val) + + def default_precision(self): + """ + Return the default precision. + + EXAMPLES:: + + sage: # needs sage.rings.finite_rings sage.rings.function_field + sage: K. = FunctionField(GF(2)); _. = K[] + sage: L. = K.extension(Y^2 + Y + x + 1/x) + sage: p = L.places_finite()[0] + sage: m = L.completion(p) + sage: m.default_precision() + 20 + """ + return self._precision def local_expansion(place,pi,f): r""" @@ -342,26 +613,7 @@ def local_expansion(place,pi,f): EXAMPLES: - sage: from vector_bundle.function_field_utility import local_expansion - sage: from vector_bundle.function_field_utility import safe_uniformizer - sage: F. = FunctionField(GF(3)) - sage: R. = F[] - sage: K. = F.extension(y^2 - x**-5 - 1) - sage: pi, place = safe_uniformizer(K) - sage: f = 1 / (1-pi) - sage: exp = local_expansion(place, pi, f) - sage: all([exp(i) == 1 for i in range(20)]) - True """ - if f == 0: - return lambda i : 0 - K = place.function_field() - der = K.higher_derivation() - k, _, to_k = place.residue_field() - val = f.valuation(place) - e = f * pi**(-val) - return lambda i : to_k(der._derive(e, i - val, pi)) if i >= val else 0 - def residue(place,pi,f): r""" Return the residue of constant répartition f at place with respect @@ -464,3 +716,310 @@ def smallest_norm_first(mat,i = 0,norms=[]): norms[i] = norms[j+i] norms[j+i] = n return norms + + +def finite_order_xgcd(left, right): + r""" + Compute `a \in left` and `b \in right` such that `a + b = 1` + + INPUT: + + -a: FunctionFieldIdeal_polymod: ideal of a finite maximal order + -b: FunctionFieldIdeal_polymod: ideal of the same maximal order + + ALGORITHM: + + [Coh00]_ Algorithm 1.3.2 + + EXAMPLES :: + + sage: from vector_bundle.function_field_utility import finite_order_xgcd + sage: F. = FunctionField(GF(7)) + sage: R. = F[] + sage: K. = F.extension(y^2 - x^3 - x) + sage: places = K.places_finite() + sage: left = places[0].prime_ideal() + sage: right = places[1].prime_ideal() + sage: a, b = finite_order_xgcd(left, right) + sage: a in left + True + sage: b in right + True + sage: a + b + 1 + """ + O = left.base_ring() + O_basis = O.basis() + m_left = left.hnf() + n = m_left.nrows() + m_right = right.hnf() + c = block_matrix([[m_left], [m_right]]) + h, u = c.hermite_form(False, True) + if not h.is_one(): + raise ValueError('The ideals left and right must be coprime.') + x = u[0,:n].list() + a = sum([c*e for c, e in zip(x,left.gens_over_base())]) + return a, 1-a + + +def euclidean_step(ideal_a, ideal_b, a, b, d=None): + r""" + Let `d = ideal_a a + ideal_b b`. Return `u \in ideal_a d^-1` and + `v \in ideal_b d^-1` such that au + bv = 1 + + ALGORITHM: + + [Coh00]_ Theorem 1.3.3 + + EXAMPLES :: + + sage: from vector_bundle.function_field_utility import euclidean_step + sage: F. = FunctionField(GF(7)) + sage: R. = F[] + sage: K. = F.extension(y^2 - x^3 - x) + sage: ideals = [P.prime_ideal() for P in K.places_finite()[:2]] + sage: a = x^2*y + 3 + sage: b = 2*y*x^5 + sage: d = a*ideals[0] + b*ideals[1] + sage: u, v = euclidean_step(ideals[0], ideals[1], a, b) + sage: u in ideals[0] * d^-1 + True + sage: v in ideals[1] * d^-1 + True + sage: a*u + b*v + 1 + """ + infinite = isinstance(ideal_a.base_ring(),FunctionFieldOrderInfinite) + if a == 0: + return 0, b^-1 + if b == 0: + return a^-1, 0 + if d == None: + d = a*ideal_a + b*ideal_b + I = a * ideal_a * d**-1 + J = b * ideal_b * d**-1 + #It would make sense to distinguish between finite and infinite order + #in a unified xgcd function, but the hnf form for infinite ideals + #needs to be refactored: we currently use opposite convention from + #that of sage. + if infinite: + s, t = infinite_order_xgcd([I,J]) + else: + s, t = finite_order_xgcd(I, J) + return s/a, t/b + + +def finite_integral_quotient(left, right): + return (left.numerator()*right.denominator()) // (left.denominator())*(right.numerator()) + + +def infinite_integral_quotient(left, right): + if left == 0: + return 0 + r = left / right + x = left.parent().gen() + return x**(r.denominator().degree() - r.numerator().degree()) + + +def hnf_reduction_mod_ideal(ideal, elem): + r""" + Reduce an element of a function field \(K\) modulo an ideal of a maximal + order of K + """ + if isinstance(ideal, FunctionFieldIdealInfinite): + quotient = infinite_integral_quotient + hnf = infinite_ideal_hnf(ideal) + else: + quotient = finite_integral_quotient + hnf = ideal.hnf() + n = hnf.ncols() + basis = ideal.base_ring().basis() + basis_matrix = matrix([e.list() for e in basis]).transpose()**-1 + y = basis_matrix * matrix(n,1,elem.list()) + for i in range(n - 1, -1, -1): + q = quotient(y[i,0],hnf[i,i]) + y -= q * hnf[:,i] + y = sum([y[i,0]*e for i, e in enumerate(basis)]) + return y + + +def pseudo_hermite_form(ideals, mat, include_zero_cols=True, transformation=False): + r""" + Return the hermite form of the pseudo-matrix ``(ideals, mat)`` with + coefficients in a function field and ideals in a maximal order. + + WARNING: + + Uses the opposite convention from sage for hermite forms, aligns with + Cohen's book instead. + + ALGORITHM: + + - Algorithm 1.4.7 from [Coh00]_ + + EXAMPLES :: + + sage: from vector_bundle.function_field_utility import pseudo_hermite_form + sage: F. = FunctionField(GF(7)) + sage: R. = F[] + sage: K. = F.extension(y^2 - x^3 - x) + sage: ideals = [P.prime_ideal() for P in K.places_finite()[:3]] + sage: mat = matrix(K, [[1, x, y],[2, x+1, y+1]]) + sage: h_ideals, h, u = pseudo_hermite_form(ideals, mat, transformation=True) + sage: h + [ 0 1 3*x^3 + 4*x] + [ 0 0 1] + sage: h == mat * u + True + sage: all([u[i,j] in ideals[i] * h_ideals[j]^-1 for i in range(3) for j in range(3)]) + True + sage: prod(ideals) == u.determinant() * prod(h_ideals) + True + """ + K = mat.base_ring() + k = mat.ncols() + n = mat.nrows() + U = identity_matrix(K,k) + h = copy(mat) + h_ideals = copy(ideals) + j = k-1 + for i in range(n-1, -1, -1): + #Check zero + if all([h[i,m] == 0 for m in range(j+1)]): + continue + m = [h[i,m] == 0 for m in range(j+1)].index(False) + h[:,m], h[:,j] = h[:,j], h[:,m] + U[:,m], U[:,j] = U[:,j], U[:,m] + h_ideals[m], h_ideals[j] = h_ideals[j], h_ideals[m] + #Put 1 on the main diagonal + a = h[i,j]**-1 + h_ideals[j] *= h[i,j] + h[:,j] *= a + U[:,j] *= a + for m in range(j-1,-1,-1): + if h[i,m] == 0: + continue + #Euclidean step + partial = h[i,m]*h_ideals[m] + h_ideals[j] + u, v = euclidean_step(h_ideals[m], h_ideals[j], + h[i,m], 1, partial) + U[:, m], U[:, j] = U[:, m] - h[i, m]*U[:, j], u*U[:, m] + v*U[:, j] + h[:, m], h[:, j] = h[:, m] - h[i, m]*h[:, j], u*h[:, m] + v*h[:, j] + h_ideals[m], h_ideals[j] = h_ideals[m] * h_ideals[j] * partial**-1, partial + #Row reduction step + for m in range(j+1,k): + ideal = h_ideals[m]**-1 * h_ideals[j] + q = h[i, m] - hnf_reduction_mod_ideal(ideal, h[i, m]) + U[:, m] -= q*U[:, j] + h[:, m] -= q*h[:, j] + j -=1 + if not include_zero_cols: + first_nonzero = [h[:,j].is_zero() for j in range(k)].index(False) + h = h[:,first_nonzero:] + h_ideals = h_ideals[first_nonzero:] + U = U[:,first_nonzero:] + if transformation: + return h_ideals, h, U + return h_ideals, h + +def hermite_form_infinite_polymod(mat, include_zero_cols=True, transformation=False): + r""" + Return the hermite normal form of mat. + + EXAMPLES :: + sage: from vector_bundle.function_field_utility import hermite_form_infinite_polymod + sage: from vector_bundle.function_field_utility import all_infinite_places + sage: F. = FunctionField(GF(7)) + sage: R. = F[] + sage: K. = F.extension(y^2 - x^3 - x) + sage: mat = matrix(K,[[1,x^-1,y^-1],[2,(x+1)^-1,(y+1)^-1]]) + sage: h, u = hermite_form_infinite_polymod(mat, transformation=True) + sage: h + [ 0 (4*x + 1)/(x^2 + x) (4*x^2 + 6)/x^2] + [ 0 0 1] + sage: h == mat * u + True + sage: O = K.maximal_order_infinite() + sage: all([c in O for c in u.list()]) + True + sage: all([u.determinant().valuation(place) == 0 for place in all_infinite_places(K)]) + True + """ + K = mat.base_ring() + O = K.maximal_order_infinite() + pis = safe_uniformizers(K) + places = all_infinite_places(K) + mins = [min([m.valuation(place) for m in mat.list()]) for place in places] + den = prod([pi**min(-m,0) for pi, m in zip(pis, mins)]) + h = den*mat + k = mat.ncols() + n = mat.nrows() + U = identity_matrix(K,k) + j = k-1 + for i in range(n-1, -1, -1): + if all([h[i,m] == 0 for m in range(j+1)]): + continue + min_vals = [min([h[i,m].valuation(place) for m in range(j+1)]) + for place in places] + gcd = prod([pi**m for pi, m in zip(pis, min_vals)]) + #put gcd on the diagonal + ideals = [O.ideal(h[i,m]/gcd) if not h[i,m].is_zero() else 0 + for m in range(j+1)] + coefs = infinite_order_xgcd(ideals) + ell = [c == 0 for c in coefs].index(False) + U[:,j], U[:, ell] = gcd*sum([(c/h[i,ell])*U[:,m] + for m,c in enumerate(coefs)]), U[:, j] + h[:,j], h[:, ell] = gcd*sum([(c/h[i,ell])*h[:,m] + for m,c in enumerate(coefs)]), h[:, j] + assert(all([u in O for u in U.list()])) + #eliminate coefficients left of diagonal + for m in range(j): + c = h[i,m]/h[i,j] + U[:,m] -= c*U[:,j] + h[:,m] -= c*h[:,j] + assert(all([u in O for u in U.list()])) + #reduce coefficients right of diagonal + for m in range(j+1,k): + ideal = O.ideal(h[i,j]) + q = (h[i, m] - hnf_reduction_mod_ideal(ideal,h[i, m]))/h[i,j] + U[:,m] -= q*U[:,j] + h[:,m] -= q*h[:,j] + j-=1 + h /= den + if not include_zero_cols: + first_nonzero = [h[:,j].is_zero() for j in range(k)].index(False) + h = h[:,first_nonzero:] + U = U[:,first_nonzero:] + if transformation: + return h, U + return h + + +def full_rank_matrix_in_completion(mat, place=None, pi=None): + r""" + Return a full rank matrix with coefficients in the constant base field. + + Its columns are concatenations of expansions of the coefficients in the + columns of mat. + """ + K = mat.base_ring() + k = K.constant_base_field() + s = mat.ncols() + r = mat.nrows() + if place is None: + if isinstance(K, RationalFunctionField): + place = K.gen().zeros()[0] + else: + place = K.get_place(1) + Kp = FunctionFieldCompletionCustom(K, place, pi, prec=infinity, name="pi", gen_name="b") + exps = [[Kp(c) for c in row] for row in mat] + vals = [min([mat[i,j].valuation(place) for j in range(s)]) + for i in range(r)] + ell = 0 + N = matrix(k,0,s) + while N.rank() < s: + for i in range(r): + row = [exps[i][j].coefficient(vals[i] + ell) for j in range(s)] + N = insert_row(N, (i+1)*(ell+1) - 1, row) + ell += 1 + return N, Kp, vals diff --git a/vector_bundle/hom_bundle.py b/vector_bundle/hom_bundle.py index fb4b12b..5d0816e 100644 --- a/vector_bundle/hom_bundle.py +++ b/vector_bundle/hom_bundle.py @@ -40,9 +40,14 @@ is the constant field of `K`:: # http://www.gnu.org/licenses/ # ########################################################################### +from sage.misc.cachefunc import cached_method from sage.matrix.constructor import matrix +from sage.matrix.special import identity_matrix from sage.modules.free_module_element import vector +from sage.categories.all import Algebras +from sage.algebras.all import FiniteDimensionalAlgebra from vector_bundle import VectorBundle +from vector_bundle import function_field_utility, algebras class HomBundle(VectorBundle): r""" @@ -231,6 +236,17 @@ class HomBundle(VectorBundle): If other is also a hom bundle, this is the hom bundle from ``self.codomain().tensor_product(other.domain())`` to ``self.domain().tensor_product(other.codomain()`` + + EXAMPLES:: + + sage: from vector_bundle import VectorBundle + sage: F. = FunctionField(GF(7)) + sage: L1 = VectorBundle(F, 2*x.zeros()[0].divisor()) + sage: L2 = VectorBundle(F, -3*x.poles()[0].divisor()) + sage: hom = L1.hom(L2) + sage: T = L1.tensor_product(L2) + sage: hom.hom(hom) == T.end() + True """ if isinstance(other, HomBundle): return self._codomain.tensor_product(other._domain)\ @@ -242,6 +258,17 @@ class HomBundle(VectorBundle): Return the tensor product of a hom bundle and a vector bundle. This is the same thing as ``self._domain.hom(self._codomain.tensor_product(other))`` + + EXAMPLES:: + + sage: from vector_bundle import VectorBundle + sage: F. = FunctionField(GF(7)) + sage: L1 = VectorBundle(F, 2*x.zeros()[0].divisor()) + sage: L2 = VectorBundle(F, -3*x.poles()[0].divisor()) + sage: hom = L1.hom(L2) + sage: hom2 = L1.hom(L2.tensor_product(L1)) + sage: hom.tensor_product(L1) == hom2 + True """ return self._domain.hom(self._codomain.tensor_product(other)) @@ -299,7 +326,246 @@ class HomBundle(VectorBundle): sage: all([all([a in O_infinity for a in (x**-2 * mat * g_infinite).list()]) for mat in h0]) True """ + if self._h0 is not None: + return self._h0 h0 = super().h0() - return [self._vector_to_matrix(v) for v in h0] + h0 = [self._vector_to_matrix(v) for v in h0] + self._h0 = h0 + return h0 + def coordinates_in_h0(self, mat): + r""" + Return the coordinates in the basis of ``self.h0()`` of the matrix + ``mat`` + EXAMPLES:: + sage: from vector_bundle import trivial_bundle, canonical_bundle + sage: F. = FunctionField(GF(3)) + sage: R. = F[] + sage: K. = F.extension(y^4 - x^-2 - 1) + sage: triv = trivial_bundle(K) + sage: can = canonical_bundle(K) + sage: V1 = triv.direct_sum(can) + sage: V2 = can.direct_sum(triv) + sage: hom = V1.hom(V2) + sage: hom.coordinates_in_h0(matrix([[0, 1], [1, 0]])) + (1, 1, 0, 0) + """ + if self._h0_matrix is None: + vec_h0 = [self._matrix_to_vector(m) for m in self.h0()] + basis_mat = matrix(self._function_field, vec_h0).transpose() + self._h0_matrix, self._h0_Kp, self._h0_vs =\ + function_field_utility.full_rank_matrix_in_completion(basis_mat) + vec_f = self._matrix_to_vector(mat) + res = VectorBundle.coordinates_in_h0(self, vec_f, check=False) + if mat == self.h0_from_vector(res): + return res + return None + + def image(self, f): + r""" + Return an image of global homomorphism ``f`` which is an element of + `H^0(\mathrm{self})`. + + That is, a vector bundle `V` together with an injective morphism of `V` + into ``self.codomain`` such that the image of `V` in ``self.codomain`` + is also the image of ``f``. + + EXAMPLES :: + + sage: from vector_bundle import trivial_bundle, canonical_bundle + sage: F. = FunctionField(GF(7)) + sage: R. = F[] + sage: K. = F.extension(y^2 - x^3 - x) + sage: triv = trivial_bundle(K) + sage: can = canonical_bundle(K) + sage: V1 = triv.direct_sum(can) + sage: V2 = can.direct_sum(triv) + sage: hom = V1.hom(V2) + sage: image, map = hom.image(matrix(K, [[0, 1], [0, 0]])) + sage: image.isomorphism_to(can) is not None + True + sage: image.hom(V2).coordinates_in_h0(map) + (1, 0) + """ + dom = self._domain + cod = self._codomain + ideals, C_fi = function_field_utility.pseudo_hermite_form( + dom._ideals, + f*dom._g_finite, + False) + C_inf = function_field_utility.hermite_form_infinite_polymod( + f*dom._g_infinite, + False) + g_fi = C_inf.solve_right(C_fi) + K = self._function_field + image = VectorBundle(K, ideals, g_fi, identity_matrix(K, g_fi.ncols())) + return image, C_inf + + def kernel(self, f): + r""" + Return a kernel of global homomorphism ``f`` which is an element of + `H^0(\mathrm{self})`. + + That is, a vector bundle `V` together with an injective morphism of `V` + into ``self.domain`` such that the image of `V` in ``self.domain`` + is the kernel of ``f``. + + EXAMPLES :: + + sage: from vector_bundle import trivial_bundle, canonical_bundle + sage: F. = FunctionField(GF(7)) + sage: R. = F[] + sage: K. = F.extension(y^2 - x^3 - x) + sage: triv = trivial_bundle(K) + sage: can = canonical_bundle(K) + sage: V1 = triv.direct_sum(can) + sage: V2 = can.direct_sum(triv) + sage: hom = V1.hom(V2) + sage: kernel, map = hom.kernel(matrix(K, [[0, 1], [0, 0]])) + sage: kernel.isomorphism_to(triv) is not None + True + sage: kernel.hom(V1).coordinates_in_h0(map) + (1, 0) + """ + dom = self._domain + cod = self._codomain + ideals, _, U_fi = function_field_utility.pseudo_hermite_form( + dom._ideals, + f*dom._g_finite, + transformation=True) + _, U_inf = function_field_utility.hermite_form_infinite_polymod( + f*dom._g_infinite, + transformation=True) + r = f.rank() + n = f.ncols() + ideals = ideals[:n-r] + C_fi = U_fi[:,:n-r] + C_inf = U_inf[:,:n-r] + g_fi = C_inf.solve_right(C_fi) + K = self._function_field + image = VectorBundle(K, ideals, g_fi, identity_matrix(K, g_fi.ncols())) + return image, C_inf + + def is_isomorphism(self, f): + r""" + Check if f is an isomorphism from ``self.domain()`` to + ``self.codomain()`` + + EXAMPLES:: + + sage: from vector_bundle import atiyah_bundle, canonical_bundle + sage: F. = FunctionField(GF(7)) + sage: R. = F[] + sage: K. = F.extension(y^2 - x^3 - x) + sage: V = atiyah_bundle(K, 2, 0) + sage: W = atiyah_bundle(K, 2, 0, canonical_bundle(K)) + sage: isom = x^2/(x^2 + 3) * identity_matrix(K, 2) + sage: V.hom(W).is_isomorphism(isom) + True + """ + dual = self.dual() + return self.is_in_h0(f) and dual.is_in_h0(f**-1) + + +class EndBundle(HomBundle): + r""" + Vector bundles representing endomorphism sheaves of vector bundles. + + EXAMPLES:: + + sage: from vector_bundle import trivial_bundle, canonical_bundle + sage: F. = FunctionField(GF(7)) + sage: R. = F[] + sage: K. = F.extension(y^4 - x^-2 - 1) + sage: triv = trivial_bundle(K) + sage: can = canonical_bundle(K) + sage: triv.direct_sum(can).end() + Endomorphism bundle of Vector bundle of rank 2 over Function field in y defined by y^4 + (6*x^2 + 6)/x^2 + """ + def __init__(self, bundle): + HomBundle.__init__(self, bundle, bundle) + self._A = None + + def __repr__(self): + return "Endomorphism bundle of %s" % (self._domain) + + def global_algebra(self): + r""" + Return ``self.h0()`` as a k-algebra in which computations may be done. + Also return maps to and from the algebra. + + EXAMPLES:: + + sage: from vector_bundle import trivial_bundle, canonical_bundle + sage: F. = FunctionField(GF(7)) + sage: R. = F[] + sage: K. = F.extension(y^4 - x^-2 - 1) + sage: triv = trivial_bundle(K) + sage: can = canonical_bundle(K) + sage: end = triv.direct_sum(can).end() + sage: A, to_A, from_A = end.global_algebra(); A + Finite-dimensional algebra of degree 4 over Finite Field of size 7 + sage: h0 = end.h0() + """ + if self._A is not None: + return self._A + k = self._function_field.constant_base_field() + category = Algebras(k).FiniteDimensional().WithBasis().Associative() + basis = self.h0() + tables = [matrix(k,[self.coordinates_in_h0(b*a) for b in basis]) + for a in basis] + algebra = FiniteDimensionalAlgebra( + k, + tables, + assume_associative=True, + category = category) + to_a = lambda mat : algebra(self.coordinates_in_h0(mat)) + from_a = lambda a : self.h0_from_vector(a.vector()) + self._A = (algebra, to_a, from_a) + return algebra, to_a, from_a + + @cached_method + def _global_algebra_split(self): + r""" + Return a splitting of ``self.global_algebra()`` + + OUTPUT: + + - ``factors`` -- List of the matrix algebras that are simple + factors of the semi-simple quotient of self + inverse. Come with maps to and from ``self.global_algebra()`` + + EXAMPLES:: + + sage: from vector_bundle import VectorBundle + sage: from sage.matrix.matrix_space import MatrixSpace + sage: F. = FunctionField(GF(7)) + sage: L1 = VectorBundle(F, 3 * x.zeros()[0].divisor()) + sage: L2 = VectorBundle(F, -2 * x.poles()[0].divisor())\ + ....: .direct_sum_repeat(2) + sage: V = L1.direct_sum(L2) + sage: T = matrix( + ....: F, 3, 3, + ....: [2*x, 4, x+1, 4*x, 4*x+6, 5*x+6, 5*x+2, 5*x+3, 2*x+6]) + sage: V = V.apply_isomorphism(T) + sage: End = V.end() + sage: _, _, from_A = End.global_algebra() + sage: S, factors = End._global_algebra_split() + sage: len(factors) + 2 + sage: deg_1_index = [f.M.nrows() for f in factors].index(1) + sage: deg_2_index = 1 - deg_1_index + sage: f = factors[deg_1_index] + sage: T^-1 * from_A(S[2](f.from_M(identity_matrix(GF(7),1)))) * T + [1 0 0] + [0 0 0] + [0 0 0] + sage: f = factors[deg_2_index] + sage: T^-1 * from_A(S[2](f.from_M(identity_matrix(GF(7),2)))) * T + [0 0 0] + [0 1 0] + [0 0 1] + """ + A, _, _ = self.global_algebra() + return algebras.full_split(A) diff --git a/vector_bundle/vector_bundle.py b/vector_bundle/vector_bundle.py index b18ee57..607a11f 100644 --- a/vector_bundle/vector_bundle.py +++ b/vector_bundle/vector_bundle.py @@ -76,10 +76,11 @@ from sage.structure.sage_object import SageObject from sage.misc.misc_c import prod from sage.arith.functions import lcm from sage.arith.misc import integer_ceil -from sage.functions.log import logb +from sage.functions.log import logb, log from sage.matrix.constructor import matrix -from sage.matrix.special import block_matrix, elementary_matrix\ - , identity_matrix +from sage.matrix.special import (block_matrix, elementary_matrix, + identity_matrix, diagonal_matrix, + zero_matrix, block_diagonal_matrix) from sage.matrix.matrix_space import MatrixSpace from sage.rings.function_field.ideal import FunctionFieldIdeal from sage.rings.function_field.function_field_rational\ @@ -138,12 +139,16 @@ class VectorBundle(SageObject): Vector bundle of rank 1 over Function field in y defined by y^2 + 2*x^3 + 2*x """ - def __init__(self,function_field, ideals,g_finite=None,g_infinite=None): + def __init__(self,function_field, ideals,g_finite=None,g_infinite=None, check=True): if g_finite is None or g_infinite is None: - self._line_bundle_from_divisor(function_field, ideals) + self._line_bundle_from_divisor(function_field, ideals, check=check) else: self._vector_bundle_from_data(function_field, ideals, - g_finite,g_infinite) + g_finite,g_infinite, check) + self._h0 = None + self._h0_matrix = None + self._h0_Kp = None + self._h0_vs = None def __hash__(self): return hash((tuple(self._ideals), @@ -164,13 +169,14 @@ class VectorBundle(SageObject): self._function_field, ) - def _line_bundle_from_divisor(self,function_field,divisor): + def _line_bundle_from_divisor(self,function_field,divisor, check=True): r""" Build a line bundle from a divisor """ - if not function_field == divisor.parent().function_field(): - raise ValueError('The divisor should be defined over the ' - + 'function field.') + if check: + if not function_field == divisor.parent().function_field(): + raise ValueError('The divisor should be defined over the ' + + 'function field.') self._function_field = function_field couples = divisor.list() finite_part = [c for c in couples if not c[0].is_infinite_place()] @@ -188,7 +194,7 @@ class VectorBundle(SageObject): self._g_infinite = matrix(function_field,[[pi]]) def _vector_bundle_from_data(self,function_field, - ideals,g_finite,g_infinite): + ideals,g_finite,g_infinite, check=True): r""" Construct a vector bundle from data. """ @@ -207,22 +213,23 @@ class VectorBundle(SageObject): g_finite.change_ring(function_field) g_infinite.change_ring(function_field) r = len(ideals) - if (g_finite.nrows() != r - or g_finite.ncols() != r - or g_infinite.nrows() != r - or g_infinite.ncols() != r): - raise ValueError('The length of the ideal list must equal the \ - size of the basis matrices') - if not g_finite.is_invertible() or not g_infinite.is_invertible(): - raise ValueError('The basis matrices must be invertible') - if not all([isinstance(I,FunctionFieldIdeal) - for I in ideals]): - raise TypeError('The second argument must be a list of \ - FunctionFieldIdeals.') - if not all([I.base_ring() == function_field.maximal_order() - for I in ideals]): - raise ValueError('All ideals must have the maximal order of\ - function_field as base ring.') + if check: + if (g_finite.nrows() != r + or g_finite.ncols() != r + or g_infinite.nrows() != r + or g_infinite.ncols() != r): + raise ValueError('The length of the ideal list must equal' + + ' the size of the basis matrices') + if not g_finite.is_invertible() or not g_infinite.is_invertible(): + raise ValueError('The basis matrices must be invertible') + if not all([isinstance(I,FunctionFieldIdeal) + for I in ideals]): + raise TypeError('The second argument must be a list of \ + FunctionFieldIdeals.') + if not all([I.base_ring() == function_field.maximal_order() + for I in ideals]): + raise ValueError('All ideals must have the maximal order of\ + function_field as base ring.') self._function_field = function_field self._ideals = ideals self._g_finite = g_finite @@ -376,10 +383,11 @@ class VectorBundle(SageObject): I = prod(self._ideals) determinant_finite = self._g_finite.determinant() determinant_infinite = self._g_infinite.determinant() - return VectorBundle(self._function_field - ,I - ,determinant_finite - ,determinant_infinite) + return VectorBundle(self._function_field, + I, + determinant_finite, + determinant_infinite, + check=False) def degree(self): r""" @@ -500,7 +508,8 @@ class VectorBundle(SageObject): [0 0], [x 0], [ 0 0], [0 1] ] """ - return self.hom(self) + from . import hom_bundle + return hom_bundle.EndBundle(self) def dual(self): r""" @@ -550,7 +559,8 @@ class VectorBundle(SageObject): ideals = self._ideals + other._ideals g_finite = block_matrix([[self._g_finite,0],[0,other._g_finite]]) g_infinite = block_matrix([[self._g_infinite,0],[0,other._g_infinite]]) - return VectorBundle(self._function_field, ideals,g_finite,g_infinite) + return VectorBundle(self._function_field, ideals, + g_finite, g_infinite, check=False) def _direct_sum_rec(self,acc,n): r""" @@ -608,7 +618,8 @@ class VectorBundle(SageObject): ideals = [I * J for I in self._ideals for J in other._ideals] g_finite = self._g_finite.tensor_product(other._g_finite) g_infinite = self._g_infinite.tensor_product(other._g_infinite) - return VectorBundle(self._function_field, ideals,g_finite,g_infinite) + return VectorBundle(self._function_field, ideals, + g_finite, g_infinite, check=False) def _tensor_power_aux(self,acc,n): r""" @@ -664,7 +675,8 @@ class VectorBundle(SageObject): """ O = K.maximal_order() ideals = [O.ideal(I.gens()) for I in self._ideals] - return VectorBundle(K, ideals,self._g_finite,self._g_infinite) + return VectorBundle(K, ideals,self._g_finite, + self._g_infinite, check=False) def restriction(self): r""" @@ -717,7 +729,8 @@ class VectorBundle(SageObject): for c in gen_infinite]) g_infinite = matrix([sum([a.list() for a in collumn],[]) for collumn in g_infinite]).transpose() - return VectorBundle(F, ideals,g_finite,g_infinite) + return VectorBundle(F, ideals,g_finite, + g_infinite, check=False) def _h0_rational(self): r""" @@ -799,7 +812,7 @@ class VectorBundle(SageObject): for i in range(self.rank()): for p in range(min([c.denominator().degree() - c.numerator().degree() - for c in mat[i, :].list()]) + 1): + for c in mat[i, :].list() if c != 0]) + 1): basis.append(self._g_infinite * vector(one.shift(p)*mat[i, :])) return basis @@ -851,6 +864,8 @@ class VectorBundle(SageObject): sage: len(L.h0()) 1 """ + if self._h0 is not None: + return self._h0 if isinstance(self._function_field,RationalFunctionField): #Compute restriction to normalize the coefficient ideals. return self.restriction()._h0_rational() @@ -863,6 +878,7 @@ class VectorBundle(SageObject): h0.append(vector([sum([y**j * v[i*deg + j] for j in range(deg)]) for i in range(self.rank())])) + self._h0 = h0 return h0 @cached_method @@ -922,17 +938,17 @@ class VectorBundle(SageObject): OUTPUT: - - ''res'' -- vector of elements of K such that the corresponding infinite répartition vectorcorresponds to form under Serre duality with respect to _safe_uniformizer(self._function_field).differential(). + - ''res'' -- vector of elements of K such that the corresponding infinite répartition vectorcorresponds to form under Serre duality with respect to ``safe_uniformizers(self._function_field)[0].differential()``. EXAMPLES :: - sage: from vector_bundle import VectorBundle + sage: from vector_bundle import trivial_bundle sage: F. = FunctionField(GF(3)) sage: R. = F[] - sage: K. = F.extension(y^4 - x**-2 - 1) - sage: trivial_ideal = K.maximal_order().ideal(1) - sage: g_finite = identity_matrix(K, 2) - sage: pi = K.places_infinite()[0].local_uniformizer() + sage: K. = F.extension(y^2 - x^3 - x) + sage: triv = trivial_bundle(K) + sage: triv.h1_element([1]) + [(x/(x^2 + 1))*y] """ K = self._function_field h1_dual, h1_dual_bundle = self.h1_dual() @@ -941,53 +957,90 @@ class VectorBundle(SageObject): form = [1] + [0] * (s-1) r = self.rank() places = function_field_utility.all_infinite_places(K) - pi_0, place_0 = function_field_utility.safe_uniformizer(K) - if not place_0 == places[0]: - raise ValueError('Something went wrong with the order of infinite' - + ' places of the function field') + pi_0 = function_field_utility.safe_uniformizers(K)[0] + place_0 = places[0] k,from_k,to_k = place_0.residue_field() form = vector([function_field_utility.invert_trace( k, K.constant_base_field(), c) for c in form]) dual_matrix = matrix([h1_dual_bundle._matrix_to_vector(mat) for mat in h1_dual]).transpose() zero_rows = [i for i, row in enumerate(dual_matrix) if row == 0] - exps = [[function_field_utility.local_expansion(place_0, - pi_0,dual_matrix[i, j]) - for j in range(s)] - for i in range(r)] - min_vals = [[min([dual_matrix[i, j].valuation(place) for j in range(s)]) - for i in range(r)] - for place in places] - ell = 0 - n_matrix = matrix(k,0,s,[]) - while n_matrix.rank() < s: - for i in range(r): - row = [exps[i][j](min_vals[0][i] + ell) for j in range(s)] - n_matrix = function_field_utility.insert_row( - n_matrix, (i+1)*(ell+1) - 1, row) - ell +=1 + n_matrix, _, _ = function_field_utility.full_rank_matrix_in_completion( + dual_matrix, + place_0, + pi_0) + ell = n_matrix.nrows() // dual_matrix.nrows() ell_pow = k.cardinality() ** integer_ceil(logb(ell,k.cardinality())) res = n_matrix.solve_left(form) + min_vals = [[min([c.valuation(place) for c in dual_matrix.list()])] + for place in places] pi = function_field_utility.infinite_approximation( places, - [1] + [integer_ceil(-min(min_val)/ell) for min_val in min_vals[1:]], + [1] + [integer_ceil(-min_val/ell) for min_val in min_vals[1:]], [1] + [0]*(len(places)-1))**ell_pow res = [function_field_utility.infinite_approximation( places, [1] + [0]*(len(places)-1), [a] + [0]*(len(places)-1))**ell_pow for a in res] + min_vals_0 = [0 if i in zero_rows + else min([dual_matrix[i,j].valuation(place_0) + for j in range(s)]) + for i in range(r)] return [pi - * pi_0**(-min_vals[0][i]-1) + * pi_0**(-min_vals_0[i]-1) * sum([pi_0**(-j) * res[i*ell + j] for j in range(ell)]) if not i in zero_rows else 0 for i in range(r)] @cached_method - def extension_group(self, other, precompute_basis=None): + def extension_group(self, other, precompute_basis=False): + r""" + Return the extension group of ``self`` by ``other``. + If ``precompute_basis`` is set to ``True``, a basis of the extension + group is precomputed. Extensions may be constructed without using + a precomputed basis, but each construction is a bit costlier this way. + You should precompute a basis if you are planning to compute + an number of extensions larger than the dimension of the ext group. + + EXAMPLES:: + + sage: from vector_bundle import trivial_bundle, canonical_bundle + sage: F. = FunctionField(GF(3)) + sage: R. = F[] + sage: K. = F.extension(y^2 - x^3 - x) + sage: trivial_bundle(K).extension_group(canonical_bundle(K)) + Extension group of Vector bundle of rank 1 over Function field in + y defined by y^2 + 2*x^3 + 2*x by Vector bundle of rank 1 over + Function field in y defined by y^2 + 2*x^3 + 2*x. + """ return ext_group.ExtGroup(self, other, precompute_basis) def non_trivial_extension(self, other): + r""" + Return any nontrivial extension of self by other. + + EXAMPLES:: + + sage: from vector_bundle import trivial_bundle, canonical_bundle + sage: F. = FunctionField(GF(3)) + sage: R. = F[] + sage: K. = F.extension(y^2 - x^3 - x) + sage: triv = trivial_bundle(K) + sage: can = canonical_bundle(K) + sage: V = triv.non_trivial_extension(can) + sage: V.rank() + 2 + sage: V.degree() + 0 + sage: V.h0() + [(1, 0)] + sage: V.end().h0() + [ + [0 1] [1 0] + [0 0], [0 1] + ] + """ ext_group = self.extension_group(other) return ext_group.extension() @@ -1000,14 +1053,6 @@ class VectorBundle(SageObject): constructions generalises to arbitrary genus if one replaces the trivial line bundle with a canonical line bundle. - WARNING: - - The implementation is hacky at the moment and relies on the - hypothesis that ``h1_dual`` returns a good basis of the space. - A more robust implementation would either use formal power series to - do linear algebra on the global sections or rely on the ``h0`` computation - using matrices in normal Popov form to have normalized output. - EXAMPLES :: sage: from vector_bundle import trivial_bundle, VectorBundle @@ -1018,7 +1063,7 @@ class VectorBundle(SageObject): sage: E = T.extension_by_global_sections() sage: E.rank() 2 - sage: E.hom(E).h0() + sage: E.end().h0() [ [0 1] [1 0] [0 0], [0 1] @@ -1043,33 +1088,366 @@ class VectorBundle(SageObject): ohm = constructions.canonical_bundle(self._function_field)\ .direct_sum_repeat(s) ext_group = self.extension_group(ohm) - ext_dual = ext_group.dual_basis() - form = [1 if any([vector(mat[:, i]) == v - for i,v in enumerate(h0)]) else 0 - for mat in ext_dual] + ext_dual = ext_group.dual_bundle() + canonical_ext = matrix(h0).transpose() + form = ext_dual.coordinates_in_h0(canonical_ext) return ext_group.extension(form) - def is_isomorphic_to(self, other): + def h0_from_vector(self, v): r""" - Checks whether self is isomorphic to other + Return an element of `H^0(\mathrm{self})` from a vector of coordinates + in the basis given by ``self.h0()`` - ALGORITHM: + EXAMPLES :: + + sage: from vector_bundle import VectorBundle + sage: F. = FunctionField(GF(7)) + sage: R. = F[] + sage: K. = F.extension(y^2 - x^3 - x) + sage: ideals = [P.prime_ideal() for P in K.places_finite()[:2]] + sage: g_finite = matrix([[1,1 / (x**5 + y)],[2, y]]) + sage: g_infinite = matrix([[x, 1], [2, y**3]]) + sage: V = VectorBundle(K, ideals, g_finite, g_infinite) + sage: h0 = V.h0() + sage: v = vector(list(range(6))) + sage: V.coordinates_in_h0(V.h0_from_vector(v)) + (0, 1, 2, 3, 4, 5) + """ + return sum([a*e for a, e in zip(v, self.h0())]) + + def coordinates_in_h0(self, f, check=True): + r""" + Return a vector of coordinates of ``f`` in the basis returned + by ``self.h0()`` - Computes the space of global homomorphisms and looks for - an invertible matrix. Can probably be improved. + If ``check`` is ``True``, it is check whether ``f`` actually lies in the + `H^0` space, and None is output if ``f`` is not in `H^0`. If ``check`` + is set to ``False`` the result may be garbage. + + EXAMPLES :: + + sage: from vector_bundle import VectorBundle + sage: F. = FunctionField(GF(7)) + sage: R. = F[] + sage: K. = F.extension(y^2 - x^3 - x) + sage: ideals = [P.prime_ideal() for P in K.places_finite()[:2]] + sage: g_finite = matrix([[1,1 / (x**5 + y)],[2, y]]) + sage: g_infinite = matrix([[x, 1], [2, y**3]]) + sage: V = VectorBundle(K, ideals, g_finite, g_infinite) + sage: h0 = V.h0() + sage: v = sum([i*e for i, e in enumerate(h0)]) + sage: V.coordinates_in_h0(v) + (0, 1, 2, 3, 4, 5) + """ + if self._h0_matrix is None: + mat = matrix(self._function_field, self.h0()).transpose() + self._h0_matrix, self._h0_Kp, self._h0_vs =\ + function_field_utility.full_rank_matrix_in_completion(mat) + ell = self._h0_matrix.nrows() // self.rank() + series = [self._h0_Kp(c) for c in f] + v = vector(sum([[s.coefficient(val + j) for j in range(ell)] + for s, val in zip(series, self._h0_vs)],[])) + res = self._h0_matrix.solve_right(v) + if not check or self.h0_from_vector(res) == f: + return res + return None + + def is_in_h0(self, v): + r""" + Check if vector ``v`` with coefficients in + ``self.function_field()`` lies in the `k`-vector space spanned + by the output of ``self.h0()`` EXAMPLES :: - sage: from vector_bundle import trivial_bundle, canonical_bundle - sage: F. = FunctionField(GF(3)) + sage: from vector_bundle import VectorBundle + sage: F. = FunctionField(GF(7)) sage: R. = F[] - sage: K. = F.extension(y^4 - x^-2 - 1) - sage: triv = trivial_bundle(K) - sage: can = canonical_bundle(K) - sage: triv.is_isomorphic_to(can) + sage: K. = F.extension(y^2 - x^3 - x) + sage: ideals = [P.prime_ideal() for P in K.places_finite()[:2]] + sage: g_finite = matrix([[1,1 / (x**5 + y)],[2, y]]) + sage: g_infinite = matrix([[x, 1], [2, y**3]]) + sage: V = VectorBundle(K, ideals, g_finite, g_infinite) + sage: V.is_in_h0(V.h0_from_vector(vector(list(range(6))))) + True + sage: V.is_in_h0(vector([x^i for i in range(6)])) + False + """ + if self.coordinates_in_h0(v) is None: + return False + return True + + def _isomorphism_to_large_field(self, other, tries=1): + r""" + Return an isomorphism from self to other if it exists and None otherwise. + + May fail to find an isomorphism with probability less than + `(\frac{s}{|k|})^\mathrm{tries}`, where k is the constant field and + s is ``len(self.hom(other).h0())``. This is only usefule if `k` has + cardinality larger than ``len(self.end().h0())``. + """ + Hom1 = self.hom(other) + hom1 = Hom1.h0() + hom2 = Hom1.dual().h0() + End = self.end() + end = End.h0() + s = len(hom1) + if s != len(hom2) or s != len(end): + return None + k = self.function_field().constant_base_field() + for _ in range(tries): + v = vector([k.random_element() for _ in range(s)]) + P = Hom1.h0_from_vector(v) + mat = matrix([End.coordinates_in_h0(P*Q) for Q in hom2]) + if mat.is_unit(): + return P + + def _isomorphism_indecomposable(self, other): + r""" + Return an isomorphism from self to other if it exists. + Assumes that self is indecomposable. + + TESTS:: + + sage: from vector_bundle import atiyah_bundle, canonical_bundle + sage: F. = FunctionField(GF(7)) + sage: R. = F[] + sage: K. = F.extension(y^2 - x^3 - x) + sage: V = atiyah_bundle(K, 2, 0) + sage: W = atiyah_bundle(K, 2, 0,canonical_bundle(K)) + sage: isom = V._isomorphism_indecomposable(W) + sage: isom is not None + True + sage: hom_1 = V.hom(W) + sage: hom_2 = W.hom(V) + sage: hom_1.is_in_h0(isom) + True + sage: hom_2.is_in_h0(isom^-1) True """ - mat_basis = self.hom(other).h0() + r = self.rank() + if other.rank() != r: + return None + V = self.direct_sum(other) + End = V.end() + A, to_A, from_A = End.global_algebra() + (S, to_S, from_S), factors = End._global_algebra_split() + if len(factors) > 1: + return None + split = factors[0] + if split.M.nrows() != 2: + return None + K = self._function_field + id_self = split.to_M(to_S(to_A( + diagonal_matrix(K, [1]*r + [0]*r)))) + im_self = id_self.transpose().echelon_form().rows() + im_self = [r for r in im_self if not r.is_zero()] + id_other = split.to_M(to_S(to_A( + diagonal_matrix(K, [0]*r + [1]*r)))) + im_other = id_other.transpose().echelon_form().rows() + im_other = [r for r in im_other if not r.is_zero()] + P = matrix(im_self + im_other).transpose() + F = split._K + s = split._n // 2 + id_s = identity_matrix(F, s) + zero_mat = zero_matrix(F, s) + isom = P * block_matrix([[zero_mat, zero_mat], [id_s, zero_mat]]) * P**-1 + isom = from_A(from_S(split.from_M(isom))) + return isom[r:,:r] + + def _indecomposable_power_split(self): + r""" + Check if self is a direct sum of copies of an indecomposable bundle. + If so, return this indecomposable `L` and an isomorphism from + ``L.direct_sum_repeat(s)`` to ``self``. + + TESTS :: + + sage: #long time (25 seconds) + sage: from vector_bundle import atiyah_bundle + sage: F. = FunctionField(GF(7)) + sage: R. = F[] + sage: K. = F.extension(y^2 - x^3 - x) + sage: V = atiyah_bundle(K, 2, 0) + sage: W = V.direct_sum_repeat(2) + sage: T = matrix(K, 4, 4, + ....: [x+1, 4, 3*x+4, 6*x+6, + ....: 4*x+5, 6*x+5, 2*x+1, 4*x+3, + ....: 3*x+1, 5*x, 3*x+1, x+6, + ....: 5*x+4, 6*x, 5*x+2, 6*x+6]) + sage: W = W.apply_isomorphism(T) + sage: ind, s, isom = W._indecomposable_power_split() + sage: s + 2 + sage: V._isomorphism_indecomposable(ind) is not None + True + sage: ind.direct_sum_repeat(s).hom(W).is_isomorphism(isom) + True + """ + End = self.end() + A, to_A, from_A = End.global_algebra() + (S, to_S, from_S), factors = End._global_algebra_split() + if len(factors) > 1: + return None, None + split = factors[0] + K = split._K + s = split._n + idems = [diagonal_matrix(K, [0]*i + [1] + [0]*(s-i-1)) + for i in range(s)] + idems = [from_A(from_S(split.from_M(idem))) for idem in idems] + images = [End.image(idem) for idem in idems] + factor = images[0][0] + isoms = [factor._isomorphism_indecomposable(im[0]) for im in images[1:]] + phis = ([images[0][1]] + + [im[1]*isom for im, isom in zip(images[1:], isoms)]) + return (factor, + len(phis), + block_matrix(self._function_field, [phis])) + + def split(self): + r""" + Return a list of indecomposable bundles ``inds``, a list of integers + `ns` and an isomorphism from + `\bigoplus_{\mathrm{ind} \in \mathrm{inds}} + \mathrm{ind}^{n_\mathrm{ind}}` to ``self``. + + EXAMPLES:: + + sage: # long time (20 seconds) + sage: from vector_bundle import atiyah_bundle, trivial_bundle + sage: F. = FunctionField(GF(7)) + sage: R. = F[] + sage: K. = F.extension(y^2 - x^3 - x) + sage: triv = trivial_bundle(K) + sage: V = atiyah_bundle(K, 2, 0) + sage: W = triv.direct_sum_repeat(2).direct_sum(V) + sage: T = matrix(K, 4, 4, + ....: [x+1, 4, 3*x+4, 6*x+6, + ....: 4*x+5, 6*x+5, 2*x+1, 4*x+3, + ....: 3*x+1, 5*x, 3*x+1, x+6, + ....: 5*x+4, 6*x, 5*x+2, 6*x+6]) + sage: W = W.apply_isomorphism(T) + sage: inds, ns, isom = W.split() + sage: b1 = [ind.rank() for ind in inds] == [1, 2] + sage: b2 = [ind.rank() for ind in inds] == [2, 1] + sage: b1 or b2 + True + sage: b1 = ns == [1, 2] + sage: b2 = ns == [2, 1] + sage: b1 or b2 + True + sage: sum = inds[0].direct_sum_repeat(ns[0]) + sage: sum = sum.direct_sum(inds[1].direct_sum_repeat(ns[1])) + sage: sum.hom(W).is_isomorphism(isom) + True + """ + K = self._function_field + End = self.end() + A, to_A, from_A = End.global_algebra() + (S, to_S, from_S), splits = End._global_algebra_split() + injections = [[from_A(from_S(s.from_M(diagonal_matrix( + s._K, + [0]*i + [1] + [0]*(s._n-1-i))))) + for i in range(s._n)] + for s in splits] + images = [[End.image(inj) for inj in injs] for injs in injections] + inds = [im[0][0] for im in images] + phis = [[ims[0][1]] + + [im[1] * ind._isomorphism_indecomposable(im[0]) + for im in ims[1:]] + for ind,ims in zip(inds, images)] + ns = [len(injs) for injs in injections] + isom = block_matrix([sum(phis, [])]) + return inds, ns, isom + + def _isomorphism_to_small_field(self, other): + r""" + Return an isomorphism from self to other if it exists and None otherwise + """ + inds_self, ns_self, isom_self = self.split() + inds_other, ns_other, isom_other = other.split() + if sorted(ns_self) != sorted(ns_other): + return None + isoms = [[ind_self._isomorphism_indecomposable(ind_other) + for ind_other in inds_other] for ind_self in inds_self] + fits = [[i for i, iso in enumerate(isos) if iso is not None] + for isos in isoms] + if any([len(fit) != 1 for fit in fits]): + return None + #We now know that self and other are isomorphic + fits = [fit[0] for fit in fits] + ranks_self = [ind.rank() for ind in inds_self] + ranks_other = [ind.rank() for ind in inds_other] + s = len(inds_self) + blocks = [block_diagonal_matrix([isoms[i][fits[i]]]*ns_self[i]) + for i in range(s)] + isom = block_matrix([[blocks[i] if j == fits[i] else 0 + for j in range(s)] + for i in range(s)]) + return isom_other * isom * isom_self**-1 + + def isomorphism_to(self, other): + r""" + Return an isomorphism from self to other if it exists and None otherwise + + EXAMPLES :: + + sage: from vector_bundle import (trivial_bundle, canonical_bundle, + ....: atiyah_bundle) + sage: F. = FunctionField(GF(3)) + sage: R. = F[] + sage: K. = F.extension(y^2 - x^3 - x) + sage: triv = trivial_bundle(K) + sage: can = canonical_bundle(K) + sage: iso = triv.isomorphism_to(can) + sage: triv.hom(can).is_isomorphism(iso) + True + sage: V = can.direct_sum(atiyah_bundle(K, 2, 0, can))\ + ....: .direct_sum(triv) + sage: W = can.direct_sum(atiyah_bundle(K, 2, 0)).direct_sum(can) + sage: iso = V.isomorphism_to(W) + sage: V.hom(W).is_isomorphism(iso) + True + + WARNING: + + Not well implemented for infinite fields: need to specify how to chose + random elements and adequatly set the sample size. + """ + s = len(self.end().h0()) k = self._function_field.constant_base_field() - return any([sum([(c*mat).is_unit() for c, mat in zip(vec, mat_basis)]) - for vec in ProjectiveSpace(len(mat_basis)-1, k)]) + if k.cardinality() > s: + return self._isomorphism_to_large_field( + other, + integer_ceil(60/log(k.cardinality()/s))) + return self._isomorphism_to_small_field(other) + + def apply_isomorphism(self, isom): + r""" + Isom is an invertible square matrix of order ``self.rank()``. + Return the image of ``self`` by ``isom``. + + EXAMPLES :: + + sage: from vector_bundle import (trivial_bundle, canonical_bundle, + ....: atiyah_bundle) + sage: F. = FunctionField(GF(3)) + sage: R. = F[] + sage: K. = F.extension(y^2 - x^3 - x) + sage: triv = trivial_bundle(K) + sage: can = canonical_bundle(K) + sage: iso = triv.isomorphism_to(can) + sage: triv.hom(can).is_isomorphism(iso) + True + sage: V = can.direct_sum(atiyah_bundle(K, 2, 0, can))\ + ....: .direct_sum(triv) + sage: W = can.direct_sum(atiyah_bundle(K, 2, 0)).direct_sum(can) + sage: iso = V.isomorphism_to(W) + sage: V.apply_isomorphism(iso) == W + True + """ + return VectorBundle(self._function_field, + self._ideals, + isom * self._g_finite, + isom * self._g_infinite, + check=False)