# Isogeny computation does not finish in Sage

I'm using the EllipticCurveIsogeny function to calculate isogenies, but what I have noticed is that the function does not finish executing (at least in a reasonable amount of time). My code segment looks like this:

proof.arithmetic(False)

params = [
5784307033157574162391672474522522983832304511218905707704962058799572462719474192769980361922537187309960524475241186527300549088533941865412874661143122262830946833377212881592965099601886901183961091839303261748866970694633,
5528941793184617364511452300962695084942165460078897881580666552736555418273496645894674314774001072353816966764689493098122556662755842001969781687909521301233517912821073526079191975713749455487083964491867894271185073160661,
4359917396849101231053336763700300892915096700013704210194781457801412731643988367389870886884336453245156775454336249199185654250159051929975600857047173121187832546031604804277991148436536445770452624367894371450077315674371,
106866937607440797536385002617766720826944674650271400721039514250889186719923133049487966730514682296643039694531052672873754128006844434636819566554364257913332237123293860767683395958817983684370065598726191088239028762772
]

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)

E0 = EllipticCurve(Fp2, [1,0])
assert E0.is_supersingular()

E = EllipticCurve(Fp2, [1,0]) # Weierstrass curve y^2=x^3+x
RR.<x> = PolynomialRing(Fp2)   # for computing isogeny kernels
phiP = E0([params[0], params[1]])
phiQ = E0([-params[0], j*params[1]])

SK = 210
eB = 239

P = E0([params[2], params[3]])
Q = E0([-params[2], j*params[3]])
R = P + SK * Q

for e in range(eB-1, 0, -1):
S = 3^e * R
ker = [x-P[0] for P in [z*S for z in [1..3]]]
ker = reduce(lambda x, y: x*y, ker)
print "ker ({}) = {}".format(e, ker)
ker_roots = ker.roots(ring=Fp2, multiplicities=False)
print "ker_roots ({}) = {}".format(e, ker_roots)
points = [E.lift_x(a) for a in ker_roots]
print "points ({}) = {}".format(e, points)
phi = EllipticCurveIsogeny(E, points)
print "phi ({}) = {}".format(e, phi)
E = phi.codomain()
R = phi(R)
phiP = phi(phiP)
phiQ = phi(phiQ)

coeffs = E.coefficients()
print "coeffs =", coeffs


I put multiple print statements inside the loop to see where the bottleneck is, and I noticed that the first iteration executes successfully, and everything is calculated, but then in the second iteration of the loop, when it comes to calculating the isogeny, it somehow doesn't finish. I actually run the loop without the isogeny computation, and it manages to find all the roots and lifts in a short period of time, so the bottleneck here is the isogeny calculations from the given points. What's causing Sage not to be able to finish the isogeny calculations, and how can I correct it so that my loop executes successfully?

edit retag close merge delete

Sort by » oldest newest most voted

We will implement the $3$--isogenies in an ad-hoc manner, using the following formulas, (Velu, special case of a three torsion point generating the kernel).

Let $E$ be an elliptic curve with short Weierstrass equation $$E\ :\qquad y^2 = x^3 +ax+b\ .$$ Let $$g(x,y) = x^3 +ax+b-y^2$$ be explicitly the function, applied on an affine point $(x,y)$, so that $E$ is given by $g=0$. We will denote by $g'_x$, $g'_y$, the derivatives of $g$ with respect to the variables. For a general point $R(x_R,y_R)$ on the curve, we write simpler $g^x_R$ the value for $g'_x(R)$. (Only one index can win and remain a lower index.)

