# Even degree large isogeny computation

I want to calculate $2^{372}$ and $3^{239}$ degree isogenies in Sage, but what I noticed is that using EllipticCurveIsogeny function does not finish executing for such large degree isogenies. There is already a partial answer by Dan Fulea here, but it does not seem to work for the even degree isogenies, and also for the above mentioned large degree odd isogeny as there is no chaining of small degree ones.

So, basically I have the following class from Dan's answer:

class EllipticCurveWithTorsionPoint(object):
"""
Class putting together in a bag:
- an elliptic curve,
- with a point R with same order as isogeny degree,
- and with a list of other points.
The main method <image> will associate an object of the same instance,
which is the "image" of the given structure w.r.t. the isogeny with kernel O, R, 2R, ... (d-1)R.
"""

def __init__(self, E, R, d=None, points=[]):
"""
E should be given in Weierstrass form.
"""
self.is_valid = True    # so far
self.E = E
if E != EllipticCurve(E.base_field(), [E.a4(), E.a6()]):
self.is_valid = False
self.points = points
self.R = R

if d:
self.d = ZZ(d)
else:
self.d = ZZ(R.order())
print "Computed: order of R is %s" % R.order().factor()

if 2.divides(self.d):
print "This class works only for an odd order d, but d=%s" % self.d
self.valid = False

self.kernel = [k*R for k in [0..(self.d-1)]]
self.d1 = ZZ((self.d-1) / 2)

def image(self):
"""
Pass to E', here denoted EE, the curve which is isogenous to E, so that
ker(E -> E') is {0, R, 2R, ... (d-1)R}. Also map to E' all the existing structure.
"""

if not self.is_valid:
print "Invalid object, no image is computed"
return self

E = self.E
F = E.base_field()
R = self.R

if R == E(0):
print "No further isogeny can be built"
return self

a, b = E.a4(), E.a6()
v, w = F(0), F(0)
data = []

for mult in [1..self.d1]:
# We split the set G = R, 2R, ... , (d-2)R, (d-1)R in two, G(+)
# and G(-) so that P in G(+) <=> -P in G(-), the multiplicities
# above correspond to those for G(+).
P        = mult*R # point in G(+)
xP , yP  = P.xy()
gxP, gyP = 3*xP^2 + a, -2*yp
vP , uP  = 2*gxP     , gyP^2

v += vP
w += uP + xP*vP

data.append([xP, yP, uP, vP])

EE = EllipticCurve(F, [a-5*v, b-7*w])

def phi(point):
print "Computing phi(%s)..." % str(point)
if point in self.kernel:
return EE(0)

x , y  = point.xy()
xx, yy = x, y

for xP, yP, uP, vP in data:
xx += vP    /(x-xP)   + uP  /(x-xP)^2
yy -= uP*2*y/(x-xP)^3 + vP*y/(x-xP)^2

return EE((xx, yy))

return EllipticCurveWithTorsionPoint(EE
, EE(0)
, d = 1
, points = [phi(point) for point in self.points])

def imageStrandard(self):
"""
This does the same thing as "image" method but using the Sage
isogeny implementation.
"""

if not self.is_valid:
return self

E = self.E
R = self.R

phi = EllipticCurveIsogeny(E, R)
EE  = phi.codomain()

return EllipticCurveWithTorsionPoint(EE
, EE(0)
, d = 1
, points = [phi(point) for point in self.points])


But when | try to calculate isogeny, using the following field and elliptic curve:

p = 10354717741769305252977768237866805321427389645549071170116189679054678940682478846502882896561066713624553211618840202385203911976522554393044160468771151816976706840078913334358399730952774926980235086850991501872665651576831
Fp = GF(p)
R.<x> = PolynomialRing(Fp)
# The quadratic extension via x^2 + 1 since p = 3 mod 4
Fp2.<j> = Fp.extension(x^2 + 1)
E = EllipticCurve(Fp2, [1,0])


and the following elements:

