vector-bundles-sagemath/vector_bundle/algebras.py

560 lines
20 KiB
Python

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