p = 3
F = GF(3^2, 'a')
Rx.<x> = PolynomialRing(F)
P1 = superelliptic(x^2 + 1, 2)
fct1 = (P1.x)^2
fct2 = fct1 + (P1.x)/(P1.y - P1.x)
fct3 = (P1.x)^4
C = heisenberg_cover(P1, [1/2*fct1, fct2, fct3], prec=300)
print(C)
B = C.holomorphic_differentials_basis()
print("Computed basis")
a1, b1, c1 = heisenberg_group_action_matrices_holo(C, basis = B)