R = (9007518108646169717637829256143902727256908604612852170262845383236308734546752870948023665818612994373405439104311563180515971827416888758364379345147971116263603311381076594736408857657724917306603115510356333363208849059629*j + 8143875544876203102731118758998042519548033956324586056599014675159430178641351639084698635604334996400279245810044652728374901773305503117205094200107841651156165349992200320758569566012680521517849444975619314122642286738078 : 5933067621852699133314119054291511797259450704514751366342623924502189539964363940282453138760679452407770908256363554867263728524097776132882938143129922521585626385240622607283870971683009886348379499590584807528366215593257*j + 3243060684426863808401390460086135176972427334603406644358264254553792355536601776050361287848102972929373427788393162068323805718475829130551068182582167101080584984966674872491237953303465307684436948645319932524248022850554 : 1)
phiP = (5784307033157574162391672474522522983832304511218905707704962058799572462719474192769980361922537187309960524475241186527300549088533941865412874661143122262830946833377212881592965099601886901183961091839303261748866970694633 : 5528941793184617364511452300962695084942165460078897881580666552736555418273496645894674314774001072353816966764689493098122556662755842001969781687909521301233517912821073526079191975713749455487083964491867894271185073160661 : 1)
phiQ = (4570410708611731090586095763344282337595085134330165462411227620255106477963004653732902534638529526314592687143599015857903362887988612527631285807628029554145760006701700452765434631350888025796273995011688240123798680882198 : 5528941793184617364511452300962695084942165460078897881580666552736555418273496645894674314774001072353816966764689493098122556662755842001969781687909521301233517912821073526079191975713749455487083964491867894271185073160661*j : 1)


such that

X = EllipticCurveWithTorsionPoint(E, R, None, [phiP, phiQ])
W = X.image()


It doesn't seem to finish execution, where R has order $3^{239}$. And the class doesn't even work for an even degree isogeny, for example, evaluation this one:

R = (5110605836845566926025654047424167572170042372900551706614198682868609124986515618280597159249062586518558366705653348177099645521579421315786637369850653444795496916683519951328090963213358196198661324642758054901662568911454*j + 4234869830555943849779568791108578084407579904872225169650118735805222570371281152259184280337943783338310393473931143980765745444039692659441232529288102314058310400426732946525612162308610371741399922257500660578782200917155 : 5100498610209823809391228850023209691234397421134363575826918390640143899484912001953261335357385806761065415074391051805515401642580444835961489218757392256505375361061973605756427459420249731757621385279542696786456958672535*j + 3945657958394557421465268691904614939064430526791563586067697568085750995914128796629365705178878140816178998799696237175564589123664889202169294820353858542756862036043886052635535912561145255328107702126311995666291927980188 : 1)
phiP = (4359917396849101231053336763700300892915096700013704210194781457801412731643988367389870886884336453245156775454336249199185654250159051929975600857047173121187832546031604804277991148436536445770452624367894371450077315674371 : 106866937607440797536385002617766720826944674650271400721039514250889186719923133049487966730514682296643039694531052672873754128006844434636819566554364257913332237123293860767683395958817983684370065598726191088239028762772 : 1)
phiQ = (5994800344920204021924431474166504428512292945535366959921408221253266209038490479113012009676730260379396436164503953186018257726363502463068559611723978695788874294047308530080408582516238481209782462483097130422588335902460 : 106866937607440797536385002617766720826944674650271400721039514250889186719923133049487966730514682296643039694531052672873754128006844434636819566554364257913332237123293860767683395958817983684370065598726191088239028762772*j : 1)


such that

X = EllipticCurveWithTorsionPoint(E, R, None, [phiP, phiQ])
W = X.image()


where R has degree $2^{372}$. How can one chain 2 and 3 isogenies compute high degree isgonies that are powers of 2 and 3?. Any ideas how to do these efficiently in Sage so that the computations finish a "reasonable" time?

edit retag close merge delete

Sort by » oldest newest most voted

Recall the posted situation, slightly changed, but the same one:

p  = 2^372 * 3^239 - 1
Fp = GF(p, proof=False)
R.<x> = PolynomialRing(Fp)
# The quadratic extension via x^2 + 1 since p = 3 mod 4
Fp2.<j> = Fp.extension(x^2 + 1)

E = EllipticCurve( Fp2, [1,0] )

