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

f1 = C_super.x^2*C_super.y
f2 = C_super.x^3
AS = as_cover(C_super, [f1, f2], prec=1000)

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

print(magma_module_decomposition(A, B))