p = 2
F = GF(p)
Rx.<x> = PolynomialRing(F)
C = superelliptic(x, 1)
AS = as_cover(C, [C.x^3, C.x^5], prec = 200)
B = AS.holo_polydifferentials_basis(3)
print(len(B) == (2 * 3 - 1) * (AS.genus() - 1)) #is the dimension as predicted by Riemann--Roch?
I2 = AS.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)
I3 = AS.canonical_ideal_polynomials(3)
J3 = R.ideal(I2 + I3)
print(J3.is_prime(), len(R.gens()) - J3.dimension() - 1)

p = 5
F = GF(p)
Rx.<x> = PolynomialRing(F)
C = superelliptic(x, 1)
AS = as_cover(C, [C.x^8], prec = 200)
I2 = AS.canonical_ideal_polynomials(2, threshold=15)
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)
I3 = AS.canonical_ideal_polynomials(3, threshold=20)
J3 = R.ideal(I2 + I3)
print(J3.is_prime(), len(R.gens()) - J3.dimension() - 1)