R = E.point( ( 9007518108646169717637829256143902727256908604612852170262845383236308734546752870948023665818612994373405439104311563180515971827416888758364379345147971116263603311381076594736408857657724917306603115510356333363208849059629*j + 8143875544876203102731118758998042519548033956324586056599014675159430178641351639084698635604334996400279245810044652728374901773305503117205094200107841651156165349992200320758569566012680521517849444975619314122642286738078, 5933067621852699133314119054291511797259450704514751366342623924502189539964363940282453138760679452407770908256363554867263728524097776132882938143129922521585626385240622607283870971683009886348379499590584807528366215593257*j + 3243060684426863808401390460086135176972427334603406644358264254553792355536601776050361287848102972929373427788393162068323805718475829130551068182582167101080584984966674872491237953303465307684436948645319932524248022850554 ) )

phiP = E.point( ( 5784307033157574162391672474522522983832304511218905707704962058799572462719474192769980361922537187309960524475241186527300549088533941865412874661143122262830946833377212881592965099601886901183961091839303261748866970694633, 5528941793184617364511452300962695084942165460078897881580666552736555418273496645894674314774001072353816966764689493098122556662755842001969781687909521301233517912821073526079191975713749455487083964491867894271185073160661 ) )

phiQ = E.point( ( 4570410708611731090586095763344282337595085134330165462411227620255106477963004653732902534638529526314592687143599015857903362887988612527631285807628029554145760006701700452765434631350888025796273995011688240123798680882198, 5528941793184617364511452300962695084942165460078897881580666552736555418273496645894674314774001072353816966764689493098122556662755842001969781687909521301233517912821073526079191975713749455487083964491867894271185073160661*j ) )


Here is a solution working for both even and odd order of the kernel. Of course, if the kernel is a "big" composite number, it is a good idea to split it into factors, and use the composition isogeny, this telescopical way to proceed is better.

If i understand the situation and the task, things are as follows:

We initialize an elliptic curve $E=E_0$ over a big field $F$, and declare some rational points, first of all $R=R_0$ is a point of order $N=3^{239}$ defining a subgroup of rational points, in our case

sage: R.order().factor()
3^239


and then we pick two further points, phiP and phiQ, on the curve.

Let $R_1=3R$, $R_2=9R$, and so on, $R_n=3^nR$. In our case, $R_N=R_{239}=0$, $N=239$, and let $R'$ be the last non-zero point in the chain $R_0, R_1,R_2,\dots\$ So $$R'=R_{N-1}\ .$$

We consider the isogeny $\phi$ from $E=E_0$ to its codomain $E'=E_1$, in code EE, which kills $R'$ (and its multiples).

In a picture: \begin{align} E_0\ &\quad R_0\overset3\to R_1\overset3\to R_2\overset3\to R_3\overset3\to \dots\overset3\to R_{N-1}\overset3\to \underbrace{R_N}_0 \\ \phi\downarrow\ & \\ E_1\ &\quad S_0\overset3\to S_1\overset3\to S_2\overset3\to S_3\overset3\to \dots\overset3\to \underbrace{S_{N-1}}_0\overset3\to \underbrace{S_N}_0 \end{align} Here, $S_0=\phi(R)$ is the image of $R_0$ through $\phi$, and its order is multiplicatively decreased by the factor $3$ compartively to $R$.

We repeat this process, till $S$ is zero.

This time the code will not use a class, since it seems to be not so easy to make transparent the computational idea. Instead, we use plain direct functional programming.

We will insert in a list the points $R_0, R_1, \dots, R_{N-1}$, then construct pointwise isogenies w.r.t. the elements in the list, as long as the list is not empty. All points in the list are mapped through the isogeny into the new curve.

Time for the code.

def pointIsogeny( E, R_list, points, degrees_list=[] ):
"""The last element, R, in R_list must have order
- case EVEN, either two
- case ODD,  or some odd number (prime or not)
The degrees list has or does not have the parallel information about the orders
of the points in the R_list. If no information is available, the code will
try to compute it.

We then construct the corresponding isogeny phi *at point level*.
E must be an elliptic curve in short Weierstrass form.
E  :: y^2 = x^3 + a*x + b

CASE EVEN: R = (s,0) has order two, and s is defined in this way.
CASE ODD:  R has two components, we will build the group generated by R.

We construct ad-hoc the isogeny map phi at point level following section 3 in
https://arxiv.org/pdf/math/0404124.pdf

This gives a wonderful unifying description, that i missed in other sources.
(It is best suited for a computer implementation (for the case of a cyclic kernel isogeny).)

We get a map (x,y) -> phi(x,y) = X, Y.
The article gives an unifying way to look at X-component for both even and odd torsion cases.
The mention of the Y-component in loc. cit., simply obtained by differentiating
the formula for X w.r.t. x.
"""

