An answer could be as follows:
'p=18446744073709551667'
then define the field
'K . < a >=GF(p^2)'
over which we define the curve
'E = EllipticCurve(K,[1,0])'. Now use the code
'E.cardinality().factor()' to get the cardinality of E, which is 24∗49992∗9225217080270832.
Use the code 'pts = E.gens()' to find the generators which turn out to be ${(3967812587828368975a + 7437411822947458617 : 16672546567636985011a + 14299206839437184441 : 1), (8681109523785822829a + 7072963597633280041 : 15515941078240688329a + 11174851755365315891 : 1)}$.
Let P1 be the first point, that is
'P1=E(13308581970510900443a + 15361864644717267429, 12035351063104026383a + 6630581890190451124)'
and use the code
'P1.order().factor()'
to get the order of P1, which is 22∗4999∗922521708027083. Now, let 'Q1=2^2 * 922521708027083 * P1', then Q1 is of the required order. We can repeat this with the other point.
Welcome to Ask Sage. Thank you for your question.
Ideally, provide a minimal concrete example, with code to define the objects, that can be copied and pasted in a fresh Sage session, leaving only the final step to answer.
p=18446744073709551667
E = EllipticCurve(GF(p^2),[1,0])
E is an Elliptic Curve defined by y^2 = x^3 + x over Finite Field in z2 of size 18446744073709551667^2 and I want to find generators P and Q for E[4999], but I can't find a suitable code.use
sage: p=E.torsion_polynomial(4999)
maybe ?