diff --git a/as_covers/as_cover_class.sage b/as_covers/as_cover_class.sage
index 20e05ea..f974ddc 100644
--- a/as_covers/as_cover_class.sage
+++ b/as_covers/as_cover_class.sage
@@ -135,11 +135,11 @@ class as_cover:
                     eta_exp = eta.expansion(pt=self.branch_points[0])
                     S += [(eta, eta_exp)]
 
-        forms = as_holomorphic_combinations(S)
+        forms = holomorphic_combinations(S)
         
         for pt in self.branch_points[1:]:
             forms = [(omega, omega.expansion(pt=pt)) for omega in forms]
-            forms = as_holomorphic_combinations(forms)
+            forms = holomorphic_combinations(forms)
         
         if len(forms) < self.genus():
             print("I haven't found all forms, only ", len(forms), " of ", self.genus())
@@ -376,7 +376,7 @@ class as_cover:
                     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 = as_holomorphic_combinations(S)
+        forms = holomorphic_combinations(S)
         if len(forms) <= self.genus():
             raise ValueError("Increase threshold!")
         for omega in forms:
@@ -398,36 +398,3 @@ class as_cover:
             result += [as_cech(self, omega, f)]
         return result
     
-def as_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
\ No newline at end of file
diff --git a/as_covers/as_form_class.sage b/as_covers/as_form_class.sage
index bcf2962..0b20167 100644
--- a/as_covers/as_form_class.sage
+++ b/as_covers/as_form_class.sage
@@ -188,75 +188,6 @@ def are_forms_linearly_dependent(set_of_forms):
     denominators = prod(denominator(omega.form) for omega in set_of_forms)
     return is_linearly_dependent([Rxyz(denominators*omega.form) for omega in set_of_forms])
 