if not R_list:
return ( E, R_list, points )

R = R_list.pop()    # R is extracted from the last place of R_list

if degrees_list:
d = degrees_list.pop()
else:
d = R.order()    # have to compute it...
kernel = [ k*R for k in [ 0 .. (d-1)] ]

if d*R != E(0) or len( [ kR for kR in kernel if kR == E(0) ] ) > 1:
print "ERROR, no d-torsion point."    # we insist to have order exactly d

# case EVEN, R is a 2-torsion point
if d == 2:

a, b = E.a4(), E.a6()
s = R.xy()[0]
t = 3*s^2 + a
w = s*t

EE = EllipticCurve( E.base_field(), [  a-5*t, b-7*w ] )

def phi(P):
"""P is a point on E, we return a point on EE
"""
if P in ( E(0), R ):
return EE(0)

x, y = P.xy()
X = x + t/(x-s)
Y = y * ( 1 - t/(x-s)^2 )

return EE.point( (X, Y) )

elif not 2.divides(d):

F = E.base_field()
a4, a6 = E.a4(), E.a6()
b2, b4, b6 = E.b2(), E.b4(), E.b6()
t, u, w = F(0), F(0), F(0)
data = []

d1 = ZZ( (d-1)/2 )
for P in kernel[ 1 : d1+1 ]:
# We split the set G = R, 2R, ... , (d-2)R, (d-1)R in two, G(+)
# and G(-) so that P in G(+) <=> -P in G(-), the multiplicities
# above correspond to those for G(+).
xP , yP = P.xy()
tP = 6*xP^2 + b2*xP + b4
uP = 4*xP^3 + b2*xP^2 + 2*b4*xP + b6
wP = uP + xP * tP

t += tP
u += uP
w += wP

data.append( [xP, yP, tP, uP] )

EE = EllipticCurve(F, [a4-5*t, a6-b2*t-7*w])

def phi(point):
if point in kernel:
return EE(0)

x , y = point.xy()
X , Y = x, y

for xP, yP, tP, uP in data:
X +=        tP / (x-xP)   + uP   / (x-xP)^2
Y += y * ( -tP / (x-xP)^2 - uP*2 / (x-xP)^3 )    # y*diff(xx,x)

return EE( (X, Y) )

else:
# torsion point of order 4,6,8,...
print "*** ERROR ***"
return

return ( EE
, [ phi(P) for P in R_list ]
, [ phi(P) for P in points ]
, degrees_list )

# Let us test the above function...
#
E = EllipticCurve( GF(19), [-1,0] )
R_list = [ E.point( (1,0) ) ]
points = E.rational_points()

# use the sage native isogeny:
R = R_list[0]
phi = EllipticCurveIsogeny( E, R )
# for P in points:
#     print "%12s -> %12s" % ( P, phi(P) )

# use degreeTwoPointIsogeny
X = ( E, R_list, points )
XX = pointIsogeny( *X )
EE, S_list, newpoints, degrees_list = XX

print "Codomain elliptic curve is (ad-hoc computation)\n%s" % EE
print "Codomain elliptic curve is (sage computation)\n%s" % phi.codomain()

print "Still usable isogeny points: %s" % S_list
print "Isogeny is mapping the points as follows:"
for k in range(len(points)):
P, PP = points[k], newpoints[k]
print "%-13s -> %-13s (ad-hoc) versus %-13s (sage)" % ( P, PP, phi(P) )


This gives so far:

