p = 5
#F.<a> = GF(p^2, 'a')
F.<a> = GF(p^8, 'a')
#F = GF(p).algebraic_closure()
#a = F.gen(2)
Rx.<x> = PolynomialRing(F)
P1 = superelliptic(x, 1)
m = 4
s = (p^8 - 1)/(p^2 - 1)
b = a^s
C = as_cover(P1, [(P1.x)^(m), b*(P1.x)^(m)], prec=300)
print(C)
A, B = group_action_matrices_dR(C)
print('matrices')
print(magma_module_decomposition(A, B, matrices=False))