Isogeny computation does not finish in Sage

edit

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.

LATER EDIT

Sorry, today, three weeks later, i saw that there is still an open issue here. So i am trying to write all details needed to construct the isogeny w.t.r. an odd order generator of the kernel. (This, combined with the sage computation for even order generators of the kernel may be enough to proceed.) Here i address only the case of an odd generator point for a cyclic kernel of odd degree. In the above article, the things are presented as follows. The formulas in the sequel can be easily implemented. (I also give a check for them. But the class covering them may not be very inspired, its only reason is to give a quick check.)

Let us fix again an elliptic curve $E$ with the standard Weierstrass equation $$E\ :\ y^2 = x^3+ax+b\ , \ a,b\in F\ . $$ ($F$ being the field of definition of $E$.) Let $R$ be a point of odd order $d$. Let $G$ be the set of $R$-multiples, $G=\{0\cdot R, 1\cdot R, \dots,(d-1)\cdot R\}$. Let $G^\times$ be this set after we eliminate $O$, the zero point (at infinity) from $G$ Let $A$ be a rational point on $E(F)$, Then the formula for the isogeny is as follows: $$\phi(x_A,y_A)=\left(\ x_A+\sum_{P\in G^\times}(x_{A+P}-x_P)\ ,\ y_A+\sum_{P\in G^\times}(y_{A+P}-y_P)\ \right)\ .$$ This formula becomes simpler in case we can write $$G^\times =G_+\sqcup G_-\ ,$$ where a point $P$ is in $G_+$ iff $-P\in G_-$, and conversely.

We do so in our case, $d$ is odd, so $d-1$ is even, we take $G_+$ to be the set of all $mP$ for $m=1,2,\dots$ up to the half $(d-1)/2$. Then $G_-$ is the set of all other multiples of $P$ (not equal to $O$).

We get the following formula, which involves a $P\in G_+$: \begin{align} g_P^x &= g'_x(P) = (x^3 +ax+b)'_x(P)= 3x_P^2 +a\ ,\\ g_P^y &= g'_y(P) = (-y^2)'_y(P) = -2y_P\ ,\\ v_P &= 2g_P^x\ ,\\ u_P &= (g_P^y)^2\ ,\\ v &= \sum_{P\in G_+}v_P\ ,\\ w &= \sum_{P\in G_+}(u_P + x_Pv_P)\ .\\ 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( \quad x+\sum_{P\in G_+}\left(\frac{v_P}{x-x_P}-\frac{u_P}{(x-x_P)^2}\right) \quad ,\quad y+\sum_{P\in G_+}\left(\frac{v_P y}{(x-x_P)^2} - \frac{2u_P y}{(x-x_P)^3}\right)\quad \right)\ ,\\ \end{align} So it is "the same formula" as in the case $d=3$, the difference is that we sum over more multiples of the generator $R$.

Let us test this, we consider the following quickly adapted code:

class EllipticCurveWithAnOddOrderTorsionPoint( object ):
    """Class putting together in a bag:
    - an elliptic curve,
    - with a point R of odd order, d,
    - 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 = []    # and we will append

        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    # and we further add
            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 EllipticCurveWithAnOddOrderTorsionPoint(
            EE
            , EE(0)
            , d = 1 
            , points = [ phi(point) for point in self.points ] )


    def imageStandard( self ):
        """This does the same 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 EllipticCurveWithAnOddOrderTorsionPoint(
            EE
            , EE(0)
            , d=1
            , points = [ phi( point ) for point in self.points ] )


# ============================== Test ==============================
F = GF(2017)
E = EllipticCurve( GF(2017), [1,0] )
ord = E.order()
print "%s\nhas order %s = %s" % ( E, ord, ord.factor() )
G = E.gens()[0]
S = E.gens()[1]
T = E.random_point()
print "One generator is G = %s, it has order %s" % ( G, G.order().factor() )
R = 100*G
print "Consider R = 100G = %s" % str( R.xy() )
print "Consider S = %s as an other point on the curve." % str( S.xy() )
print "Consider T = %s as an other point on the curve." % str( T.xy() )

X = EllipticCurveWithAnOddOrderTorsionPoint( E, R, None, [S, T] )

print "\nIsogeny computation using ad-hoc formula ..." 
XX = X.image()
print "The codomain of isogeny is", XX.E
print "The mapped points are:", XX.points

print "\nIsogeny computation using sage constructions ..." 
XX = X.imageStandard()
print "The codomain of isogeny is", XX.E
print "The mapped points are:", XX.points

Results this time (T is a random point):

Elliptic Curve defined by y^2 = x^3 + x over Finite Field of size 2017
has order 2000 = 2^4 * 5^3
One generator is G = (1615 : 1881 : 1), it has order 2^2 * 5^3
Consider R = 100G = (656, 961)
Consider S = (1141, 1768) as an other point on the curve.
Consider T = (19, 1236) as an other point on the curve.
Computed: Order of R is 5

Isogeny computation using ad-hoc formula ...
Computing phi( (1141 : 1768 : 1) )....
Computing phi( (19 : 1236 : 1) )....
The codomain of isogeny is Elliptic Curve defined by y^2 = x^3 + 1455*x over Finite Field of size 2017
The mapped points are: [(644 : 284 : 1), (345 : 685 : 1)]

Isogeny computation using sage constructions ...
The codomain of isogeny is Elliptic Curve defined by y^2 = x^3 + 1455*x over Finite Field of size 2017
The mapped points are: [(644 : 284 : 1), (345 : 685 : 1)]

A final note: It is thus simple to compute images of points using python functions, i suppose sage needs such a long time for the computations, since there may be also computed a the "algebraic version" of the python function.