DeRhamComputation/as_covers/as_polyforms.sage

class as_polyform:
    def __init__(self, form, mult):
        '''Elements of H^0(Ω^⊗n). Usage: mult is n, form should be as_function.'''
        self.form = form
        self.curve = form.curve
        self.mult =  mult
        
    def __add__(self, other):
        return as_polyform(self.form + other.form, self.mult)
    
    def __repr__(self):
        return '(' + str(self.form) + ') dx⊗' + str(self.mult)
    
    def expansion_at_infty(self, place):
        return self.form.expansion_at_infty(place=place)*(self.curve.dx.expansion_at_infty(place=place))^(self.mult)
    
    def reduce(self):
        AS = self.curve
        reduced = self.form.reduce()
        return as_polyform(reduced, self.mult)
    
    def coordinates(self, basis = 0):
        """Find coordinates of the given holomorphic form self in terms of the basis forms in a list holo."""
        AS = self.curve
        if basis == 0:
            basis = AS.holo_polydifferentials_basis(mult = self.mult)
        RxyzQ, Rxyz, x, y, z = AS.fct_field
        # We need to have only polynomials to use monomial_coefficients in linear_representation_polynomials,
        # and sometimes basis elements have denominators. Thus we multiply by them.
        denom = LCM([denominator(omega.form.function) for omega in basis])
        basis = [denom*omega.form.function for omega in basis]
        self_reduced = self.reduce()
        self_with_no_denominator = denom*self_reduced.form.function
        return linear_representation_polynomials(Rxyz(self_with_no_denominator), [Rxyz(omega) for omega in basis])
    
    def group_action(self, elt):
        return as_polyform(self.form.group_action(elt), self.mult)
    
    
def as_holo_polydifferentials_basis(AS, mult, threshold = 8):
    '''Give the basis of H^0(Ω^⊗n) for n = mult.'''
    v = AS.dx.valuation()
    result = AS.at_most_poles(mult*v, threshold=threshold)
    result = [as_polyform(omega, mult) for omega in result]
    if mult == 1 and len(result) < AS.genus():
        raise ValueError('Increase threshold, not all forms found.')
    if mult > 1 and len(result) < (2*mult - 1)*(AS.genus() - 1):
        raise ValueError('Increase threshold, not all forms found.')
    return result

as_cover.holo_polydifferentials_basis = as_holo_polydifferentials_basis

def as_symmetric_power_basis(AS, n, threshold = 8):
    '''Give the basis of H^0(Ω)^⊙n for n = mult.'''
    g = AS.genus()
    B0 = AS.holomorphic_differentials_basis(threshold=threshold)
    from itertools import product
    indices = [list(range(g)) for i in range(n)]
    indices_nonrepeating = []
    for i in product(*indices):
        if non_decreasing(i):
            indices_nonrepeating += [i]
    result = []
    for i in indices_nonrepeating:
        tensor_form = [B0[i[j]] for j in range(n)]
        tensor_form = [1] + tensor_form
        result += [tensor_form]
    result = [as_symmetric_product_forms([a]) for a in result]
    return result

def as_canonical_ideal(AS, n, threshold=8):
    '''Return the n-th homogeneous part of the canonical ideal.'''
    B0 = AS.holomorphic_differentials_basis(threshold=threshold)
    F = AS.base_ring
    g = AS.genus()
    RxyzQ, Rxyz, x, y, z = AS.fct_field
    from itertools import product
    B1 = as_symmetric_power_basis(AS, n, threshold = threshold)
    B2 = AS.holo_polydifferentials_basis(n, threshold = threshold)
    g = AS.genus()
    r = len(B2)
    M = matrix(F, len(B1), r)
    for i in range(0, len(B1)):
        c = B1[i].multiply()
        #return M, c.coordinates(basis=B2)
        M[i, :] = vector(c.coordinates(basis=B2))
    K = M.kernel().basis()
    result = []
    for v in K:
        kernel_vector = 0*B1[0]
        for i in range(len(v)):
            kernel_vector += v[i]*B1[i]
        result += [kernel_vector]
    return result

as_cover.canonical_ideal = as_canonical_ideal

def as_canonical_ideal_polynomials(AS, n, threshold=8):
    '''Return the polynomials defining n-th homogeneous part of the canonical ideal.'''
    return [a.polynomial() for a in AS.canonical_ideal(n, threshold=threshold)]

as_cover.canonical_ideal_polynomials = as_canonical_ideal_polynomials

