p = 3
F = GF(3)
Rx.<x> = PolynomialRing(F)
P1 = superelliptic(x^2 + 1, 1)
fct1 = (P1.x)^2
fct2 = fct1 + (P1.x)/(P1.y - P1.x)
fct3 = (P1.x)^4
C = heisenberg_cover(P1, [fct1, fct2, fct3], prec=300)
print(C)
a, b, c = heisenberg_group_action_matrices_holo(C)