Let now $R$ be a point of order three on $E$, i.e. $R\ne O$, $3R=O$. Its components will be denoted with $$(x_R, y_R)\ .$$ Then the group spanned by $R$ has elements $O,R,-R$. We define then the following quantities and objects: \begin{align} g_R^x &= g'_x(R) = (x^3 +ax+b)'_x(R)= 3x_R^2 +a\ ,\\ g_R^y &= g'_y(R) = (-y^2)'_y(R) = -2y_R\ ,\\ v_R &= 2g_R^x\ ,\\ u_R &= (g_R^y)^2\ ,\\ v &= v_R\ ,\\ w &= u_R + x_Rv_R\ .\\ E&:\ y^2 =x^3 +ax+b\ ,\\ E'&:\ y^2 =x^3 +a'x+b'\ ,\\\ &\qquad\text{ with } a'=a-5v\ ,\ b'=b-7w\ ,\\ \phi&:\ E\to E'\ ,\\ \phi(x,y)&= \left( \ x+\frac{v_R}{x-x_R}-\frac{u_R}{(x-x_R)^2} \ ,\ y +\frac{v_R y}{(x-x_R)^2} - \frac{2u_R y}{(x-x_R)^3}\ \right)\ ,\\ \end{align} The notations have been borrowed from the following article: Analogues of Vélu's formulas for..., Dustin Moody, Daniel Shumow, pages 3-4

