# an implementation of algorithm 1


def does_p_divide_class_number(p,K,N=100):
    m=K.disc().sqrt()
    if ZZ((m-1)/3) % p == 0 or m % p == 0:
        return true # cannot apply criterion
    if not cyclotomic_polynomial((m-1)//2).change_ring(GF(p)).is_irreducible():
        return true # p is not inert in Q(zeta_{(m-1)/2})
    i = 0
    for ell in range(1,N*m^2*p^2,m*p):
        if not ZZ(ell).is_prime():
            continue
        Fellx.<x>=GF(ell)[]
        if ((x+1)^((ell-1)//p) -1 ) % cyclotomic_polynomial(m) != 0:
            return false # could prove that p doesn't divide h
        i += 1
        if i >= N:
            break
    return true # couldn't prove that p doesn't divide the class number