example 5.2.1 from Pairings for beginners

sage_5-2-1.sage

q = 5
Fq = GF(q)
R.<x> = Fq[]
Fq2.<i> = Fq.extension(x^2 + 2)
a = 0
b = -3

# Create elliptic curves over Fq and Fq2
E = EllipticCurve(Fq, [a, b])
E2 = EllipticCurve(Fq2, [a, b])

# debug
print(E.order().factor())
print(E2.order().factor())
print(E.is_supersingular())

# ?
r = 3 # torsion factor
h = E2.order() // r^2  # cofactor

# Coset representatives (using list comprehensions)
main_coset = set([r * P for P in E2.points()])
all_cosets = [main_coset]
points_seen = {}
for P in main_coset:
    points_seen[P] = True

# Generate cosets until all_cosets points are covered
print("looking for ", E2.order(), " cosets")
while len(all_cosets) < r * r: # r * r no?
    R = E2.random_point()
    point_seen = 0
    for P in main_coset:
        if P + R in points_seen:
            point_seen += 1
    # either we haven't seen all the points before
    # or we've seen all of them
    assert(point_seen == 0 or point_seen == len(main_coset))
            
    # add to all_cosets if not seen before
    if point_seen == 0:
        print("new point never seen before!")
        new_coset = set()
        for P in main_coset:
            new_point = P + R
            print("new point: ", new_point)
            new_coset.add(new_point)
            points_seen[new_point] = True
        all_cosets.append(new_coset)
        print("we now have ", len(all_cosets), " cosets")

# Torsion points (using list comprehensions)
TorsPts = [P for P in E2.points() if r * P == E2.points()[0]]
assert len(TorsPts) == r * r # isomorphic to Z_r * Z_R

# Flower generator function
def FlowerGenerator(TorsPts):
    V = []
    S = TorsPts.copy()
    petals = ZZ(len(S)).sqrt() + 1  # Convert to SageMath integer
    ptsInPet = petals - 1
    T = []

    for _ in range(petals):
        S = [P for P in S if P not in V]  # Set difference
        P = S[randint(0, len(S) - 1)]   # Choose random point 
        while P[2] == 0:                # Check for non-zero z-coordinate 
            P = S[randint(0, len(S) - 1)]
        V = [j * P for j in range(1, ptsInPet + 1)]
        T.append(V)
    return T

print(FlowerGenerator(TorsPts))