DeRhamComputation/sage/superelliptic/superelliptic_cech_class.sage

251 lines
8.0 KiB
Python
Raw Normal View History

2022-11-18 15:00:34 +01:00
class superelliptic_cech:
def __init__(self, C, omega, fct):
self.omega0 = omega
self.omega8 = omega - fct.diffn()
self.f = fct
self.curve = C
def __add__(self, other):
C = self.curve
return superelliptic_cech(C, self.omega0 + other.omega0, self.f + other.f)
def __sub__(self, other):
C = self.curve
return superelliptic_cech(C, self.omega0 - other.omega0, self.f - other.f)
def __rmul__(self, constant):
C = self.curve
w1 = self.omega0.form
f1 = self.f.function
w2 = superelliptic_form(C, constant*w1)
f2 = superelliptic_function(C, constant*f1)
return superelliptic_cech(C, w2, f2)
def __repr__(self):
return "(" + str(self.omega0) + ", " + str(self.f) + ", " + str(self.omega8) + ")"
def verschiebung(self):
C = self.curve
omega = self.omega0
F = C.base_ring
Rx.<x> = PolynomialRing(F)
return superelliptic_cech(C, omega.cartier(), superelliptic_function(C, Rx(0)))
def frobenius(self):
C = self.curve
fct = self.f.function
p = C.characteristic
Rx.<x> = PolynomialRing(F)
return superelliptic_cech(C, superelliptic_form(C, Rx(0)), superelliptic_function(C, fct^p))
def coordinates(self):
C = self.curve
F = C.base_ring
m = C.exponent
Rx.<x> = PolynomialRing(F)
Fx = FractionField(Rx)
FxRy.<y> = PolynomialRing(Fx)
g = C.genus()
degrees_holo = C.degrees_holomorphic_differentials()
degrees_holo_inv = {b:a for a, b in degrees_holo.items()}
degrees0 = C.degrees_de_rham0()
degrees0_inv = {b:a for a, b in degrees0.items()}
degrees1 = C.degrees_de_rham1()
degrees1_inv = {b:a for a, b in degrees1.items()}
basis = C.de_rham_basis()
omega = self.omega0
fct = self.f
if fct.function == Rx(0) and omega.form != Rx(0):
for j in range(1, m):
omega_j = Fx(omega.jth_component(j))
if omega_j != Fx(0):
d = degree_of_rational_fctn(omega_j, F)
index = degrees_holo_inv[(d, j)]
a = coeff_of_rational_fctn(omega_j, F)
a1 = coeff_of_rational_fctn(basis[index].omega0.jth_component(j), F)
elt = self - (a/a1)*basis[index]
return elt.coordinates() + a/a1*vector([F(i == index) for i in range(0, 2*g)])
for j in range(1, m):
fct_j = Fx(fct.jth_component(j))
if (fct_j != Rx(0)):
d = degree_of_rational_fctn(fct_j, p)
if (d, j) in degrees1.values():
index = degrees1_inv[(d, j)]
a = coeff_of_rational_fctn(fct_j, F)
elt = self - (a/m)*basis[index]
return elt.coordinates() + a/m*vector([F(i == index) for i in range(0, 2*g)])
if d<0:
a = coeff_of_rational_fctn(fct_j, F)
h = superelliptic_function(C, FxRy(a*y^j*x^d))
elt = superelliptic_cech(C, self.omega0, self.f - h)
return elt.coordinates()
if (fct_j != Rx(0)):
G = superelliptic_function(C, y^j*x^d)
a = coeff_of_rational_fctn(fct_j, F)
elt =self - a*superelliptic_cech(C, diffn(G), G)
return elt.coordinates()
return vector(2*g*[0])
def is_cocycle(self):
w0 = self.omega0
w8 = self.omega8
fct = self.f
if not w0.is_regular_on_U0() and not w8.is_regular_on_Uinfty():
return('w0 & w8')
if not w0.is_regular_on_U0():
return('w0')
if not w8.is_regular_on_Uinfty():
return('w8')
if w0.is_regular_on_U0() and w8.is_regular_on_Uinfty():
return 1
return 0
#Auxilliary. If f = f1/f2 is a rational function, return deg f_1 - deg f_2.
def degree_of_rational_fctn(f, F):
Rx.<x> = PolynomialRing(F)
Fx = FractionField(Rx)
f = Fx(f)
f1 = f.numerator()
f2 = f.denominator()
d1 = f1.degree()
d2 = f2.degree()
return(d1 - d2)
#Auxilliary. If f = f1/f2 is a rational function, return (leading coeff of f1)/(leading coeff of f2).
def coeff_of_rational_fctn(f, F):
Rx.<x> = PolynomialRing(F)
Fx = FractionField(Rx)
f = Fx(f)
if f == Rx(0):
return 0
f1 = f.numerator()
f2 = f.denominator()
d1 = f1.degree()
d2 = f2.degree()
a1 = f1.coefficients(sparse = false)[d1]
a2 = f2.coefficients(sparse = false)[d2]
return(a1/a2)
#Auxilliary. Given polynomial f(x) and integer d, return
#coefficient of x^d in f (and 0 is deg(f)<d).
def coff(f, d):
lista = f.coefficients(sparse = false)
if len(lista) <= d:
return 0
return lista[d]
#Auxilliary. Given polynomial f(x) = \sum_i a_i x^i and integer i, return
#only \sum_{j >= i+1} a_j x^{j - i -1}
def cut(f, i):
R = f.parent()
coeff = f.coefficients(sparse = false)
return sum(R(x^(j-i-1)) * coeff[j] for j in range(i+1, f.degree() + 1))
def polynomial_part(p, h):
F = GF(p)
Rx.<x> = PolynomialRing(F)
h = Rx(h)
result = Rx(0)
for i in range(0, h.degree()+1):
if (i%p) == p-1:
power = Integer((i-(p-1))/p)
result += Integer(h[i]) * x^(power)
return result
#Find delta-th root of unity in field F
def root_of_unity(F, delta):
Rx.<x> = PolynomialRing(F)
cyclotomic = x^(delta) - 1
for root, a in cyclotomic.roots():
powers = [root^d for d in delta.divisors() if d!= delta]
if 1 not in powers:
return root
def preimage(U, V, M): #preimage of subspace U under M
basis_preimage = M.right_kernel().basis()
imageU = U.intersection(M.transpose().image())
basis = imageU.basis()
for v in basis:
w = M.solve_right(v)
basis_preimage = basis_preimage + [w]
return V.subspace(basis_preimage)
def image(U, V, M):
basis = U.basis()
basis_image = []
for v in basis:
basis_image += [M*v]
return V.subspace(basis_image)
def flag(F, V, p, test = 0):
dim = F.dimensions()[0]
space = VectorSpace(GF(p), dim)
flag_subspaces = (dim+1)*[0]
flag_used = (dim+1)*[0]
final_type = (dim+1)*['?']
flag_subspaces[dim] = space
flag_used[dim] = 1
while 1 in flag_used:
index = flag_used.index(1)
flag_used[index] = 0
U = flag_subspaces[index]
U_im = image(U, space, V)
d_im = U_im.dimension()
final_type[index] = d_im
U_pre = preimage(U, space, F)
d_pre = U_pre.dimension()
if flag_subspaces[d_im] == 0:
flag_subspaces[d_im] = U_im
flag_used[d_im] = 1
if flag_subspaces[d_pre] == 0:
flag_subspaces[d_pre] = U_pre
flag_used[d_pre] = 1
if test == 1:
print('(', final_type, ')')
for i in range(0, dim+1):
if final_type[i] == '?' and final_type[dim - i] != '?':
i1 = dim - i
final_type[i] = final_type[i1] - i1 + dim/2
final_type[0] = 0
for i in range(1, dim+1):
if final_type[i] == '?':
prev = final_type[i-1]
if prev != '?' and prev in final_type[i+1:]:
final_type[i] = prev
for i in range(1, dim+1):
if final_type[i] == '?':
final_type[i] = min(final_type[i-1] + 1, dim/2)
if is_final(final_type, dim/2):
return final_type[1:dim/2 + 1]
print('error:', final_type[1:dim/2 + 1])
def is_final(final_type, dim):
n = len(final_type)
if final_type[0] != 0:
return 0
if final_type[n-1] != dim:
return 0
for i in range(1, n):
if final_type[i] != final_type[i - 1] and final_type[i] != final_type[i - 1] + 1:
return 0
return 1