p = 3
m = 2
#F = GF(p)
F = GF(p)
Rx.<x> = PolynomialRing(F)
f = x^5 + 1
g = f(x^p - x)
Cf = superelliptic(f, m)
C = superelliptic(g, m)

#class superelliptic_automorphism:
#    def __init__

def autom(omega):
    C1 = omega.curve
    f = omega.form
    F = C1.base_ring
    Rxy.<x, y> = PolynomialRing(F, 2)
    RxyQ = FractionField(Rxy)
    f = RxyQ(f)
    return superelliptic_form(C1, f(x=x+1, y=y))

def automdR(omega):
    C1 = omega.curve
    om0 = omega.omega0.form
    f = omega.f.function
    F = C1.base_ring
    Rxy.<x, y> = PolynomialRing(F, 2)
    RxyQ = FractionField(Rxy)
    f = RxyQ(f)
    om0 = RxyQ(om0)
    f = f(x=x+1, y=y)
    om0 = om0(x=x+1, y=y)
    f = superelliptic_function(C1, f)
    om0 = superelliptic_form(C1, om0)
    return superelliptic_cech(C1, om0, f)

A = C.de_rham_basis()
gC = C.genus()
M = matrix(C.base_ring, 2*gC, 2*gC)
for i in range(2*gC):
    M[:, i] = vector(automdR(A[i]).coordinates())
print(M, M^p == identity_matrix(2*gC))

N = C.verschiebung_matrix().transpose()
print(magmathis(M, N))