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