However, the formula for $\phi$ found there did not work, had to calculate with bare hands and find the right one. (The formula for the target isogenous curve $E'$ was ok, it was a big present to have the numbers $a'=a-5v$ and $b'=b-7w$, so i stayed with the article.)

(Please also check the formulas. This was the main reason for implementing also a comparison.)

The code now. We will design a self made class EllipticCurveWithAThreeTorsionPoint with two methods. offering a chance to follow both roads:

• go with the computations in the posted code, which uses EllipticCurveIsogeny in the method imageAsPosted, (and constructs the isogeny algebraically, but we need less...) Note that it is no reason to pass from $R$ to the kernel polynomial, than back to $R$ in the case of interest, where $R$ is a three-torsion point.
• and go the ad-hoc direction, using the explicit formulas from loc. cit.. This is definitively simpler, and does exactly what is needed, namely recursively get the target isogenous curve, and map some known points from the source curve to the target curve. (We do not need the isogeny map $\phi$ as a map between curves, we only need to know how to map points. The above formula is then easily implemented.) The ad-hoc method is image.

First. we will compare the results in some special cases, where "smaller" primes are involved. If this works fine, and we get more confidence by using the ad-hoc solution, then we pass to the case of interest, which uses

• the posted "big" prime $$p = 2^{372} \cdot 3^{239} -1\ ,$$
• and the posted curve $E$ defined over $\mathbb F_p$ of order $|\;E(\mathbb F_p)\; |=(p+1)$,
• and a point $P$ defined over $\mathbb F_p$ of order corresponding to the $3$-power inside $(p+1)$.

First of all, let us search for some "small" primes $p$, so that we can construct a similar situation.

myprimes = []
for k3 in range( 10, 20 ):
p3 = 3^k3
for k2 in range( k3+1, 30)  :
N = ZZ(p3 * 2^k2)
if (N-1).is_prime():
myprime = N-1
E = EllipticCurve( GF(myprime), [1,0] )
if E.order() == N:
print N.factor()
myprimes.append( myprime )
myprimes . sort()
for myprime in myprimes[:10]:
print "%s -1 = %s " % ( (myprime+1).factor(), myprime )


This gives:

2^13 * 3^12 -1 = 4353564671
2^20 * 3^11 -1 = 185752092671
2^18 * 3^13 -1 = 417942208511
2^23 * 3^10 -1 = 495338913791
2^17 * 3^16 -1 = 5642219814911
2^22 * 3^15 -1 = 60183678025727
2^27 * 3^13 -1 = 213986410758143
2^29 * 3^12 -1 = 285315214344191
2^23 * 3^16 -1 = 361102068154367


and we pick in the sequel one case, the last one. The class now:

class EllipticCurveWithAThreeTorsionPoint( object ):
"""Class putting together in a bag:
- an elliptic curve,
- with a point R of order three,
- which comes from a point of order a power of three,
- 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.
"""

def __init__( self, E, points, threePowerTorsionPoint, power=None ):
"""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.threePowerTorsionPoint = threePowerTorsionPoint

if power:
self.power = power
else:
self.power = threePowerTorsionPoint.order()

self.power = ZZ(self.power)
if self.power.prime_factors() not in ( [], [3] ):
self.is_valid = False

if self.power == 1:
self.exponent = 0
else:
# quick+dirty, get the exponent of the first (and only) (maximal) prime power factor
self.exponent = self.power.factor()[0][1]

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 } .
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()

if self.exponent == 0:
print "...no further isogeny can be built"
return self

R = ZZ(self.power/3) * self.threePowerTorsionPoint

a, b     = E.a4(), E.a6()
xR, yR   = R.xy()
gxR, gyR = 3*xR^2 + a, -2*yR
vR, uR   = 2*gxR, gyR^2
v, w     = vR, uR + xR*vR

EE = EllipticCurve( F, [a-5*v, b-7*w] )
def phi(point):
if point in ( E(0), R, -R ):
return EE(0)
# else
x,y = point.xy()
return EE( ( x + vR/(x-xR) + uR/(x-xR)^2 ,
y - uR*2*y/(x-xR)^3 - vR*y/(x-xR)^2 ) )

return EllipticCurveWithAThreeTorsionPoint( EE
, [ phi(point) for point in self.points ]
, phi( self.threePowerTorsionPoint )
, ZZ( self.power/3 ) )

def imageAsPosted( self ):
"""This does exactly what the post does,
the step of computing the kernel polynomial is omitted,
since we already have the kernel points.
"""
if not self.is_valid:
return self
E = self.E
R = ZZ(self.power/3) * self.threePowerTorsionPoint

phi = EllipticCurveIsogeny( E, R )
return EllipticCurveWithAThreeTorsionPoint( phi.codomain()
, [ phi( point ) for point in self.points ]
, phi( self.threePowerTorsionPoint )
, ZZ( self.power/3 ) )

def __str__(self):
"""implement a simple string / print routine on an instance of this class
"""
s =  ( 'IElliptic (urve y^2 = x^3 + a4 x + a6 with\n    a4 = %s\n    a6 = %s\n'
% ( self.E.a4(), self.E.a6() ) )
s += ( '  +orsion point of order 3^%s:\n    %s\n'
% ( self.exponent, self.threePowerTorsionPoint ) )
if self.exponent:
s += ( '  +he isogeny will be built w.r.t. the point\n    %s\n\n'
% ( ( 3^self.exponent-1 ) * self.threePowerTorsionPoint ) )
s += '  Further points on the curve are:\n'
for point in self.points:
s += '    %s\n' % str(point)

return s


This is the class, now let us see it in action. The first examples are over $\mathbb F_p$, the field with $p= 2^{23} \cdot 3^{16} -1 = 361102068154367$ (instead of the later $p^2$) elements.

p     = 2^23 * 3^16 - 1    # 361102068154367
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)

E1 = EllipticCurve( Fp , [1,0] )
E2 = EllipticCurve( Fp2, [1,0] )

P0 = (2^23) * E1.gens()[0]
x0, y0 = P0.xy()
P0 = E2( (  x0,   y0 ) )    #   P0,  but over Fp2
Q0 = E2( ( -x0, j*y0 ) )    # j*P0, also over Fp2
R0 = P0 + 210*Q0            # some linear combination of the two

w0 = EllipticCurveWithAThreeTorsionPoint( E2, [P0, Q0], R0 )

w1_as_posted = w0.imageAsPosted()
w1           = w0.image()


And we print the information for the two instances:

print w1_as_posted
print w1


and get twice the same information on the image object:

IElliptic (urve y^2 = x^3 + a4 x + a6 with
a4 = 51991424516405
a6 = 106070704802680
+orsion point of order 3^15:
(302203874623956*j + 20241918426644 : 322186935621588*j + 312811758520168 : 1)
+he isogeny will be built w.r.t. the point
(302203874623956*j + 20241918426644 : 38915132532779*j + 48290309634199 : 1)

Further points on the curve are:
(160564418087326 : 134713587156739 : 1)
(71255284428256 : 70954627395972*j : 1)


(The print has some abused characters, this suggests that we are not using a native sage class.)

The result encourages us to construct the chain of isogenies, as long as the point stored in the attribute threePowerTorsionPoint has an exponent $>1$.

We use first the isogeny building method:

w = w0
while w.exponent > 0:
print ( "exponent = %2s :: building image for y^2 = x^3 + (%s) x + (%s)"
% ( w.exponent, w.E.a4(), w.E.a6() ) )
w = w.imageAsPosted()
print w


This gives:

exponent = 16 :: building image for y^2 = x^3 + (1) x + (0)
exponent = 15 :: building image for y^2 = x^3 + (51991424516405) x + (106070704802680)
exponent = 14 :: building image for y^2 = x^3 + (101255065698811*j + 101998360721687) x + (330431961294422*j + 25640531191555)
exponent = 13 :: building image for y^2 = x^3 + (352041190233198*j + 16843401014417) x + (119053808071932*j + 53355491076829)
exponent = 12 :: building image for y^2 = x^3 + (112624919445751*j + 82843834368224) x + (124262375340213*j + 302574096852259)
exponent = 11 :: building image for y^2 = x^3 + (50168658132117*j + 67647704823496) x + (170114631740667*j + 36753695819527)
exponent = 10 :: building image for y^2 = x^3 + (30928695707655*j + 183195208012603) x + (265120391455965*j + 347451132891807)
exponent =  9 :: building image for y^2 = x^3 + (277462814931074*j + 213640330656493) x + (323487652702660*j + 204825182375950)
exponent =  8 :: building image for y^2 = x^3 + (173601874168728*j + 76247491125478) x + (25871313617758*j + 356092938050180)
exponent =  7 :: building image for y^2 = x^3 + (250624169642861*j + 230365614922372) x + (252436205330755*j + 327691346106738)
exponent =  6 :: building image for y^2 = x^3 + (33619484794617*j + 158677295637010) x + (308521227821834*j + 149713210013437)
exponent =  5 :: building image for y^2 = x^3 + (348558022695077*j + 288578526046131) x + (253490762074226*j + 317654848250943)
exponent =  4 :: building image for y^2 = x^3 + (82379787001615*j + 138028939196786) x + (231911717052998*j + 161383009984030)
exponent =  3 :: building image for y^2 = x^3 + (194536414255940*j + 314065693503502) x + (237925635030217*j + 234644912327985)
exponent =  2 :: building image for y^2 = x^3 + (300393944202888*j + 41823399542207) x + (206277166539057*j + 273449608961593)
exponent =  1 :: building image for y^2 = x^3 + (83910679013649*j + 38669304005764) x + (149931244633737*j + 174603264088713)
IElliptic (urve y^2 = x^3 + a4 x + a6 with
a4 = 126073008059375*j + 104049819145236
a6 = 359176734976232*j + 166096128313403
+orsion point of order 3^0:
(0 : 1 : 0)
Further points on the curve are:
(269391644360339*j + 37967495842715 : 187578019259485*j + 238323252588410 : 1)
(344452009602701*j + 88416939333687 : 105343006059252*j + 302867275180483 : 1)


The same is obtained by using the ad-hoc method image:

w = w0
while w.exponent > 0:
print ( "exponent = %2s :: building image for y^2 = x^3 + (%s) x + (%s)"
% ( w.exponent, w.E.a4(), w.E.a6() ) )
w = w.image()
print w


This should be enough for the confidence. We finally apply the method image recursively for the situation in the posted question:

p     = 2^372 * 3^239 - 1
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)

x0 = 4359917396849101231053336763700300892915096700013704210194781457801412731643988367389870886884336453245156775454336249199185654250159051929975600857047173121187832546031604804277991148436536445770452624367894371450077315674371
y0 = 106866937607440797536385002617766720826944674650271400721039514250889186719923133049487966730514682296643039694531052672873754128006844434636819566554364257913332237123293860767683395958817983684370065598726191088239028762772

E  = EllipticCurve( Fp2, [1,0] )
P0 = E( (  x0,   y0 ) )
Q0 = E( ( -x0, j*y0 ) )
R0 = P0 + 210*Q0

w0 = EllipticCurveWithAThreeTorsionPoint( E, [P0, Q0], R0 )

w = w0
while w.exponent > 0:
print ( "exponent = %2s :: building image for y^2 = x^3 + (%s) x + (%s)"
% ( w.exponent, w.E.a4(), w.E.a6() ) )
w = w.image()
print w


This finally builds:

IElliptic (urve y^2 = x^3 + a4 x + a6 with
a4 = 9941905577249112621469246065425110821350719912107274783202186679465880266108084272870578652169622033033287734603565106172021579672428567992207417984928001915443792506581083414819631162355315474660473917599782239436793726843839*j + 2946724050834478985321445898239727252831818838289417899276985485853434756704829638077752203048021007103400137388580807566509509465255459102293317130002209435917122438527365276268239226236808877746332435971858250480260987080303
a6 = 4542539462589491958039849159619489262298226045275930640172064661787827291349433884124124339137609032724837343401235162947672566178714780903635487084785206544990561968525198653520460001313626004064990324340749794219991011726709*j + 6539896041012828166415736540021157859063938080758824271318838813777024260451932027047549259029251959252686859384332609437478578564216583785613603428933665799255028118170072559061461339831670487442990678076875312932414988112166
+orsion point of order 3^0:
(0 : 1 : 0)
Further points on the curve are:
(840330075266097098012742632864003386484260836813506546489726156765269120850766886866702513728233478613303948408691476297176547452485103585998818504139217010240338448809626375064720694578374062040485739348323048513110745385248*j + 3631117883451682334564713665100760377954652702481381082145400116289119393182640553809752129975570766526660183243812887136161277347553039054565653848432209570658617836874423389756886496530827346326347812938715566043877640930134 : 6752074772817769174079436029355327077672427260077485064638913088045654990683937363103311433711692705103529919542578009074989001433108466176516106774717123871276123524166867491669634446650573483815646122269973602608584045745385*j + 5300429866790144453496118748964075869721327676415610914943214711039869462486555654575277845491673873113512648830699054430011366710424869623809450389946736293221040732299394727069618963598635154541072518931890214684120653405947 : 1)
(9504884181665882619221119333795526046079464038402854413734436783207374710864936946472828883342585080762306046522952381796806852003567363266196038656783536742634373224984983661213652542712394821902578936678779599785110962732924*j + 6774194310353160483626699390014464172817940562986390128871705048641081983605314292195991865678096793894309639849495491019251134893033290390022357832137438886247925463306315515096565278042249805466199618653007280483904791566515 : 9204574484896333902716788631110480691421118707703587981927032022515666045032512213688539029367801791890984612552398053445102850985194349538412058985258007202160147628762221975794700299371203389505876094942576447710531845411107*j + 5442466631823826726067031217288108231018778098280360211573310329856059797916949670046581410065913468909932930890429505261541023989713829505640917502609310463204092834265761312631970799287827944796524047962553671995719042559911 : 1)


Note: We can check for the last two points $P, Q$, the expected equality: $P+210\cdot Q=O$, recursively inherited from $P_0+201\cdot Q_0=R_0$, since after the chain of isogenies $R_0$ is mapped to the origin of the final curve:

P, Q = w.points
print (P + 210*Q ).is_zero()


This evaluates to True.

Note: Thanks for the question! Special examples like this one are important to play and work with the structure! I learned a lot in the given framework.

Note: It was relatively simple to build the ad-hoc code, a slight adaptation may / should solve the general case, i will take a look to see where the present sage path of computing the isogeny gets frozen with the big numbers.

more

Thank you very much for the long and beautiful answer. Though, your answers works for 3-isogeny (I also tested), I do not see what I should modify to make it work for other isogenies, let's say for 2-isogeny, 4-isogeny or 5-isogeny? While I see you have also used different formula for $\phi$ than the one in the paper.

( 2018-02-02 15:02:53 -0600 )edit

## Stats

Seen: 95 times

Last updated: Jan 25