EC4

from random import randint

#Check if h is between 30 and 70
def CheckClassNumber(D):
    h = QuadraticField(D).class_number()
    if (h in Primes()) and h > 30 and h < 70:
        return h
    else:
        return 0

#Generate degrees of isogeny and j_invariant
def GetDegreesAndJ(D, K):
    j = hilbert_class_polynomial(D).roots(ring = K, multiplicities = False)
    if j == []:
        return 0, [], 0
    else:
        isogenyDegrees = list()
        for pr in Primes():
            if kronecker_symbol(D, pr) == 1:
                isogenyDegrees.append(pr)
            if len(isogenyDegrees) == 5:
                return j[0], isogenyDegrees, 1

#Elliptic Curve generation
def GenerateEC():
    while 1:
        p = next_prime(randint(10**4, 10**5))
        K = GF(p)
        A, B = randint(1, p), randint(1, p)
        p, A, B = 53831, 20363, 31029
        K = GF(p)
        E = EllipticCurve(K, [A, B])
        T = E.trace_of_frobenius()
        D = T**2 - 4*p
        h = CheckClassNumber(D)
        if h == 0:
            continue
        else:
            j, isogenyDegrees, flag = GetDegreesAndJ(D, K)
            if flag != 1:
                continue
            else:
                return E, isogenyDegrees, j
#
def GetCurveInvariant(E, j, degree):
    PHI = ClassicalModularPolynomialDatabase()
    K = GF(E.base_field().characteristic())
    phi = (PHI[degree]).change_ring(K)
    PR.<x> = PolynomialRing(K)
    curveInvariants = list()
    curveInvariants.append(j)
    buf = j
    while 1:
        root = phi(buf, x).roots(ring = K, multiplicities = False)
        if curveInvariants.count(root[0]) != 0:
            if curveInvariants.count(root[1]) != 0:
                break
            else:
                curveInvariants.append(root[1])
                buf = root[1]
        else:
            curveInvariants.append(root[0])
            buf = root[0]
    return curveInvariants


#Printing 5 graps
def PrintGraphs(E, isogenyDegrees, j):
    bufGraph = Graph()
    commonGraph = Graph()
    for degree in isogenyDegrees:
        isogenyCurveInvariants = GetCurveInvariant(E, j, degree)
        print("For degree {} invariants are next: {} " . format(degree, isogenyCurveInvariants))
        buf = isogenyCurveInvariants[len(isogenyCurveInvariants) - 1]
        for j_invariant in isogenyCurveInvariants:
            bufGraph.add_edge(buf, j_invariant)
            commonGraph.add_edge(buf, j_invariant)
            buf = j_invariant
        graph = bufGraph.graphplot(vertex_labels = True, vertex_size = 1200, layout = 'circular', vertex_color = 'white')
        graph.show(figsize=[12,12])
    print("Last one: ")
    graph = commonGraph.graphplot(vertex_labels = True, vertex_size = 1200, layout = 'circular', vertex_color = 'white')
    graph.show(figsize=[12,12])





#main
E, isogenyDegrees, j = GenerateEC()
print("E: {}" . format(E))
print("Degrees: {}" . format(isogenyDegrees))
print("j-invariant: {}" . format(j))
PrintGraphs(E, isogenyDegrees, j)