Codomain elliptic curve is (ad-hoc computation)
Elliptic Curve defined by y^2 = x^3 + 8*x + 5 over Finite Field of size 19
Codomain elliptic curve is (sage computation)
Elliptic Curve defined by y^2 = x^3 + 8*x + 5 over Finite Field of size 19
Still usable isogeny points: []
Isogeny is mapping the points as follows:
(0 : 0 : 1)   -> (17 : 0 : 1)  (ad-hoc) versus (17 : 0 : 1)  (sage)
(0 : 1 : 0)   -> (0 : 1 : 0)   (ad-hoc) versus (0 : 1 : 0)   (sage)
(1 : 0 : 1)   -> (0 : 1 : 0)   (ad-hoc) versus (0 : 1 : 0)   (sage)
(2 : 5 : 1)   -> (4 : 14 : 1)  (ad-hoc) versus (4 : 14 : 1)  (sage)
(2 : 14 : 1)  -> (4 : 5 : 1)   (ad-hoc) versus (4 : 5 : 1)   (sage)
(3 : 9 : 1)   -> (4 : 14 : 1)  (ad-hoc) versus (4 : 14 : 1)  (sage)
(3 : 10 : 1)  -> (4 : 5 : 1)   (ad-hoc) versus (4 : 5 : 1)   (sage)
(5 : 5 : 1)   -> (15 : 2 : 1)  (ad-hoc) versus (15 : 2 : 1)  (sage)
(5 : 14 : 1)  -> (15 : 17 : 1) (ad-hoc) versus (15 : 17 : 1) (sage)
(6 : 1 : 1)   -> (14 : 7 : 1)  (ad-hoc) versus (14 : 7 : 1)  (sage)
(6 : 18 : 1)  -> (14 : 12 : 1) (ad-hoc) versus (14 : 12 : 1) (sage)
(9 : 6 : 1)   -> (14 : 7 : 1)  (ad-hoc) versus (14 : 7 : 1)  (sage)
(9 : 13 : 1)  -> (14 : 12 : 1) (ad-hoc) versus (14 : 12 : 1) (sage)
(11 : 3 : 1)  -> (15 : 17 : 1) (ad-hoc) versus (15 : 17 : 1) (sage)
(11 : 16 : 1) -> (15 : 2 : 1)  (ad-hoc) versus (15 : 2 : 1)  (sage)
(12 : 5 : 1)  -> (7 : 9 : 1)   (ad-hoc) versus (7 : 9 : 1)   (sage)
(12 : 14 : 1) -> (7 : 10 : 1)  (ad-hoc) versus (7 : 10 : 1)  (sage)
(15 : 4 : 1)  -> (7 : 9 : 1)   (ad-hoc) versus (7 : 9 : 1)   (sage)
(15 : 15 : 1) -> (7 : 10 : 1)  (ad-hoc) versus (7 : 10 : 1)  (sage)
(18 : 0 : 1)  -> (17 : 0 : 1)  (ad-hoc) versus (17 : 0 : 1)  (sage)


so we have a good check for the routine above. The formulae were taken for the first component, X in the above function, from section 3 of