-def holomorphic_combinations_fcts(S, pole_order):
-    '''given a set S of (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([Rt(g[1]).valuation() for g in S])
-    if minimal_valuation >= -pole_order:
-        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 + Rt(eta_exp).valuation()
-        list_coeffs = a*[0] + Rt(eta_exp).list() + (-minimal_valuation)*[0]
-        list_coeffs = list_coeffs[:-minimal_valuation - pole_order]
-        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_function(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
-
-#given a set S of (form, corresponding Laurent series at some pt), find their combinations holomorphic at that pt
-def holomorphic_combinations_forms(S, pole_order):
-    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([Rt(g[1]).valuation() for g in S])
-    if minimal_valuation >= -pole_order:
-        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 + Rt(eta_exp).valuation()
-        list_coeffs = a*[0] + Rt(eta_exp).list() + (-minimal_valuation)*[0]
-        list_coeffs = list_coeffs[:-minimal_valuation - pole_order]
-        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
-
-#print only forms that are log at the branch pts, but not holomorphic
 def only_log_forms(C_AS):
     list1 = AS.at_most_poles_forms(0)
     list2 = AS.at_most_poles_forms(1)
diff --git a/as_covers/as_transform.sage b/as_covers/as_transform.sage
index 6aec5d9..c1506c3 100644
--- a/as_covers/as_transform.sage
+++ b/as_covers/as_transform.sage
@@ -2,7 +2,6 @@ def artin_schreier_transform(power_series, prec = 10):
     """Given a power_series, find correction such that power_series - (correction)^p +correction has valuation
        -jump non divisible by p. Also, express t (the variable) in terms of the uniformizer at infty on the curve
        z^p - z = power_series, where z = 1/t_new^(jump) and express z in terms of the new uniformizer."""
-    print('as transform for ', power_series)
     correction = 0
     F = power_series.parent().base()
     p = F.characteristic()
@@ -16,7 +15,6 @@ def artin_schreier_transform(power_series, prec = 10):
         aux = t^p - t
         z = new_reverse(aux, prec = prec)
         z = z(t = power_series)
-        print("a")
         return(0, 0, t, z)
     
     while(power_series.valuation() % p == 0 and power_series.valuation() < 0):
@@ -41,5 +39,4 @@ def artin_schreier_transform(power_series, prec = 10):
         z = new_reverse(aux, prec = prec)
         z = z(t = power_series)
         z = z + correction
-        print('corr', correction, 'ps', power_series, 'z', z)
         return(0, correction, t, z)
\ No newline at end of file
diff --git a/as_covers/holomorphic_combinations.sage b/as_covers/holomorphic_combinations.sage
new file mode 100644
index 0000000..b5fcf34
--- /dev/null
+++ b/as_covers/holomorphic_combinations.sage
@@ -0,0 +1,106 @@
+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 = 0*C_AS.dx
+        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
+
+def holomorphic_combinations_fcts(S, pole_order):
+    '''given a set S of (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([Rt(g[1]).valuation() for g in S])
+    if minimal_valuation >= -pole_order:
+        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 + Rt(eta_exp).valuation()
+            if eta_exp !=0:
+                list_coeffs = a*[0] + Rt(eta_exp).list() + (-minimal_valuation)*[0]
+                list_coeffs = list_coeffs[:-minimal_valuation - pole_order]
+            else:
+                list_coeffs = (-minimal_valuation - pole_order)*[0]
+            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 = 0*C_AS.x
+        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
+
+def holomorphic_combinations_forms(S, pole_order):
+    '''given a set S of (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([Rt(g[1]).valuation() for g in S])
+    if minimal_valuation >= -pole_order:
+        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 + Rt(eta_exp).valuation()
+        list_coeffs = a*[0] + Rt(eta_exp).list() + (-minimal_valuation)*[0]
+        list_coeffs = list_coeffs[:-minimal_valuation - pole_order]
+        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 = 0*C_AS.dx
+        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
+
+#print only forms that are log at the branch pts, but not holomorphic
\ No newline at end of file
diff --git a/heisenberg_covers/dual_element.sage b/heisenberg_covers/dual_element.sage
new file mode 100644
index 0000000..2b81a29
--- /dev/null
+++ b/heisenberg_covers/dual_element.sage
@@ -0,0 +1,28 @@
+def heisenberg_dual_elt(AS, zmag):
+    '''Find the trace dual of a given elt zmag in the function field of an Artin-Schreier cover AS.'''
+    p = AS.characteristic
+    n = 3
+    G = AS.group
+    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)
+    M = matrix(RxyzQ, p^n, p^n)
+    for i in range(p^n):
+        for j in range(p^n):
+            print('creating matrix', i, j)
+            elt = tr2((zmag.group_action(G[i])*zmag.group_action(G[j]))).function
+            elt = Rxyz(elt.numerator())/Rxyz(elt.denominator())
+            M[i, j] = RxyzQ(elt)
+    main_det = -1 #M.determinant()
+    zvee = as_function(AS, 0)
+    for i in range(p^n):
+        print(i)
+        Mprim = matrix(RxyzQ, M)
+        Mprim[:, i] = vector([(j == 0) for j in range(p^n)])
+        fi = Mprim.determinant()/main_det
+        zvee += heisenberg_function(AS, fi*(zmag.group_action(G[i]).function))
+    return zvee
\ No newline at end of file
diff --git a/heisenberg_covers/heisenberg_covers.sage b/heisenberg_covers/heisenberg_covers.sage
index ef58250..e3348f9 100644
--- a/heisenberg_covers/heisenberg_covers.sage
+++ b/heisenberg_covers/heisenberg_covers.sage
@@ -137,14 +137,14 @@ class heisenberg_cover:
                     eta_exp = eta.expansion(pt=self.branch_points[0])
                     S += [(eta, eta_exp)]
 
-        forms = heisenberg_holomorphic_combinations(S)
+        forms = holomorphic_combinations(S)
         print("next iteration")
         for pt in self.branch_points:
-            for g in [(0, i, 0) for i in range(p)]:
+            for g in self.fiber(place = pt):
                 print("next iteration")
                 if pt != self.branch_points[0] or g != (0, 0, 0):
                     forms = [(omega, omega.group_action(g).expansion(pt=pt)) for omega in forms]
-                    forms = heisenberg_holomorphic_combinations(forms)
+                    forms = holomorphic_combinations(forms)
                     
         return forms
         if len(forms) < self.genus():
@@ -422,7 +422,8 @@ class heisenberg_cover:
                 result += [g]
         return result
     
-    def stabilizer_coset_reps(self, place = 0):
+    def fiber(self, place = 0):
+        'Gives representatives for the quotient G/G_P for given place. Those are in bijection with the fiber.'
         result = [(0, 0, 0)]
         p = self.characteristic
         H = self.stabilizer(place = place)
@@ -466,46 +467,11 @@ def at_most_poles2(self, pole_orders, threshold = 8):
                 S += [(eta, eta_exp)]
 
     forms = holomorphic_combinations_fcts(S, pole_orders[(0, (0, 0, 0))])
-
+    print('iteration')
     for i in range(delta):
-        for g in self.stabilizer_coset_reps(place = i):
+        for g in self.fiber(place = i):
             if i!=0 or g != (0, 0, 0):
-                forms = [(omega, omega.expansion_at_infty(place = i)) for omega in forms]
+                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 heisenberg_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 = heisenberg_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
\ No newline at end of file
diff --git a/heisenberg_covers/heisenberg_form_class.sage b/heisenberg_covers/heisenberg_form_class.sage
index 8448b36..4fc4c03 100644
--- a/heisenberg_covers/heisenberg_form_class.sage
+++ b/heisenberg_covers/heisenberg_form_class.sage
@@ -213,74 +213,6 @@ def are_forms_linearly_dependent(set_of_forms):
     denominators = prod(denominator(omega.form) for omega in set_of_forms)
     return is_linearly_dependent([Rxyz(denominators*omega.form) for omega in set_of_forms])
 
-def holomorphic_combinations_fcts(S, pole_order):
-    '''given a set S of (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([Rt(g[1]).valuation() for g in S])
-    if minimal_valuation >= -pole_order:
-        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 + Rt(eta_exp).valuation()
-        list_coeffs = a*[0] + Rt(eta_exp).list() + (-minimal_valuation)*[0]
-        list_coeffs = list_coeffs[:-minimal_valuation - pole_order]
-        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 = heisenberg_function(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
-
-#given a set S of (form, corresponding Laurent series at some pt), find their combinations holomorphic at that pt
-def holomorphic_combinations_forms(S, pole_order):
-    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([Rt(g[1]).valuation() for g in S])
-    if minimal_valuation >= -pole_order:
-        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 + Rt(eta_exp).valuation()
-        list_coeffs = a*[0] + Rt(eta_exp).list() + (-minimal_valuation)*[0]
-        list_coeffs = list_coeffs[:-minimal_valuation - pole_order]
-        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 = heisenberg_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
-
 #print only forms that are log at the branch pts, but not holomorphic
 def only_log_forms(C_AS):
     list1 = AS.at_most_poles_forms(0)
diff --git a/heisenberg_covers/heisenberg_function_class.sage b/heisenberg_covers/heisenberg_function_class.sage
index f20a520..e578757 100644
--- a/heisenberg_covers/heisenberg_function_class.sage
+++ b/heisenberg_covers/heisenberg_function_class.sage
@@ -148,6 +148,7 @@ class heisenberg_function:
             return self.group_action(elt1).group_action((0, 0, 1))
 
     def trace(self):
+        print('trace ')
         C = self.curve
         C_super = C.quotient
         n = C.height
@@ -166,7 +167,8 @@ class heisenberg_function:
         result = result.function
         Rxy.<x, y> = PolynomialRing(F, 2)
         Qxy = FractionField(Rxy)
-        result = heisenberg_reduction(AS, result)
+        result = heisenberg_reduction(C, result)
+        #print(result)
         return superelliptic_function(C_super, Qxy(result))
     
     def coordinates(self, prec = 100, basis = 0):
diff --git a/heisenberg_covers/heisenberg_group.sage b/heisenberg_covers/heisenberg_group.sage
new file mode 100644
index 0000000..f903a4e
--- /dev/null
+++ b/heisenberg_covers/heisenberg_group.sage
@@ -0,0 +1,23 @@
+def heisenberg_mult(v1, v2, p):
+    i1, j1, k1 = v1
+    i2, j2, k2 = v2
+    i3 = (i1 + i2)%p
+    j3 = (j1 + j2)%p
+    k3 = (-i2*j1 + k1 + k2)%p
+    return(i3, j3, k3)
+    
+def heisenberg_inv(v1, p):
+    i1, j1, k1 = v1
+    i2 = p-i1
+    j2 = p-j1
+    k2 = p-k1
+    k2 = (k2 - i2*j2)%p
+    return (i2, j2, k2)
+
+def heisenberg_power(v1, n, p):
+    i1, j1, k1 = v1
+    if n < 0:
+        return heisenberg_inv(heisenberg_power((i1, j1, k1), -n, p), p)
+    if n == 0:
+        return (0, 0, 0)
+    return heisenberg_mult((i1, j1, k1), heisenberg_power((i1, j1, k1), n-1, p), p)
\ No newline at end of file
diff --git a/heisenberg_covers/heisenberg_reduction.sage b/heisenberg_covers/heisenberg_reduction.sage
index 8665c2d..d98cb44 100644
--- a/heisenberg_covers/heisenberg_reduction.sage
+++ b/heisenberg_covers/heisenberg_reduction.sage
@@ -10,7 +10,7 @@ def heisenberg_reduction(AS, fct):
     fct1 = numerator(fct)
     fct2 = denominator(fct)
     denom = heisenberg_function(AS, fct2)
-    denom_norm = prod(heisenberg_function(AS, fct2).group_action(g) for g in AS.group if list(g) != n*[0])
+    denom_norm = prod(heisenberg_function(AS, fct2).group_action(g) for g in AS.group if g != (0, 0, 0))
     fct1 = Rxyz(fct1*denom_norm.function)
     fct2 = Rxyz(fct2*denom_norm.function)
     if fct2 != 1:
diff --git a/heisenberg_covers/ith_magical_component.sage b/heisenberg_covers/ith_magical_component.sage
new file mode 100644
index 0000000..c67ba0d
--- /dev/null
+++ b/heisenberg_covers/ith_magical_component.sage
@@ -0,0 +1,6 @@
+def heisenberg_ith_magical_component(omega, zvee, g):
+    '''Given a form omega on AS cover, element g of group AS.group and normal basis element zmag, find the decomposition
+    sum_g g(zmag) omega_g and return omega_g.'''
+    z_vee_g = zvee.group_action(g)
+    new_form = z_vee_g*omega
+    return new_form.trace()
\ No newline at end of file
diff --git a/init.sage b/init.sage
index 7858480..73b57d2 100644
--- a/init.sage
+++ b/init.sage
@@ -4,6 +4,7 @@ load('superelliptic/superelliptic_form_class.sage')
 load('superelliptic/superelliptic_cech_class.sage')
 load('superelliptic/frobenius_kernel.sage')
 load('as_covers/as_transform.sage')
+load('as_covers/holomorphic_combinations.sage')
 load('as_covers/as_cover_class.sage')
 load('as_covers/as_function_class.sage')
 load('as_covers/as_form_class.sage')
diff --git a/quaternion_covers/quaternion_covers.sage b/quaternion_covers/quaternion_covers.sage
index 92535b7..70bc0ca 100644
--- a/quaternion_covers/quaternion_covers.sage
+++ b/quaternion_covers/quaternion_covers.sage
@@ -137,11 +137,11 @@ class quaternion_cover:
                     eta_exp = eta.expansion(pt=self.branch_points[0])
                     S += [(eta, eta_exp)]
 
-        forms = quaternion_holomorphic_combinations(S)
+        forms = holomorphic_combinations(S)
         
         for pt in self.branch_points[1:]:
             forms = [(omega, omega.expansion(pt=pt)) for omega in forms]
-            forms = quaternion_holomorphic_combinations(forms)
+            forms = holomorphic_combinations(forms)
         
         if len(forms) < self.genus():
             print("I haven't found all forms, only ", len(forms), " of ", self.genus())
@@ -378,7 +378,7 @@ class quaternion_cover:
                     eta = quaternion_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 = quaternion_holomorphic_combinations(S)
+        forms = holomorphic_combinations(S)
         if len(forms) <= self.genus():
             raise ValueError("Increase threshold!")
         for omega in forms:
@@ -398,38 +398,4 @@ class quaternion_cover:
         for f in cohomology_basis:
             omega = self.lift_to_de_rham(f, threshold = threshold)
             result += [quaternion_cech(self, omega, f)]
-        return result
-    
-def quaternion_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 = quaternion_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
\ No newline at end of file
+        return result
\ No newline at end of file
diff --git a/quaternion_covers/quaternion_form_class.sage b/quaternion_covers/quaternion_form_class.sage
index 1ccccca..465f539 100644
--- a/quaternion_covers/quaternion_form_class.sage
+++ b/quaternion_covers/quaternion_form_class.sage
@@ -207,74 +207,6 @@ def are_forms_linearly_dependent(set_of_forms):
     denominators = prod(denominator(omega.form) for omega in set_of_forms)
     return is_linearly_dependent([Rxyz(denominators*omega.form) for omega in set_of_forms])
 
-def holomorphic_combinations_fcts(S, pole_order):
-    '''given a set S of (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([Rt(g[1]).valuation() for g in S])
-    if minimal_valuation >= -pole_order:
-        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 + Rt(eta_exp).valuation()
-        list_coeffs = a*[0] + Rt(eta_exp).list() + (-minimal_valuation)*[0]
-        list_coeffs = list_coeffs[:-minimal_valuation - pole_order]
-        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 = quaternion_function(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
-
-#given a set S of (form, corresponding Laurent series at some pt), find their combinations holomorphic at that pt
-def holomorphic_combinations_forms(S, pole_order):
-    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([Rt(g[1]).valuation() for g in S])
-    if minimal_valuation >= -pole_order:
-        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 + Rt(eta_exp).valuation()
-        list_coeffs = a*[0] + Rt(eta_exp).list() + (-minimal_valuation)*[0]
-        list_coeffs = list_coeffs[:-minimal_valuation - pole_order]
-        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 = quaternion_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
-
 #print only forms that are log at the branch pts, but not holomorphic
 def only_log_forms(C_AS):
     list1 = AS.at_most_poles_forms(0)