class as_cover:
    def __init__(self, C, list_of_fcts, prec = 10):
        self.quotient = C
        self.functions = list_of_fcts
        self.height = len(list_of_fcts)
        F = C.base_ring
        self.base_ring = F
        p = C.characteristic
        self.characteristic = p
        self.prec = prec
        f = C.polynomial
        m = C.exponent
        r = f.degree()
        delta = GCD(m, r)
        self.nb_of_pts_at_infty = delta
        Rxy.<x, y> = PolynomialRing(F, 2)
        Rt.<t> = LaurentSeriesRing(F, default_prec=prec)

        all_x_series = []
        all_y_series = []
        all_z_series = []
        all_dx_series = []
        all_jumps = []

        for i in range(delta):
            x_series = superelliptic_function(C, x).expansion_at_infty(i = i, prec=prec)
            y_series = superelliptic_function(C, y).expansion_at_infty(i = i, prec=prec)
            z_series = []
            jumps = []
            n = len(list_of_fcts)
            list_of_power_series = [g.expansion_at_infty(i = i, prec=prec) for g in list_of_fcts]
            for i in range(n):
                power_series = list_of_power_series[i]
                jump, correction, t_old, z = artin_schreier_transform(power_series, prec = prec)
                x_series = x_series(t = t_old)
                y_series = y_series(t = t_old)
                z_series = [zi(t = t_old) for zi in z_series]
                z_series += [z]
                jumps += [jump]
                list_of_power_series = [g(t = t_old) for g in list_of_power_series]
            
            all_jumps += [jumps]
            all_x_series += [x_series]
            all_y_series += [y_series]
            all_z_series += [z_series]
            all_dx_series += [x_series.derivative()]
        self.jumps = all_jumps
        self.x = all_x_series
        self.y = all_y_series
        self.z = all_z_series
        self.dx = all_dx_series
        
    def __repr__(self):
        n = self.height
        p = self.characteristic
        if n==1:
            return "(Z/p)-cover of " + str(self.quotient)+" with the equation:\n z^" + str(p) + " - z = " + str(self.functions[0])
        
        result = "(Z/p)^"+str(self.height)+ "-cover of " + str(self.quotient)+" with the equations:\n"
        for i in range(n):
            result += 'z' + str(i) + "^" + str(p) + " - z" + str(i) + " = " + str(self.functions[i]) + "\n"
        return result
        
    def genus(self):
        jumps = self.jumps
        gY = self.quotient.genus()
        n = self.height
        delta = self.nb_of_pts_at_infty
        p = self.characteristic
        return p^n*gY + (p^n - 1)*(delta - 1) + sum(p^(n-j-1)*(jumps[i][j]-1)*(p-1)/2 for j in range(n) for i in range(delta))
    
    def exponent_of_different(self, i = 0):
        jumps = self.jumps
        n = self.height
        delta = self.nb_of_pts_at_infty
        p = self.characteristic
        return sum(p^(n-j-1)*(jumps[i][j]+1)*(p-1) for j in range(n))

    def exponent_of_different_prim(self, i = 0):
        jumps = self.jumps
        n = self.height
        delta = self.nb_of_pts_at_infty
        p = self.characteristic
        return sum(p^(n-j-1)*(jumps[i][j])*(p-1) for j in range(n))
    
    def holomorphic_differentials_basis(self, threshold = 8):
        from itertools import product
        x_series = self.x
        y_series = self.y
        z_series = self.z
        dx_series = self.dx
        delta = self.nb_of_pts_at_infty
        p = self.characteristic
        n = self.height
        prec = self.prec
        C = self.quotient
        F = self.base_ring
        m = C.exponent
        r = C.polynomial.degree()
        variable_names = 'x, y'
        for i in range(n):
            variable_names += ', z' + str(i)
        Rxyz = PolynomialRing(F, n+2, variable_names)
        x, y = Rxyz.gens()[:2]
        z = Rxyz.gens()[2:]
        RxyzQ = FractionField(Rxyz)
        Rt.<t> = LaurentSeriesRing(F, default_prec=prec)
        #Tworzymy zbiór S form z^i x^j y^k dx/y o waluacji >= waluacja z^(p-1)*dx/y
        S = []
        RQxyz = FractionField(Rxyz)
        pr = [list(GF(p)) for _ in range(n)]
        for i in range(0, threshold*r):
            for j in range(0, m):
                for k in product(*pr):
                    eta = as_form(self, x^i * prod(z[i1]^(k[i1]) for i1 in range(n))/y^j)
                    eta_exp = eta.expansion_at_infty()
                    S += [(eta, eta_exp)]

        forms = holomorphic_combinations(S)
        
        for i in range(1, delta):
            forms = [(omega, omega.expansion_at_infty(i = i)) for omega in forms]
            forms = holomorphic_combinations(forms)
        
        if len(forms) < self.genus():
            print("I haven't found all forms.")
            return holomorphic_differentials_basis(self, threshold = threshold + 1)
        if len(forms) > self.genus():
            print("Increase precision.")
        return forms

    def at_most_poles(self, pole_order, threshold = 8):
        """ Find fcts with pole order in infty's at most pole_order. Threshold gives a bound on powers of x in the function.
           If you suspect that you haven't found all the functions, you may increase it."""
        from itertools import product
        x_series = self.x
        y_series = self.y
        z_series = self.z
        delta = self.nb_of_pts_at_infty
        p = self.characteristic
        n = self.height
        prec = self.prec
        C = self.quotient
        F = self.base_ring
        m = C.exponent
        r = C.polynomial.degree()
        variable_names = 'x, y'
        for i in range(n):
            variable_names += ', z' + str(i)
        Rxyz = PolynomialRing(F, n+2, variable_names)
        x, y = Rxyz.gens()[:2]
        z = Rxyz.gens()[2:]
        RxyzQ = FractionField(Rxyz)
        Rt.<t> = LaurentSeriesRing(F, default_prec=prec)
        #Tworzymy zbiór S form z^i x^j y^k dx/y o waluacji >= waluacja z^(p-1)*dx/y
        S = []
        RQxyz = FractionField(Rxyz)
        pr = [list(GF(p)) for _ in range(n)]
        for i in range(0, threshold*r):
            for j in range(0, m):
                for k in product(*pr):
                    eta = as_function(self, x^i * prod(z[i1]^(k[i1]) for i1 in range(n))*y^j)
                    eta_exp = eta.expansion_at_infty()
                    S += [(eta, eta_exp)]

        forms = holomorphic_combinations_fcts(S, pole_order)

        for i in range(1, delta):
            forms = [(omega, omega.expansion_at_infty(i = i)) for omega in forms]
            forms = holomorphic_combinations_fcts(forms, pole_order)
                  
        return forms
    
    def magical_element(self, threshold = 8):
        list_of_elts = self.at_most_poles(self.exponent_of_different_prim(), threshold)
        result = []
        for a in list_of_elts:
            if a.trace().function != 0:
                result += [a]
        return result

    def pseudo_magical_element(self, threshold = 8):
        list_of_elts = self.at_most_poles(self.exponent_of_different(), threshold)
        result = []
        for a in list_of_elts:
            if a.trace().function != 0:
                result += [a]
        return result

    def at_most_poles_forms(self, pole_order, threshold = 8):
        """Find forms with pole order in all the points at infty equat at most to pole_order. Threshold gives a bound on powers of x in the form.
           If you suspect that you haven't found all the functions, you may increase it."""
        from itertools import product
        x_series = self.x
        y_series = self.y
        z_series = self.z
        delta = self.nb_of_pts_at_infty
        p = self.characteristic
        n = self.height
        prec = self.prec
        C = self.quotient
        F = self.base_ring
        m = C.exponent
        r = C.polynomial.degree()
        variable_names = 'x, y'
        for i in range(n):
            variable_names += ', z' + str(i)
        Rxyz = PolynomialRing(F, n+2, variable_names)
        x, y = Rxyz.gens()[:2]
        z = Rxyz.gens()[2:]
        RxyzQ = FractionField(Rxyz)
        Rt.<t> = LaurentSeriesRing(F, default_prec=prec)
        #Tworzymy zbiór S form z^i x^j y^k dx/y o waluacji >= waluacja z^(p-1)*dx/y
        S = []
        RQxyz = FractionField(Rxyz)
        pr = [list(GF(p)) for _ in range(n)]
        for i in range(0, threshold*r):
            for j in range(0, m):
                for k in product(*pr):
                    eta = as_form(self, x^i * prod(z[i1]^(k[i1]) for i1 in range(n))/y^j)
                    eta_exp = eta.expansion_at_infty()
                    S += [(eta, eta_exp)]

        forms = holomorphic_combinations_forms(S, pole_order)

        for i in range(1, delta):
            forms = [(omega, omega.expansion_at_infty(i = i)) for omega in forms]
            forms = holomorphic_combinations_forms(forms, pole_order)
                  
        return forms

