diff --git a/sage/superelliptic_drw/decomposition_into_g0_g8.sage b/sage/superelliptic_drw/decomposition_into_g0_g8.sage
new file mode 100644
index 0000000..ed85141
--- /dev/null
+++ b/sage/superelliptic_drw/decomposition_into_g0_g8.sage
@@ -0,0 +1,66 @@
+def decomposition_g0_g8(fct, prec = 50):
+    '''Writes fct as a difference g0 - g8 + f, with g0 regular on the affine patch and g8 at the points in infinity
+    and f is combination of basis of H^1(X, OX). Output is (g0, g8, f).'''
+    C = fct.curve
+    g = C.genus()
+    coord = fct.coordinates()
+    nontrivial_part = 0*C.x
+    for i, a in enumerate(C.cohomology_of_structure_sheaf_basis()):
+        nontrivial_part += coord[i]*a
+    fct -= nontrivial_part
+    
+    Fxy, Rxy, x, y = C.fct_field
+    fct = Fxy(fct.function)
+    num = fct.numerator()
+    den = fct.denominator()
+    aux_den = superelliptic_function(C, Rxy(den))
+    g0 = superelliptic_function(C, 0)
+    g8 = superelliptic_function(C, 0)
+    for monomial in num.monomials():
+        aux = superelliptic_function(C, monomial)
+        if aux.expansion_at_infty().valuation() >= aux_den.expansion_at_infty().valuation():
+            g8 -= num.monomial_coefficient(monomial)*aux/aux_den
+        else:
+            g0 += num.monomial_coefficient(monomial)*aux/aux_den
+    return (g0, g8, nontrivial_part)
+
+def decomposition_omega0_omega8(omega, prec=50):
+    '''Writes omega as a difference omega0 - omega8, with omega0 regular on the affine patch and omega8 at the points in infinity.'''
+    C = omega.curve
+    omega.form = reduction(C, omega.form)
+    F = C.base_ring
+    delta = C.nb_of_pts_at_infty
+    m = C.exponent
+    if sum(omega.residue(place = i, prec = 50) for i in range(delta)) != 0:
+        raise ValueError(str(omega) + " has non zero residue!")
+    Fxy, Rxy, x, y = C.fct_field
+    Rx.<x> = PolynomialRing(F)
+    Fx = FractionField(Rx)
+    fct = Fxy(omega.form)
+    num = fct.numerator()
+    den = fct.denominator()
+    aux_den = superelliptic_function(C, Rxy(den))
+    g0 = superelliptic_function(C, 0)
+    g8 = superelliptic_function(C, 0)
+    for j in range(0, m):
+        component = Fx(omega.jth_component(j))
+        q, r = component.numerator().quo_rem(component.denominator())
+        g0 += (C.y)^(-j)*superelliptic_function(C, Rxy(q))
+        if ((C.y)^(-j)*superelliptic_function(C, Fxy(r/component.denominator()))*C.dx).expansion_at_infty().valuation() < 0:
+            raise ValueError("Something went wrong for "+str(omega))
+        g8 -= (C.y)^(-j)*superelliptic_function(C, Fxy(r/component.denominator()))
+    g0, g8  = g0*C.dx, g8*C.dx
+    if g0.is_regular_on_U0():
+        return (g0, g8)
+    #Rx.<x> = PolynomialRing(F)
+    #Rx.<x> = PolynomialRing(F)
+    #aux_fct = (g0.form)*y
+    else:
+        raise ValueError("Something went wrong for "+str(omega) +". Result would be "+str(g0)+ " and " + str(g8))
+        
+        
+def decomposition_g0_g8_pth_power(fct):
+    '''Decompose fct as g0 - g8 + A^p, if possible. Output: (g0, g8, A).'''
+    coor = fct.coordinates()
+    C = fct.curve
+    return 0
\ No newline at end of file
diff --git a/sage/superelliptic_drw/second_patch.sage b/sage/superelliptic_drw/second_patch.sage
new file mode 100644
index 0000000..7a85d09
--- /dev/null
+++ b/sage/superelliptic_drw/second_patch.sage
@@ -0,0 +1,30 @@
+def patch(C):
+    if C.exponent != 2:
+        raise ValueError("Not implemented yet!")
+    Fxy, Rxy, x, y = C.fct_field
+    F = C.base_ring
+    Rx.<x> = PolynomialRing(F)
+    f = C.polynomial
+    g = C.genus()
+    f_star = Rx(x^(2*g+2)*f(1/x))
+    return superelliptic(f_star, 2)
+
+def second_patch(argument):
+    C = argument.curve
+    C1 = patch(C)
+    Fxy, Rxy, x, y = C.fct_field
+    g = C.genus()
+    if isinstance(argument, superelliptic_function):
+        fct = Fxy(argument.function)
+        fct1 = Fxy(fct.subs({x : 1/x, y : y/x^(g+1)}))
+        return superelliptic_function(C1, fct1)
+    if isinstance(argument, superelliptic_form):
+        fct = Fxy(argument.form)
+        fct1 = Fxy(fct.subs({x : 1/x, y : y/x^(g+1)}))
+        fct1 *= -x^(-2)
+        return superelliptic_form(C1, fct1)
+    
+def lift_form_to_drw(omega):
+    A, B = regular_form(omega)
+    A, B = A.change_ring(QQ), B.change_ring(QQ)
+    print("%s dx + %s dy"%(A, B))
\ No newline at end of file
diff --git a/sage/superelliptic_drw/superelliptic_drw.sage b/sage/superelliptic_drw/superelliptic_drw.sage
new file mode 100644
index 0000000..985ea86
--- /dev/null
+++ b/sage/superelliptic_drw/superelliptic_drw.sage
@@ -0,0 +1,419 @@
+class superelliptic_witt:
+    def __init__(self, t, f):
+        ''' Define Witt function on C of the form [t] + V(f). '''
+        self.curve = t.curve
+        C = t.curve
+        p = C.characteristic
+        self.t = t #superelliptic_function
+        self.f = f #superelliptic_function
+    
+    def __repr__(self):
+        f = self.f
+        t = self.t
+        if f.function == 0:
+            return "[" + str(t) + "]"
+        if t.function == 0:
+            return "V(" + str(f) + ")"
+        return "[" + str(t) + "] + V(" + str(f) + ")"
+    
+    def __neg__(self):
+        f = self.f
+        t = self.t
+        return superelliptic_witt(-t, -f)
+    
+    def __add__(self, other):
+        C = self.curve
+        second_coor = 0*C.x
+        X = self.t
+        Y = other.t
+        for i in range(1, p):
+            second_coor -= binomial_prim(p, i)*X^i*Y^(p-i)
+        return superelliptic_witt(self.t + other.t, self.f + other.f + second_coor)
+    
+    def __sub__(self, other):
+        return self + (-other)
+    
+    def __rmul__(self, other):
+        p = self.curve.characteristic
+        if other in ZZ:
+            if other == 0:
+                return superelliptic_witt(0*C.x, 0*C.x)
+            if other > 0:
+                return self + (other-1)*self
+            if other < 0:
+                return (-other)*(-self)
+        if other in QQ:
+            other_integer = Integers(p^2)(other)
+            return other_integer*self
+    
+    def __mul__(self, other):
+        C = self.curve
+        p = C.characteristic
+        if isinstance(other, superelliptic_witt):
+            t1 = self.t
+            f1 = self.f
+            t2 = other.t
+            f2 = other.f
+            return superelliptic_witt(t1*t2, t1^p*f2 + t2^p*f1)
+        if isinstance(other, superelliptic_drw_form):
+            h1 = other.h1
+            h2 = other.h2
+            omega = other.omega
+            t = self.t
+            f = self.f
+            aux_form = t^p*omega - h2*t^(p-1)*t.diffn() + f*h1^p*(C.x)^(p-1)*C.dx
+            return superelliptic_drw_form(t*h1, aux_form, t^p*h2)
+    
+    def __eq__(self, other):
+        return self.t == other.t and self.f == other.f
+    
+    def diffn(self, dy_w = 0):
+        if dy_w == 0:
+            dy_w = self.curve.dy_w()
+        C = self.curve
+        t = self.t
+        f = self.f
+        fC = C.polynomial
+        F = C.base_ring
+        Rxy.<x, y> = PolynomialRing(F, 2)
+        if t.function == 0:
+            return superelliptic_drw_form(0*C.x, 0*C.dx, f)
+        t_polynomial = t.function
+        num = t_polynomial.numerator()
+        den = t_polynomial.denominator()
+        num_t_fct = superelliptic_function(C, num)
+        den_t_fct = superelliptic_function(C, den)
+        inv_den_t_fct = superelliptic_function(C, 1/den)
+        if den != 1:
+            # d([N/D] + V(f)) = [1/D]*d([N]) - [N]*[D^(-2)]*d([D]) + dV(f)
+            return ((den_t_fct)^(-1)).teichmuller()*num_t_fct.teichmuller().diffn() - ((den_t_fct)^(-2)).teichmuller()*num_t_fct.teichmuller()*den_t_fct.teichmuller().diffn() + superelliptic_drw_form(0*C.x, 0*C.dx, f)
+        t_polynomial = Rxy(t_polynomial)
+        M = t_polynomial.monomials()[0]
+        a = t_polynomial.monomial_coefficient(M)
+        #[P] = [aM] + Q, where Q = ([P] - [aM])
+        aM_fct = superelliptic_function(C, a*M)
+        Q = self - aM_fct.teichmuller()
+        exp_x = M.exponents()[0][0]
+        exp_y = M.exponents()[0][1]
+        return Q.diffn() + exp_x*superelliptic_drw_form(aM_fct/C.x, 0*C.dx, 0*C.x) + exp_y*(aM_fct/C.y).teichmuller()*dy_w
+
+def binomial_prim(p, i):
+    return binomial(p, i)/p
+    
+def reduce_rational_fct(fct, p):
+    Rxy.<x, y> = PolynomialRing(QQ)
+    Fxy = FractionField(Rxy)
+    fct = Fxy(fct)
+    num = Rxy(fct.numerator())
+    den = Rxy(fct.denominator())
+    num1 = Rxy(0)
+    for m in num.monomials():
+        a = num.monomial_coefficient(m)
+        num1 += (a%p^2)*m
+    den1 = Rxy(0)
+    for m in den.monomials():
+        a = den.monomial_coefficient(m)
+        den1 += (a%p^2)*m
+    return num1/den1
+
+def teichmuller(fct):
+    C = fct.curve
+    return superelliptic_witt(fct, 0*C.x)
+
+superelliptic_function.teichmuller = teichmuller
+
+#dy = [f(x)]'/2*y dx
+#[f1 + M] = [f1] + [M] + V(cos)
+#d[f1 + M] = d[f1] + d[M] + dV(f1*M)
+#M = b x^a
+#d[M] = a*[b x^(a-1)]
+
+def auxilliary_derivative(P):
+    '''Return "derivative" of P, where P depends only on x. In other words d[P(x)].'''
+    P0 = P.t.function
+    P1 = P.f.function
+    C = P.curve
+    F = C.base_ring
+    Rx.<x> = PolynomialRing(F)
+    P0 = Rx(P0)
+    P1 = Rx(P1)
+    if P0 == 0:
+        return superelliptic_drw_form(0*C.x, 0*C.dx, P.f)
+    M = P0.monomials()[0]
+    a = P0.monomial_coefficient(M)
+    #[P] = [aM] + Q, where Q = ([P] - [aM])
+    aM_fct = superelliptic_function(C, a*M)
+    Q = P - aM_fct.teichmuller()
+    exp = M.exponents()[0]
+    return auxilliary_derivative(Q) + exp*superelliptic_drw_form(aM_fct/C.x, 0*C.dx, 0*C.x)
+    
+class superelliptic_drw_form:
+    def __init__(self, h1, omega, h2):
+        '''Form [h1] d[x] + V(omega) + dV([h])'''
+        self.curve = h1.curve
+        self.h1 = h1
+        self.omega = omega
+        self.h2 = h2
+        
+    def r(self):
+        C = self.curve
+        h1 = self.h1
+        return superelliptic_form(C, h1.function)
+        
+    def __eq__(self, other):
+        eq1 = (self.h1 == self.h1)
+        try:
+            H = (self.h2 - other.h2).pth_root()
+        except:
+            return False
+        eq2 = ((self.omega - other.omega).cartier() - H.diffn()) == 0*self.curve.dx
+        if eq1 and eq2:
+            return True
+        return False
+
+    def __repr__(self):
+        C = self.curve
+        h1 = self.h1
+        omega = self.omega
+        h2 = self.h2
+        result = ""
+        if h1.function != 0:
+            result += "[" + str(h1) + "] d[x]"
+        if h1.function !=0 and omega.form != 0:
+            result += " + "
+        if omega.form != 0:
+            result += "V(" + str(omega) + ")"
+        if h2.function !=0 and omega.form != 0:
+            result += " + "
+        if h2.function != 0:
+            result += "dV([" + str(h2) +"])"
+        if result == "":
+            result += "0"
+        return result
+
+    def __rmul__(self, other):
+        h1 = self.h1
+        h2 = self.h2
+        omega = self.omega
+        p = self.curve.characteristic
+        if other in ZZ:
+            if other == 0:
+                return superelliptic_drw_form(0*C.x, 0*C.dx, 0*C.x)
+            if other > 0:
+                return self + (other-1)*self
+            if other < 0:
+                return (-other)*(-self)
+        if other in QQ:
+            other_integer = Integers(p^2)(other)
+            return other_integer*self
+        t = other.t
+        f = other.f
+        aux_form = t^p*omega - h2*t^(p-1)*t.diffn() + f*h1^p*(C.x)^(p-1)*C.dx
+        return superelliptic_drw_form(t*h1, aux_form, t^p*h2)
+
+    def __neg__(self):
+        C = self.curve
+        h1 = self.h1
+        h2 = self.h2
+        omega = self.omega
+        return superelliptic_drw_form(-h1, -omega, -h2)
+    
+    def __sub__(self, other):
+        return self + (-other)
+    
+    def __add__(self, other):
+        C = self.curve
+        h1 = self.h1
+        h2 = self.h2
+        omega = self.omega
+        H1 = other.h1
+        H2 = other.h2
+        OMEGA = other.omega
+        aux = (teichmuller(h1) + teichmuller(H1))*superelliptic_drw_form(C.one, 0*C.dx, 0*C.x)
+        h1_aux = aux.h1
+        h2_aux = aux.h2
+        omega_aux = aux.omega
+        return superelliptic_drw_form(h1_aux, omega + OMEGA + omega_aux, h2 + H2 + h2_aux)
+    
+    def frobenius(self):
+        C = self.curve
+        h1 = self.h1
+        h2 = self.h2
+        p = C.characteristic
+        return h1^p*C.x^(p-1)*C.dx + h2.diffn()
+
+def mult_by_p(omega):
+    C = omega.curve
+    fct = omega.form
+    Fxy, Rxy, x, y = C.fct_field
+    omega = superelliptic_form(C, fct^p * x^(p-1))
+    result = superelliptic_drw_form(0*C.x, omega, 0*C.x)
+    return result
+
+def verschiebung(elt):
+    C = elt.curve
+    if isinstance(elt, superelliptic_function):
+        return superelliptic_witt(0*C.x, elt)
+
+    if isinstance(elt, superelliptic_form):
+        return superelliptic_drw_form(0*C.x, elt, 0*C.x)
+
+superelliptic_form.verschiebung = verschiebung
+superelliptic_function.verschiebung = verschiebung
+
+
+class superelliptic_drw_cech:
+    def __init__(self, omega0, f):
+        self.curve = omega0.curve
+        self.omega0 = omega0
+        self.omega8 = omega0 - f.diffn()
+        self.f = f
+        
+        
+    def reduce(self):
+        C = self.curve
+        fct = self.f
+        f_first_comp = fct.t
+        f_second_comp = fct.f
+        decomp_first_comp = decomposition_g0_g8(f_first_comp)
+        decomp_second_comp = decomposition_g0_g8(f_second_comp)
+        new = self
+        new.omega0 -= decomposition_g0_g8(f_first_comp)[0].teichmuller().diffn()
+        new.omega0 -= decomposition_g0_g8(f_second_comp)[0].verschiebung().diffn()
+        new.f = decomposition_g0_g8(f_first_comp)[2].teichmuller() + decomposition_g0_g8(f_second_comp)[2].verschiebung()
+        new.omega8 = new.omega0 - new.f.diffn()
+        return new
+        
+    def __repr__(self):
+        return("(" + str(self.omega0) + ", "+ str(self.f) + ", " + str(self.omega8) + ")")
+
+    def __add__(self, other):
+        C = self.curve
+        omega0 = self.omega0
+        f = self.f
+        omega0_1 = other.omega0
+        f_1 = other.f
+        return superelliptic_drw_cech(omega0 + omega0_1, f + f_1)
+    
+    def __sub__(self, other):
+        C = self.curve
+        omega0 = self.omega0
+        f = self.f
+        omega0_1 = other.omega0
+        f_1 = other.f
+        return superelliptic_drw_cech(omega0 - omega0_1, f - f_1)
+        
+    def __neg__(self):
+        C = self.curve
+        omega0 = self.omega0
+        f = self.f
+        return superelliptic_drw_cech(-omega0, -f)
+    
+    def __rmul__(self, other):
+        omega0 = self.omega0
+        f = self.f
+        return superelliptic_drw_cech(other*omega0, other*f)
+    
+    def r(self):
+        omega0 = self.omega0
+        f = self.f
+        C = self.curve
+        return superelliptic_cech(C, omega0.h1*C.dx, f.t)
+    
+    def coordinates(self, basis = 0):
+        C = self.curve
+        g = C.genus()
+        coord_mod_p = self.r().coordinates()
+        print(coord_mod_p)
+        coord_lifted = [lift(a) for a in coord_mod_p]
+        if basis == 0:
+            basis = C.crystalline_cohomology_basis()
+        aux = self
+        for i, a in enumerate(basis):
+            aux -= coord_lifted[i]*a
+        print('aux before reduce', aux)
+        #aux = aux.reduce() # Now aux = p*cech class.
+        # Now aux should be of the form (V(smth), V(smth), V(smth))
+        print('aux V(smth)', aux)
+        aux_divided_by_p = superelliptic_cech(C, aux.omega0.omega.cartier(), aux.f.f.pth_root())
+        print('aux.omega0.omega.cartier()', aux.omega0.omega.cartier())
+        coord_aux_divided_by_p = aux_divided_by_p.coordinates()
+        coord_aux_divided_by_p = [ZZ(a) for a in coord_aux_divided_by_p]
+        coordinates = [ (coord_lifted[i] + p*coord_aux_divided_by_p[i])%p^2 for i in range(2*g)]
+        return coordinates
+    
+    def is_regular(self):
+        print(self.omega0.r().is_regular_on_U0(), self.omega8.r().is_regular_on_Uinfty(), self.omega0.frobenius().is_regular_on_U0(),     self.omega8.frobenius().is_regular_on_Uinfty())
+        eq1 = self.omega0.r().is_regular_on_U0() and self.omega8.r().is_regular_on_Uinfty()
+        eq2 = self.omega0.frobenius().is_regular_on_U0() and self.omega8.frobenius().is_regular_on_Uinfty()
+        return eq1 and eq2
+        
+
+def de_rham_witt_lift(cech_class, prec = 50):
+    C = cech_class.curve
+    g = C.genus()
+    omega0 = cech_class.omega0
+    omega8 = cech_class.omega8
+    fct = cech_class.f
+    omega0_regular = regular_form(omega0) #Present omega0 in the form P dx + Q dy
+    print('omega0_regular', omega0_regular)
+    omega0_lift = omega0_regular[0].teichmuller()*(C.x.teichmuller().diffn()) + omega0_regular[1].teichmuller()*(C.y.teichmuller().diffn())
+    #Now the obvious lift of omega0 = P dx + Q dy to de Rham-Witt is [P] d[x] + [Q] d[y]
+    print('omega8', omega8, 'second_patch(omega8)', second_patch(omega8))
+    omega8_regular = regular_form(second_patch(omega8)) # The same for omega8.
+    print('omega8_regular 1', omega8_regular)
+    omega8_regular = (second_patch(omega8_regular[0]), second_patch(omega8_regular[1]))
+    print('omega8_regular 2', omega8_regular)
+    u = (C.x)^(-1)
+    v = (C.y)/(C.x)^(g+1)
+    omega8_lift = omega8_regular[0].teichmuller()*(u.teichmuller().diffn()) + omega8_regular[1].teichmuller()*(v.teichmuller().diffn())
+    aux = omega0_lift - omega8_lift - fct.teichmuller().diffn() # now aux is of the form (V(smth) + dV(smth), V(smth))
+    return aux
+    if aux.h1.function != 0:
+        raise ValueError('Something went wrong - aux is not of the form (V(smth) + dV(smth), V(smth)).')
+    decom_aux_h2 = decomposition_g0_g8(aux.h2, prec=prec) #decompose dV(smth) in aux as smth regular on U0 - smth regular on U8.
+    aux_h2 = decom_aux_h2[0]
+    aux_f = decom_aux_h2[2]
+    aux_omega0 = decomposition_omega0_omega8(aux.omega, prec=prec)[0]
+    result = superelliptic_drw_cech(omega0_lift - aux_h2.verschiebung().diffn() - aux_omega0.verschiebung(), fct.teichmuller() + aux_f.verschiebung())
+    return result.reduce()
+
+def crystalline_cohomology_basis(self, prec = 50):
+    result = []
+    for a in self.de_rham_basis():
+        result += [de_rham_witt_lift(a, prec = prec)]
+    return result
+
+superelliptic.crystalline_cohomology_basis = crystalline_cohomology_basis
+
+def autom(self):
+    C = self.curve
+    F = C.base_ring
+    Rxy.<x, y> = PolynomialRing(F, 2)
+    Fxy = FractionField(Rxy)
+    if isinstance(self, superelliptic_function):
+        result = superelliptic_function(C, Fxy(self.function).subs({x:x+1, y:y}))
+        return result
+    if isinstance(self, superelliptic_form):
+        result = superelliptic_form(C, Fxy(self.form).subs({x:x+1, y:y}))
+        return result
+    if isinstance(self, superelliptic_witt):
+        result = superelliptic_witt(autom(self.t), autom(self.f))
+        return result
+    if isinstance(self, superelliptic_drw_form):
+        result = superelliptic_drw_form(0*C.x, autom(self.omega), autom(self.h2))
+        result += autom(self.h1).teichmuller()*(C.x + C.one).teichmuller().diffn()
+        return result
+    if isinstance(self, superelliptic_drw_cech):
+        result = superelliptic_drw_cech(autom(self.omega0), autom(self.f))
+        return result
+    
+    
+def dy_w(C):
+    '''Return d[y].'''
+    fC = C.polynomial
+    fC = superelliptic_function(C, fC)
+    fC = fC.teichmuller()
+    dy_w = 1/2* ((C.y)^(-1)).teichmuller()*auxilliary_derivative(fC)
+    return dy_w
+superelliptic.dy_w = dy_w
\ No newline at end of file