p = 7
m = 2
F = GF(p)
Rx.<x> = PolynomialRing(F)
f = x^3 + 1
C_super = superelliptic(f, m)

Rxy.<x, y> = PolynomialRing(F, 2)
f1 = superelliptic_function(C_super, x^2*y)
f2 = superelliptic_function(C_super, x^3)
AS = as_cover(C_super, [f1, f2], prec=1000)

A, B = group_action_matrices(AS)
n = A.dimensions()[0]
print(A*B == B*A)
print(A^p == identity_matrix(n))
print(B^p == identity_matrix(n))