def holomorphic_combinations(S):
    """Given a list S of pairs (form, corresponding Laurent series at some pt), find their combinations holomorphic at that pt."""
    C_AS = S[0][0].curve
    p = C_AS.characteristic
    F = C_AS.base_ring
    prec = C_AS.prec
    Rt.<t> = LaurentSeriesRing(F, default_prec=prec)
    RtQ = FractionField(Rt)
    minimal_valuation = min([g[1].valuation() for g in S])
    if minimal_valuation >= 0:
        return [s[0] for s in S]
    list_of_lists = [] #to będzie lista złożona z list współczynników część nieholomorficznych rozwinięcia form z S
    for eta, eta_exp in S:
        a = -minimal_valuation + eta_exp.valuation()
        list_coeffs = a*[0] + eta_exp.list() + (-minimal_valuation)*[0]
        list_coeffs = list_coeffs[:-minimal_valuation]
        list_of_lists += [list_coeffs]
    M = matrix(F, list_of_lists)
    V = M.kernel() #chcemy wyzerować części nieholomorficzne, biorąc kombinacje form z S


    # Sprawdzamy, jakim formom odpowiadają elementy V.
    forms = []
    for vec in V.basis():
        forma_holo = as_form(C_AS, 0)
        forma_holo_power_series = Rt(0)
        for vec_wspolrzedna, elt_S in zip(vec, S):
            eta = elt_S[0]
            #eta_exp = elt_S[1]
            forma_holo += vec_wspolrzedna*eta
            #forma_holo_power_series += vec_wspolrzedna*eta_exp
        forms += [forma_holo]
    return forms