# Revision history [back]

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()
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
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
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:
XX = pointIsogeny( *XX )
count += 1
# print "Step %3s :: Passing to %s" % ( count, XX )
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.