(Everest, King, arXiv, 0404124.pdf)[https://arxiv.org/pdf/math/0404124.pdf]

which give an excellent unifying point of view for both even and odd cases.

"Small" digression. The above article gives only the formula for the $x$--component, to get the one for the $y$--component, one can proceed as follows. (The following applies only for the two-torsion case.) Let us start with the elliptic curve $E$ given by $y^2=x^3 +ax+b$, which should have a $2$-torsion point $(s,0)$ (over the tacitly chosen definition field). Then we get rid of $b$, using $s$ instead, $-b=s^3+as$. We also know the equation for the codomain of the to-be-constructed isogeny $\phi$, all the data put together: \begin{align} E\ :\ y^2 &= x^3+ax+b\ ,\\ y^2 &= (x^3-s^3) + a(x-s)\ ,\\ t &= 3s^2 +a\ ,\\ w &= st\ ,\\ E'\ :\ Y^2 &= X^3 + (a-5t)X +(b-7w)\ ,\\ (x,y)&\to (X,Y)=\phi(x,y) \ ,\\ X &= x+\frac t{x-s}\ ,\\ Y &=?\ . \end{align} To obtain the value for $Y$, for the second (affine) component in the definition of $\phi$, we have the Ansatz that $Y$ should be of the shape $Y = yP/(x-s)^2$, where $P$ is a rational function of $x$. So we only have to detect the $P$. The following sage code does the job:

var( 'x,a,s' );

b  = -s^3 -a*s
yy = x^3 + a*x + b
t  = 3*s^2 + a
w  = s*t
X  = x + t/(x-s)
YY = X^3 + (a-5*t)*X + (b-7*w)
# we search an expression of the shape Y = yP/(x-s)^2, i.e. YY = yy PP / (x-s)^4
# so let us determine P by factorizing...
factor( YY / yy * (x-s)^4 )


This gives:

(2*s^2 + 2*s*x - x^2 + a)^2


From here, after considerying rather -(2*s^2 + 2*s*x - x^2 + a) as the "true" square root of the above, after writing it in the form $(x^2-2sx+s^2)-(3s^2+a)=(x-s)^2-t$ we obtain the beautiful formula for the two-Torsion isogeny: \begin{align} (X,Y)&= \phi(x,y)\ ,\\ X &= x + \frac t{x-s}\ ,\\ Y &= y\left(1-\frac t{(x-s)^2}\right)=y\frac{\partial X}{\partial x}\ .\\ \end{align}

Note that taking the derivative of X to get Y also works for the odd-degree-torsion case! (This was already somehow in my blood before, but i could not find any explicit reference, maybe i have to write a one page article exactly on this small point, it is too important structurally, mnemotechnically, and for an implementation as here. Same applies in the odd degree case.)

So let us also test the odd degree torsion case. Same curve, an other $R$.

# test for a 5-torsion point R
E = EllipticCurve( GF(19), [-1,0] )
R_list = [ E.point( (5,5) ) ]    # a point of odd order 5
points = E.rational_points()

# use the sage native isogeny:
R = R_list[0]
phi = EllipticCurveIsogeny( E, R )

# use pointIsogeny
X = ( E, R_list, points )
XX = pointIsogeny( *X )
EE, S_list, newpoints, degrees_list = XX
print "Codomain elliptic curve is (ad-hoc computation)\n%s" % EE
print "Codomain elliptic curve is (sage computation)\n%s" % phi.codomain()

print "Still usable isogeny points: %s" % S_list
print "Isogeny is mapping the points as follows:"
for k in range(len(points)):
P, PP = points[k], newpoints[k]
print "%-13s -> %-13s (ad-hoc) versus %-13s (sage)" % ( P, PP, phi(P) )


This gives a perfect match:

Codomain elliptic curve is (ad-hoc computation)
Elliptic Curve defined by y^2 = x^3 + 13*x over Finite Field of size 19
Codomain elliptic curve is (sage computation)
Elliptic Curve defined by y^2 = x^3 + 13*x over Finite Field of size 19
Still usable isogeny points: []
Isogeny is mapping the points as follows:
(0 : 0 : 1)   -> (14 : 0 : 1)  (ad-hoc) versus (14 : 0 : 1)  (sage)
(0 : 1 : 0)   -> (0 : 1 : 0)   (ad-hoc) versus (0 : 1 : 0)   (sage)
(1 : 0 : 1)   -> (0 : 0 : 1)   (ad-hoc) versus (0 : 0 : 1)   (sage)
(2 : 5 : 1)   -> (5 : 0 : 1)   (ad-hoc) versus (5 : 0 : 1)   (sage)
(2 : 14 : 1)  -> (5 : 0 : 1)   (ad-hoc) versus (5 : 0 : 1)   (sage)
(3 : 9 : 1)   -> (14 : 0 : 1)  (ad-hoc) versus (14 : 0 : 1)  (sage)
(3 : 10 : 1)  -> (14 : 0 : 1)  (ad-hoc) versus (14 : 0 : 1)  (sage)
(5 : 5 : 1)   -> (0 : 1 : 0)   (ad-hoc) versus (0 : 1 : 0)   (sage)
(5 : 14 : 1)  -> (0 : 1 : 0)   (ad-hoc) versus (0 : 1 : 0)   (sage)
(6 : 1 : 1)   -> (0 : 1 : 0)   (ad-hoc) versus (0 : 1 : 0)   (sage)
(6 : 18 : 1)  -> (0 : 1 : 0)   (ad-hoc) versus (0 : 1 : 0)   (sage)
(9 : 6 : 1)   -> (0 : 0 : 1)   (ad-hoc) versus (0 : 0 : 1)   (sage)
(9 : 13 : 1)  -> (0 : 0 : 1)   (ad-hoc) versus (0 : 0 : 1)   (sage)
(11 : 3 : 1)  -> (0 : 0 : 1)   (ad-hoc) versus (0 : 0 : 1)   (sage)
(11 : 16 : 1) -> (0 : 0 : 1)   (ad-hoc) versus (0 : 0 : 1)   (sage)
(12 : 5 : 1)  -> (5 : 0 : 1)   (ad-hoc) versus (5 : 0 : 1)   (sage)
(12 : 14 : 1) -> (5 : 0 : 1)   (ad-hoc) versus (5 : 0 : 1)   (sage)
(15 : 4 : 1)  -> (14 : 0 : 1)  (ad-hoc) versus (14 : 0 : 1)  (sage)
(15 : 15 : 1) -> (14 : 0 : 1)  (ad-hoc) versus (14 : 0 : 1)  (sage)
(18 : 0 : 1)  -> (5 : 0 : 1)   (ad-hoc) versus (5 : 0 : 1)   (sage)


So we have now a good framework to attack the situation in the post. Let us do this. I will copy the posted data again.

p  = 2^372 * 3^239 - 1
Fp = GF(p, proof=False)
R.<x> = PolynomialRing(Fp)
# The quadratic extension via x^2 + 1 since p = 3 mod 4
Fp2.<j> = Fp.extension(x^2 + 1)

E = EllipticCurve( Fp2, [1,0] )

R = E.point( ( 9007518108646169717637829256143902727256908604612852170262845383236308734546752870948023665818612994373405439104311563180515971827416888758364379345147971116263603311381076594736408857657724917306603115510356333363208849059629*j + 8143875544876203102731118758998042519548033956324586056599014675159430178641351639084698635604334996400279245810044652728374901773305503117205094200107841651156165349992200320758569566012680521517849444975619314122642286738078, 5933067621852699133314119054291511797259450704514751366342623924502189539964363940282453138760679452407770908256363554867263728524097776132882938143129922521585626385240622607283870971683009886348379499590584807528366215593257*j + 3243060684426863808401390460086135176972427334603406644358264254553792355536601776050361287848102972929373427788393162068323805718475829130551068182582167101080584984966674872491237953303465307684436948645319932524248022850554 ) )

phiP = E.point( ( 5784307033157574162391672474522522983832304511218905707704962058799572462719474192769980361922537187309960524475241186527300549088533941865412874661143122262830946833377212881592965099601886901183961091839303261748866970694633, 5528941793184617364511452300962695084942165460078897881580666552736555418273496645894674314774001072353816966764689493098122556662755842001969781687909521301233517912821073526079191975713749455487083964491867894271185073160661 ) )

phiQ = E.point( ( 4570410708611731090586095763344282337595085134330165462411227620255106477963004653732902534638529526314592687143599015857903362887988612527631285807628029554145760006701700452765434631350888025796273995011688240123798680882198, 5528941793184617364511452300962695084942165460078897881580666552736555418273496645894674314774001072353816966764689493098122556662755842001969781687909521301233517912821073526079191975713749455487083964491867894271185073160661*j ) )

points = [ phiP, phiQ ]

R_list = []
degrees_list = []
S = R
while S != E(0):
R_list      .append( S )
degrees_list.append( 3 )
S = 3*S
X = ( E, R_list, points, degrees_list )
XX = X

count = 0
while XX[1]:
XX = pointIsogeny( *XX )
count += 1
# print "Step %3s :: Passing to %s" % ( count, XX[0] )
print "Step %3s DONE" % count


In this moment, the object XX contains the image curve E239, and the image points points239 as follows:

E239, _, points239, _ = XX


It remains to print and use the two objects E239 and points239. I think, the question wanted them.

more

@dan_fulea Thank you very much for the beautiful reply. The paper also looks quite interesting.

( 2018-02-28 17:01:10 -0500 )edit