class as_cover:
    def __init__(self, C, cover_template, list_of_fcts, branch_points = [], prec = 10):
        self.quotient = C
        self.cover_template = cover_template
        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
        #group acting
        self.height = cover_template.height
        self.group = cover_template.group
        #########
        f = C.polynomial
        m = C.exponent
        r = f.degree()
        delta = GCD(m, r)
        self.nb_of_pts_at_infty = delta
        self.branch_points = list(range(delta)) + branch_points
        Rxy.<x, y> = PolynomialRing(F, 2)
        Rt.<t> = LaurentSeriesRing(F, default_prec=prec)
        Rzf, zgen, fgen = cover_template.fct_field
        all_x_series = {}
        all_y_series = {}
        all_z_series = {}
        all_dx_series = {}
        all_jumps = {}

        for pt in self.branch_points:
            x_series = superelliptic_function(C, x).expansion(pt=pt, prec=prec)
            y_series = superelliptic_function(C, y).expansion(pt=pt, prec=prec)
            z_series = []
            jumps = []
            n = len(list_of_fcts)
            list_of_power_series = [g.expansion(pt=pt, prec=prec) for g in list_of_fcts]
            for j in range(n):
                ####
                #
                power_series = Rzf(cover_template.fcts[j]).subs({zgen[i] : z_series[i] for i in range(j)} | {zgen[i] : 0 for i in range(j, n)} | {fgen[i] : list_of_power_series[i] for i in range(n)})
                ####
                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[pt] = jumps
            all_x_series[pt] = x_series
            all_y_series[pt] = y_series
            all_z_series[pt] = z_series
            all_dx_series[pt] = x_series.derivative()
        self.jumps = all_jumps
        self.x_series = all_x_series
        self.y_series = all_y_series
        self.z_series = all_z_series
        self.dx_series = all_dx_series
        ##############
        #Function field
        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)
        self.fct_field = (RxyzQ, Rxyz, x, y, z)
        self.x = as_function(self, x)
        self.y = as_function(self, y)
        self.z = [as_function(self, z[j]) for j in range(n)]
        self.dx = as_form(self, 1)
        self.one = as_function(self, 1)
        Rzf, zgen, fgen = cover_template.fct_field
        subs_fs = {zgen[i] : z[i]}| {fgen[i] : RxyzQ(list_of_fcts[i].function) for i in range(n)}
        self.rhs = [RxyzQ(cover_template.fcts[i].subs(subs_fs)) for i in range(n)]
        #####
        ##### We compute now the differentials dz[i]
        y_super = superelliptic_function(C, y)
        dy_super = y_super.diffn().form
        dz = []
        for i in range(n):
            aux_fct = self.rhs[i]
            result = 0
            for j in range(i):
                result += aux_fct.derivative(z[j])*dz[j]
            result += aux_fct.derivative(x)
            result += aux_fct.derivative(y)*dy_super
            dz += [-result]
        self.dz = dz

    def __repr__(self):
        n = self.height
        p = self.characteristic
        result = "("+self.group.short_name+")"
        result += "-cover of " + str(self.quotient)+" with the equations: \n"
        for i in range(n):
            result += 'z'+str(i)+'^p - z'+str(i) + ' = '
            aux = str(self.cover_template.fcts[i])
            for t in range(n):
                aux = aux.replace("f"+str(t), "(" + str(self.functions[t]) + ")")
            result += aux + '\n'
        return result
        
    def genus(self):
        jumps = self.jumps
        gY = self.quotient.genus()
        n = self.height
        p = self.characteristic
        return p^n*(gY-1) + 1 + 1/2*sum(self.exponent_of_different(place)*len(self.fiber(place)) for place in self.branch_points)
    
    def exponent_of_different(self, place = 0):
        jumps = self.jumps
        n = self.height
        p = self.characteristic
        dd = [0]
        for i in range(1, n+1):
            if jumps[place][i-1] == 0:
                dd += [dd[i-1]]
            else:
                dd += [(jumps[place][i-1]+1)*(p-1) + p*dd[i-1]]
        return dd[n]

    def exponent_of_different_prim(self, place = 0):
        jumps = self.jumps
        n = self.height
        p = self.characteristic
        dd = [0]
        for i in range(1, n+1):
            if jumps[place][i-1] == 0:
                dd += [dd[i-1]]
            else:
                dd += [jumps[place][i-1]*(p-1) + p*dd[i-1]]
        return dd[n]
    
    def exponent_of_different_bis(self, place = 0):
        jumps = self.jumps
        n = self.height
        p = self.characteristic
        dd = [0]
        for i in range(1, n+1):
            if jumps[place][i-1] == 0:
                dd += [dd[i-1]]
            else:
                dd += [(jumps[place][i-1]-1)*(p-1) + p*dd[i-1]]
        return dd[n]
    
    def holomorphic_differentials_basis(self, threshold = 8):
        from itertools import product
        x_series = self.x_series
        y_series = self.y_series
        z_series = self.z_series
        dx_series = self.dx_series
        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()
        RxyzQ, Rxyz, x, y, z = self.fct_field
        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 = []
        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(pt=self.branch_points[0])
                    S += [(eta, eta_exp)]

        forms = holomorphic_combinations(S)
        
        for pt in self.branch_points[1:]:
            forms = [(omega, omega.expansion(pt=pt)) for omega in forms]
            forms = holomorphic_combinations(forms)
        
        if len(forms) < self.genus():
            print("I haven't found all forms, only ", len(forms), " of ", self.genus())
            raise ValueError("Increase threshold.")
            #return holomorphic_differentials_basis(self, threshold = threshold + 1)
        if len(forms) > self.genus():
            raise ValueError("Increase precision.")
        return forms
    
    def cartier_matrix(self, prec=50):
        g = self.genus()
        F = self.base_ring
        M = matrix(F, g, g)
        for i, omega in enumerate(self.holomorphic_differentials_basis()):
            M[:, i] = vector(omega.cartier().coordinates())
        return M

    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_series
        y_series = self.y_series
        z_series = self.z_series
        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()
        RxyzQ, Rxyz, x, y, z = self.fct_field
        F = C.base_ring
        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(place = 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_series
        y_series = self.y_series
        z_series = self.z_series
        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()
        RxyzQ, Rxyz, x, y, z = self.fct_field
        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 pt in self.branch_points[1:]:
            forms = [(omega, omega.expansion(pt=pt)) for omega in forms]
            forms = holomorphic_combinations_forms(forms, pole_order)
                  
        return forms
    
    def uniformizer(self, place = 0, threshold = 10):
        '''Return uniformizer of curve self at place-th place at infinity.'''
        p = self.characteristic
        n = self.height
        F = self.base_ring
        list_of_fcts = self.at_most_poles(threshold)
        list_of_fcts2 = self.z + [self.x, self.y]
        for f1 in list_of_fcts:
            for f2 in list_of_fcts2:
                d, a, b = xgcd(f1.valuation(place), f2.valuation(place))
                if d == 1:
                    return f1^a*f2^b
        raise ValueError("Increase threshold.")
        
    def quasiuniformizer(self, place = 0, threshold = 10):
        '''Return an element with valuation coprime to p.'''
        p = self.characteristic
        n = self.height
        F = self.base_ring
        rr_space = self.at_most_poles(threshold)
        list_of_fcts = [ff for ff in rr_space if ff.valuation(place)%p != 0]
        list_of_fcts2 = [len(str(ff)) for ff in list_of_fcts]
        i_min = list_of_fcts2.index(min(list_of_fcts2))
        result = list_of_fcts.pop(i_min)
        flag = 1
        while flag == 1:
            flag = 0
            for g in self.group.elts:
                if result.group_action(g) == result and g != self.group.one:
                    flag = 1
            if flag == 1:
                list_of_fcts2 = [len(str(ff)) for ff in list_of_fcts]
                i_min = list_of_fcts2.index(min(list_of_fcts2))
                result = list_of_fcts.pop(i_min)
        return result

    def stabilizer(self, place = 0):
        result = []
        for g in self.group.elts:
            flag = 1
            for i in range(self.height):
                if self.z[i].valuation(place = place) > 0:
                    fct = self.z[i]
                elif self.z[i].valuation(place = place) < 0:
                    fct = self.one/self.z[i]
                if fct.group_action(g).valuation(place = place) <= 0:
                    flag = 0
            if flag:
                result += [g]
        return result
    
    def fiber(self, place = 0):
        'Gives representatives for the quotient G/G_P for given place. Those are in bijection with the fiber.'
        result = [AS.group.one]
        p = self.characteristic
        G = self.group
        H = self.stabilizer(place = place)
        for g in self.group.elts:
            if g != AS.group.one:
                flag = 1
                for v in result:
                    if (G.elt(g)*(-G.elt(v))).as_tuple in H:
                        flag = 0
                if flag:
                    result += [g]
        return result

    def ith_ramification_gp(self, i, place = 0, quasiuniformizer = 0, threshold = 20):
        '''Find ith ramification group at place at infty of nb place.'''
        G = self.group.elts
        p = self.characteristic
        Gi = [G[0]]
#        list_of_fcts = self.z #[self.one/zz for zz in AS.z if zz.valuation(place) < 0]
#        list_of_fcts += [zz^2 for zz in self.z]
        if isinstance(quasiuniformizer, int) or isinstance(quasiuniformizer, Integer):
            t = self.quasiuniformizer(place, threshold=threshold)
        else:
            t = quasiuniformizer
        for g in G:
            if g != G[0]:
                tg = t.group_action(g)
                v = (tg - t).valuation(place)
                if v >= i + t.valuation(place):
                    Gi += [g]
        return Gi
    
    def ramification_jumps(self, place = 0, quasiuniformizer = 0, threshold = 20):
        '''Return list of lower ramification jumps at at place at infty of nb place.'''
        G = self.stabilizer(place=place)
        p = self.characteristic
        Gi = [G[0]]
#        list_of_fcts = self.z #[self.one/zz for zz in AS.z if zz.valuation(place) < 0]
#        list_of_fcts += [zz^2 for zz in self.z]
        if isinstance(quasiuniformizer, int) or isinstance(quasiuniformizer, Integer):
            t = self.quasiuniformizer(place, threshold=threshold)
        else:
            t = quasiuniformizer
        result = []
        i = 0
        while len(G) > 1:
            Gi = [G[0]]
            for g in G:
                if g != G[0]:
                    tg = t.group_action(g)
                    v = (tg - t).valuation(place)
                    if v >= i + t.valuation(place):
                        Gi += [g]
            if len(Gi) < len(G):
                result += [i-1]
            G = Gi
            i+=1
        return result
    
    def upper_ramification_jumps(self, place = 0, quasiuniformizer = 0, threshold = 20):
        lj = self.ramification_jumps(place = place, quasiuniformizer = quasiuniformizer, threshold = threshold)
        result = []
        result += [lj[0]]
        for j in range(1, len(lj)):
            aux = len(self.stabilizer(place=place))//len(self.ith_ramification_gp(lj[j], place = place, quasiuniformizer = quasiuniformizer, threshold = threshold))
            result += [result[j-1] + (lj[j] - lj[j-1])//aux]
        return result
    
    def a_number(self):
        g = self.genus()
        return g - self.cartier_matrix().rank()

    def cohomology_of_structure_sheaf_basis(self, holo_basis = 0, threshold = 8):
        if holo_basis == 0:
            holo_basis = self.holomorphic_differentials_basis(threshold = threshold)
        from itertools import product
        x_series = self.x_series
        y_series = self.y_series
        z_series = self.z_series
        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()
        RxyzQ, Rxyz, x, y, z = self.fct_field
        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
        result_fcts = []
        V = VectorSpace(F,self.genus())
        S = V.subspace([])
        RQxyz = FractionField(Rxyz)
        pr = [list(GF(p)) for _ in range(n)]
        i = 0
        while len(result_fcts) < self.genus():
            for j in range(0, m):
                for k in product(*pr):
                    f = as_function(self, prod(z[i1]^(k[i1]) for i1 in range(n))/x^i*y^j)
                    f_products = [omega.serre_duality_pairing(f) for omega in holo_basis]
                    if vector(f_products) not in S:
                        S = S+V.subspace([V(f_products)])
                        result_fcts += [f]
            i += 1
        return result_fcts

    def lift_to_de_rham(self, fct, threshold = 30):
        '''Given function fct, find form eta regular on affine part such that eta - d(fct) is regular in infty. (Works for one place at infty now)'''
        from itertools import product
        x_series = self.x_series
        y_series = self.y_series
        z_series = self.z_series
        dx_series = self.dx_series
        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()
        RxyzQ, Rxyz, x, y, z = self.fct_field
        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 = [(fct.diffn(), fct.diffn().expansion_at_infty())]
        pr = [list(GF(p)) for _ in range(n)]
        holo = self.holomorphic_differentials_basis(threshold = threshold)
        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)
        if len(forms) <= self.genus():
            raise ValueError("Increase threshold!")
        for omega in forms:
            for a in F:
                if (a*omega + fct.diffn()).is_regular_on_U0():
                    return a*omega + fct.diffn()
        raise ValueError("Unknown.")
        
    def lift_to_de_rham_form(self, eta, threshold = 20):
        '''Given form eta regular on affine part find fct such that eta - d(fct) is regular in infty. (Works for one place at infty now)'''
        from itertools import product
        x_series = self.x_series
        y_series = self.y_series
        z_series = self.z_series
        dx_series = self.dx_series
        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()
        RxyzQ, Rxyz, x, y, z = self.fct_field
        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 = [(eta, eta.expansion_at_infty())]
        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):
                    ff = as_function(self, prod(z[i1]^(k[i1]) for i1 in range(n))/x^i*y^j)
                    ff_exp = ff.diffn().expansion_at_infty()
                    if ff_exp != 0*ff_exp:
                        S += [(ff, ff_exp)]
        for i in range(0, threshold*r):
            for j in range(0, m):
                for k in product(*pr):
                    ff = as_function(self, prod(z[i1]^(k[i1]) for i1 in range(n))*x^i*y^j)
                    ff_exp = ff.diffn().expansion_at_infty()
                    if ff_exp != 0*ff_exp:
                        S += [(ff, ff_exp)]              
        forms = holomorphic_combinations_mixed(S)
        if len(forms) <= self.genus():
            raise ValueError("Increase threshold!")
        result = []
        for ff in forms:
            if ff[0] != 0*self.dx:
                result += [as_cech(self, ff[0], -ff[1])]
        return result
        raise ValueError("Unknown.")

    def de_rham_basis(self, holo_basis = 0, cohomology_basis = 0, threshold = 30):
        if holo_basis == 0:
            holo_basis = self.holomorphic_differentials_basis(threshold = threshold)
        if cohomology_basis == 0:
            cohomology_basis = self.cohomology_of_structure_sheaf_basis(holo_basis = holo_basis, threshold = threshold)
        result = []
        for omega in holo_basis:
            result += [as_cech(self, omega, as_function(self, 0))]
        for f in cohomology_basis:
            omega = self.lift_to_de_rham(f, threshold = threshold)
            result += [as_cech(self, omega, f)]
        return result
    
    def group_action_matrices_holo(AS, basis=0, threshold=10):
        n = AS.height
        if basis == 0:
            basis = AS.holomorphic_differentials_basis(threshold=threshold)
        F = AS.base_ring
        return as_group_action_matrices(F, basis, AS.group.gens, basis = basis)
    
    def group_action_matrices_dR(AS, basis = 0, threshold=8):
        n = AS.height
        if basis == 0:
            holo_basis = AS.holomorphic_differentials_basis(threshold = threshold)
            str_basis = AS.cohomology_of_structure_sheaf_basis(holo_basis = holo_basis, threshold = threshold)
            dr_basis = AS.de_rham_basis(holo_basis = holo_basis, cohomology_basis = str_basis, threshold=threshold)
            basis = [holo_basis, str_basis, dr_basis]
        return as_group_action_matrices(F, basis[2], AS.group.gens, basis = basis)
    
    def group_action_matrices_log_holo(AS):
        n = AS.height
        F = AS.base_ring
        return as_group_action_matrices(F, AS.at_most_poles_forms(1), AS.group.gens, basis = AS.at_most_poles_forms(1))
    
    def riemann_roch_space(self, pole_orders, threshold = 8):
        """ Find fcts with pole order in infty's at most pole_order from the given dictionary. The keys of the dictionary are pairs (place_at_infty, group element). The items are the poles orders at those places.
        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_series
        y_series = self.y_series
        z_series = self.z_series
        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()
        RxyzQ, Rxyz, x, y, z = self.fct_field
        F = C.base_ring
        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_orders[(0, self.group.one)])
        print('iteration')
        for i in range(delta):
            for g in self.fiber(place = i):
                if i!=0 or g != self.group.one:
                    print('iteration')
                    forms = [(omega, omega.group_action(g).expansion_at_infty(place = i)) for omega in forms]
                    forms = holomorphic_combinations_fcts(forms, pole_orders[(i, g)])
        return forms

    def riemann_roch_space_forms(self, pole_orders, threshold = 8):
        """ Find fcts with pole order in infty's at most pole_order from the given dictionary. The keys of the dictionary are pairs (place_at_infty, group element). The items are the poles orders at those places.
        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_series
        y_series = self.y_series
        z_series = self.z_series
        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()
        RxyzQ, Rxyz, x, y, z = self.fct_field
        F = C.base_ring
        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_orders[(0, self.group.one)])
        for i in range(delta):
            for g in self.fiber(place = i):
                if i!=0 or g != self.group.one:
                    forms = [(omega, omega.group_action(g).expansion_at_infty(place = i)) for omega in forms]
                    forms = holomorphic_combinations_forms(forms, pole_orders[(i, g)])
        return forms