p = 2
F.<a> = GF(p^2)
Rx.<x> = PolynomialRing(F)
C = superelliptic(x, 1)
Q = quaternion_cover(C, [C.x^3, a*C.x^3, 0*C.x], prec = 300)
print(Q.genus() == len(Q.holomorphic_differentials_basis()))
I2 = Q.canonical_ideal_polynomials(2)
R = I2[0].parent() #ring, to which the polynomials belong
J2 = R.ideal(I2) #ideal defined by the set I
print(J2.is_prime(), len(R.gens()) - J2.dimension() - 1)
A, B = quaternion_matrices_group_action_canonical_ideal(Q, 2)
print(A * B == B^3*A, A^2 == B^2, A^4 == identity_matrix(A.dimensions()[0]))