p = 3
m = 2
F = GF(p)
Rx.<x> = PolynomialRing(F)
f = x^3 - x + 1
C = superelliptic(f, m)
C1 = patch(C)
print(C1.crystalline_cohomology_basis())