p = 5
m = 1
F = GF(p)
Rx.<x> = PolynomialRing(F)
f = x
C_super = superelliptic(f, m)

Rxy.<x, y> = PolynomialRing(F, 2)
for a in range(3, 13):
    for b in range(3, 13):
        if a %p != 0 and b%p != 0 and a !=b:
            try:
                f1 = superelliptic_function(C_super, x^a+x)
                f2 = superelliptic_function(C_super, x^b)
                AS = as_cover(C_super, [f1, f2], prec=1000)
                #print(AS.at_most_poles(AS.exponent_of_different_prim()))
                print(AS.magical_element(threshold = 20))
            except:
                pass