polycyclic presentation

edit

Actually I found it much more efficient to use a direct computation in the class group. All my attempts to interact with GAP were quite slow.

def get_bqf(D,p):
    d = fundamental_discriminant(D)
    v = isqrt(D//d)
    if v % p == 0:
        return None
    ks = kronecker_symbol(D,p)
    if ks == -1:
        return None
    c= floor(-D/(4*p))
    while True:
        c += 1
        b2 = D+4*p*c
        if isqrt(b2)**2 == b2:
            b = isqrt(b2)
            return  BinaryQF([p,b,c])
    return None

def get_prime_forms(D):
    P = Primes()
    p = 2
    B = isqrt(abs(D)/3)
    gens = []
    while p < B:
        bq = get_bqf(D,p)
        if  bq is not None :
            gens.append(bq)
        p = P.next(p)
    return gens

# p.16 of https://arxiv.org/pdf/0903.2785

def TableInsert(T,g):
    T.append(g.reduced_form())
    return None

def TableLookup(T,beta):
    b = beta.reduced_form()
    if b in T:
        return T.index(b)
    else:
        return None

def poly_rep(gens):
    G = (gens[0].form_class()).parent()
    n = len(gens)
    T = []
    r = []
    s = []
    TableInsert(T,G.zero().form())
    for i in range(1,n+1):
        beta = gens[i-1].reduced_form()
        ri = 1
        N = len(T)
        si = TableLookup(T,beta)
        while si == None:
            for j in range(N):
                TableInsert(T,beta*T[j])
            beta = beta*gens[i-1]
            si = TableLookup(T,beta)
            ri += 1
        r.append(ri)
        s.append(si)
    return T,r,s

D = -5291
gens = get_prime_forms(D)  
T,r,s=poly_rep(gens)

This does not include improvements from p.17 of Computing Hilbert Class Polynomials With Chinese Remainder Theorem from Andrew Sutherland.