class as_symmetric_product_forms:
    def __init__(self, forms_and_coeffs):
        '''Elements of forms_and_coeffs are of the form [coeff, form_1, ..., form_n]'''
        self.n = len(forms_and_coeffs[0]) - 1
        forms_and_coeffs1 = []
        for atuple in forms_and_coeffs:
            if atuple[0] != 0 and atuple[1:] not in [a[1:] for a in forms_and_coeffs1]:
                forms_and_coeffs1 += [atuple]
            elif atuple[1:] in [a[1:] for a in forms_and_coeffs1]:
                i = [a[1:] for a in forms_and_coeffs1].index(atuple[1:])
                forms_and_coeffs1[i][0] += atuple[0]
        if len(forms_and_coeffs1) == 0:
            forms_and_coeffs1 = [[0] + forms_and_coeffs[0][1:]]
        self.tuples = forms_and_coeffs1
        self.curve = forms_and_coeffs1[0][1].curve
        
    def coordinates(self, basis = 0):
        AS = self.curve
        g = AS.genus()
        n = self.n
        F = AS.base_ring
        if basis == 0:
            basis = AS.holomorphic_differentials_basis()
        from itertools import product
        indices = [list(range(g)) for i in range(self.n)]
        result = {i : 0 for i in product(*indices)}
        for atuple in self.tuples:
            coors = [omega.coordinates(basis = basis) for omega in atuple[1:]]
            for i in product(*indices):
                aux_product = 1
                for j in range(n):
                    aux_product *= coors[j][i[j]]
                result[i] += atuple[0]*aux_product
        return result
    
    def __repr__(self):
        result = ''
        for atuple in self.tuples:
            if atuple[0] !=0:
                result += str(atuple[0]) + ' * '
                for j in range(1, self.n+1):
                    result += "(" + str(atuple[j]) + ")"
                    if j != self.n:
                        result + '⊗'
            result += ' + '
        return result
    
    def __add__(self, other):
        return as_symmetric_product_forms(self.tuples + other.tuples)
    
    def __rmul__(self, other):
        AS = self.curve
        F = AS.base_ring
        if isinstance(other, int) or other in ZZ or other in F:
            aux_tuples = []
            for atuple in self.tuples:
                atuple1 = [other * atuple[0]] + atuple[1:]
                aux_tuples += [atuple1]
            return as_symmetric_product_forms(aux_tuples)
    
    def multiply(self):
        '''The map H^0(Ω)^⊙n ---> H^0(Ω^⊗n)'''
        n = self.n
        AS = self.curve
        RxyzQ, Rxyz, x, y, z = AS.fct_field
        result = as_polyform(0*AS.x, n)
        for atuple in self.tuples:
            aux_product = Rxyz(1)
            for fct in atuple[1:]:
                aux_product = aux_product * fct.form
            aux_product = as_function(AS, aux_product)
            aux_product = as_reduction(AS, aux_product)
            aux_product = as_function(AS, aux_product)
            result += as_polyform(atuple[0]*aux_product, n)
        return result

    def polynomial(self):
        '''Return the associated polynomial.'''
        AS = self.curve
        F = AS.base_ring
        g = AS.genus()
        M = self.coordinates()
        n = self.n
        Rg = PolynomialRing(F, 'X', g)
        X = Rg.gens()
        from itertools import product
        indices = [list(range(g)) for i in range(self.n)]
        result = Rg(0)
        for i in product(*indices):
            aux_product = Rg(1)
            for j in range(n):
                aux_product *= X[i[j]]
            result += M[i] * aux_product
        return result
    
    def group_action(self, elt):
        p = self.curve.base_ring.characteristic()
        n = self.curve.height
        aux_tuples = []
        for atuple in self.tuples:
            aux_tuple = [atuple[0]] + [a.group_action(elt) for a in atuple[1:]]
            aux_tuples += [aux_tuple]
        return as_symmetric_product_forms(aux_tuples)
    
def non_decreasing(L):
    return all(x<=y for x, y in zip(L, L[1:]))

def as_matrices_group_action_canonical_ideal(AS, mult, threshold = 8):
    '''Return the group action matrices for the n-th homogeneous part of the canonical ideal.'''
    K = as_canonical_ideal(AS, mult, threshold = threshold)
    n = AS.height
    F = AS.base_ring
    K_polynomials = [a.polynomial() for a in K]
    r = len(K)
    matrices = []
    for i in range(n):
        M = matrix(F, r, r)
        K_group_action_polynomials = [a.group_action([j == i for j in range(n)]).polynomial() for a in K]
        for i in range(r):
            M[i, :] = vector(linear_representation_polynomials(K_group_action_polynomials[i], K_polynomials))
        matrices += [M]
    return matrices

as_cover.group_action_canonical_ideal = as_matrices_group_action_canonical_ideal