source: sage/schemes/elliptic_curves/ell_generic.py @ 7058:667d6986cd44

Revision 7058:667d6986cd44, 49.3 KB checked in by William Stein <wstein@…>, 6 years ago (diff)

trac #250 -- strange elliptic curve printing
Line 
1"""
2Elliptic curves over a general ring
3
4EXAMPLES:
5We construct an elliptic curve over an elaborate base ring:
6    sage: p = 97; a=1; b=3
7    sage: R, u = PolynomialRing(GF(p), 'u').objgen()
8    sage: S, v = PolynomialRing(R, 'v').objgen()
9    sage: T = S.fraction_field()
10    sage: E = EllipticCurve(T, [a, b]); E
11    Elliptic Curve defined by y^2  = x^3 + x + 3 over Fraction Field of Univariate Polynomial Ring in v over Univariate Polynomial Ring in u over Finite Field of size 97
12    sage: latex(E)
13    y^2  = x^3 + x + 3
14"""
15
16#*****************************************************************************
17#       Copyright (C) 2005 William Stein <wstein@gmail.com>
18#
19#  Distributed under the terms of the GNU General Public License (GPL)
20#
21#    This code is distributed in the hope that it will be useful,
22#    but WITHOUT ANY WARRANTY; without even the implied warranty of
23#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
24#    General Public License for more details.
25#
26#  The full text of the GPL is available at:
27#
28#                  http://www.gnu.org/licenses/
29#*****************************************************************************
30
31
32import math
33
34from sage.rings.all import MPolynomialRing, PolynomialRing
35
36import sage.plot.all as plot
37
38import sage.rings.arith as arith
39import sage.rings.all as rings
40import sage.rings.number_field as number_field
41from sage.rings.all import is_Infinite
42import sage.misc.misc as misc
43import sage.misc.latex as latex
44import sage.modular.modform as modform
45import sage.functions.transcendental as transcendental
46
47# Schemes
48import sage.schemes.generic.projective_space as projective_space
49import sage.schemes.generic.homset as homset
50
51import ell_point
52import constructor
53import formal_group
54
55factor = arith.factor
56sqrt = math.sqrt
57exp = math.exp
58mul = misc.mul
59next_prime = arith.next_prime
60
61oo = rings.infinity       # infinity
62O = rings.O         # big oh
63
64import sage.schemes.plane_curves.projective_curve as plane_curve
65
66def is_EllipticCurve(x):
67    """
68    EXAMPLES:
69        sage: E = EllipticCurve([1,2,3/4,7,19])
70        sage: is_EllipticCurve(E)
71        True
72        sage: is_EllipticCurve(0)
73        False
74    """
75    return isinstance(x, EllipticCurve_generic)
76
77class EllipticCurve_generic(plane_curve.ProjectiveCurve_generic):
78    """
79    Elliptic curve over a generic base ring.
80
81    EXAMPLES:
82        sage: E = EllipticCurve([1,2,3/4,7,19]); E
83        Elliptic Curve defined by y^2 + x*y + 3/4*y = x^3 + 2*x^2 + 7*x + 19 over Rational Field
84        sage: loads(E.dumps()) == E
85        True
86        sage: E = EllipticCurve([1,3])
87        sage: P = E([-1,1,1])
88        sage: -5*P
89        (179051/80089 : -91814227/22665187 : 1)
90    """
91    def __init__(self, ainvs, extra=None):
92        if extra != None:   # possibility of two arguments
93            K, ainvs = ainvs, extra
94        else:
95            K = ainvs[0].parent()
96        assert len(ainvs) == 2 or len(ainvs) == 5
97        self.__base_ring = K
98        ainvs = [K(x) for x in ainvs]
99        if len(ainvs) == 2:
100            ainvs = [K(0),K(0),K(0)] + ainvs
101        self.__ainvs = ainvs
102        if self.discriminant() == 0:
103            raise ArithmeticError, \
104                  "Invariants %s define a singular curve."%ainvs
105        PP = projective_space.ProjectiveSpace(2, K, names='xyz');
106        x, y, z = PP.coordinate_ring().gens()
107        a1, a2, a3, a4, a6 = ainvs
108        f = y**2*z + (a1*x + a3*z)*y*z \
109            - (x**3 + a2*x**2*z + a4*x*z**2 + a6*z**3)
110        plane_curve.ProjectiveCurve_generic.__init__(self, PP, f)
111        # TODO: cleanup, are these two point classes redundant?
112        if K.is_field():
113            self._point_morphism_class = self._point_class = ell_point.EllipticCurvePoint_field
114        else:
115            self._point_morphism_class = self._point_class = ell_point.EllipticCurvePoint
116
117    def _defining_params_(self):
118        return (self.__base_ring, self.__ainvs)
119
120    def _repr_(self):
121        """
122        String representation of elliptic curve.
123
124        EXAMPLES:
125            sage: EllipticCurve([1,2,3,4,5])
126            Elliptic Curve defined by y^2 + x*y + 3*y = x^3 + 2*x^2 + 4*x + 5 over Rational Field
127
128            sage: R.<x> = QQ['x']
129            sage: K.<a> = NumberField(x^3-17)
130            sage: EllipticCurve([a^2-3, -2/3*a + 3])
131            Elliptic Curve defined by y^2  = x^3 + (a^2-3)*x + (-2/3*a+3) over Number Field in a
132            with defining polynomial x^3 - 17           
133        """
134        #return "Elliptic Curve with a-invariants %s over %s"%(self.ainvs(), self.base_ring())
135        b = self.ainvs()
136        #return "y^2 + %s*x*y + %s*y = x^3 + %s*x^2 + %s*x + %s"%\
137        #       (a[0], a[2], a[1], a[3], a[4])
138        a = [z._coeff_repr() for z in b]
139        s = "Elliptic Curve defined by "
140        s += "y^2 "
141        if a[0] == "-1":
142            s += "- x*y "
143        elif a[0] == '1':
144            s += "+ x*y "
145        elif b[0]:
146            s += "+ %s*x*y "%a[0]
147        if a[2] == "-1":
148            s += " - y"
149        elif a[2] == '1':
150            s += "+ y"
151        elif b[2]:
152            s += "+ %s*y"%a[2]
153        s += " = x^3 "
154        if a[1] == "-1":
155            s += "- x^2 "
156        elif a[1] == '1':
157            s += "+ x^2 "
158        elif b[1]:
159            s += "+ %s*x^2 "%a[1]
160        if a[3] == "-1":
161            s += "- x "
162        elif a[3] == '1':
163            s += "+ x "
164        elif b[3]:
165            s += "+ %s*x "%a[3]
166        if a[4] == '-1':
167            s += "-1 "
168        elif a[4] == '1':
169            s += "+1 "
170        elif b[4]:
171            s += "+ %s "%a[4]
172        s = s.replace("+ -","- ")
173        s += "over %s"%self.base_ring()
174        return s
175
176    def _latex_(self):
177        b = self.ainvs()
178        a = [z._latex_coeff_repr() for z in b]
179        s = "y^2 "
180        if a[0] == '-1':
181            s += "- xy "
182        elif a[0] == '1':
183            s += "+ xy "
184        elif b[0]:
185            s += "+ %sxy "%a[0]
186        if a[2] == '-1':
187            s += " - y"
188        elif a[2] == '1':
189            s += "+ y"
190        elif b[2]:
191            s += "+ %sy"%a[2]
192        s += " = x^3 "
193        if a[1] == '-1':
194            s += "- x^2 "
195        elif a[1] == '1':
196            s += "+ x^2 "
197        elif b[1]:
198            s += "+ %sx^2 "%a[1]
199        if a[3] == '-1':
200            s += "- x "
201        elif a[3] == '1':
202            s += "+ x "
203        elif b[3]:
204            s += "+ %sx "%a[3]
205        if a[4] == '-1':
206            s += "-1 "
207        elif a[4] == '1':
208            s += "+1 "
209        elif b[4]:
210            s += "+ %s "%a[4]
211        s = s.replace("+ -","- ")
212        return s
213
214    def _pari_init_(self):
215        return 'ellinit([%s])'%(','.join([x._pari_init_() for x in self.ainvs()]))
216
217    def _magma_init_(self):
218        return 'EllipticCurve([%s])'%(','.join([x._magma_init_() for x in self.ainvs()]))
219
220    def _symbolic_(self, SR):
221        r"""
222        Many elliptic curves can be converted into a symbolic expression
223        using the \code{symbolic_expression} command.
224       
225        EXAMPLES:
226        We find a torsion point on 11a.
227            sage: E = EllipticCurve('11a')
228            sage: E.torsion_subgroup().gens()
229            ((5 : 5 : 1),)
230
231        We find the corresponding symbolic equality:
232            sage: eqn = symbolic_expression(E); eqn
233            y^2 + y == x^3 - x^2 - 10*x - 20
234            sage: print eqn
235                                      2        3    2
236                                     y  + y == x  - x  - 10 x - 20
237
238        We verify that the given point is on the curve:
239            sage: eqn(x=5,y=5)
240            30 == 30
241            sage: bool(eqn(x=5,y=5))
242            True
243
244        We create a single expression:
245            sage: F = eqn.lhs() - eqn.rhs(); print F
246                                      2        3    2
247                                     y  + y - x  + x  + 10 x + 20
248            sage: y = var('y')
249            sage: print F.solve(y)
250            [
251                                  3      2
252                        - sqrt(4 x  - 4 x  - 40 x - 79) - 1
253                    y == -----------------------------------
254                                         2,
255                                 3      2
256                         sqrt(4 x  - 4 x  - 40 x - 79) - 1
257                     y == ---------------------------------
258                                         2
259            ]
260
261        You can also solve for x in terms of y, but the result is horrendous.
262        Continuing with the above example, we can explicitly find points
263        over random fields by substituting in values for x:
264       
265            sage: v = F.solve(y)[0].rhs()   
266            sage: print v
267                                            3      2
268                                  - sqrt(4 x  - 4 x  - 40 x - 79) - 1
269                                  -----------------------------------
270                                                   2
271            sage: v(3)
272            (-sqrt(127)*I - 1)/2
273            sage: v(7)
274            (-sqrt(817) - 1)/2
275            sage: v(-7)
276            (-sqrt(1367)*I - 1)/2
277            sage: v(sqrt(2))
278            (-sqrt(-32*sqrt(2) - 87) - 1)/2
279
280        We can even do arithmetic with them, as follows:
281            sage: E2 = E.change_ring(SR); E2
282            Elliptic Curve defined by y^2 + y = x^3 - x^2 - 10*x - 20 over Symbolic Ring
283            sage: P = E2.point((3, v(3), 1), check=False)
284            sage: P
285            (3 : (-sqrt(127)*I - 1)/2 : 1)
286            sage: P + P
287            (-756/127 : (sqrt(127)*I + 1)/2 + 12507*I/(127*sqrt(127)) - 1 : 1)
288
289        We can even throw in a transcendental:
290            sage: w = E2.point((pi,v(pi),1), check=False); w
291            (pi : (-sqrt(4*pi^3 - 4*pi^2 - 40*pi - 79) - 1)/2 : 1)
292            sage: 2*w
293            ((3*pi^2 - 2*pi - 10)^2/(4*pi^3 - 4*pi^2 - 40*pi - 79) - 2*pi + 1 : (sqrt(4*pi^3 - 4*pi^2 - 40*pi - 79) + 1)/2 - (3*pi^2 - 2*pi - 10)*(-(3*pi^2 - 2*pi - 10)^2/(4*pi^3 - 4*pi^2 - 40*pi - 79) + 3*pi - 1)/sqrt(4*pi^3 - 4*pi^2 - 40*pi - 79) - 1 : 1)
294        """
295        a = [SR(x) for x in self.a_invariants()]
296        x, y = SR.var('x, y')
297        return y**2 + a[0]*x*y + a[2]*y == x**3 + a[1]*x**2 + a[3]*x + a[4]
298       
299    def __cmp__(self, other):
300        if not isinstance(other, EllipticCurve_generic):
301            return -1
302        t = cmp(self.base_ring(), other.base_ring())
303        if t:
304            return t
305        return cmp(self.ainvs(), other.ainvs())
306
307    def __contains__(self, P):
308        """
309        Returns True if and only if P defines is a point on the
310        elliptic curve.  P just has to be something that can be
311        coerced to a point.
312
313        EXAMPLES:
314            sage: E = EllipticCurve([0, 0, 1, -1, 0])
315            sage: (0,0) in E
316            True
317            sage: (1,3) in E
318            False
319            sage: E = EllipticCurve([GF(7)(0), 1])
320            sage: [0,0] in E
321            False
322            sage: [0,8] in E
323            True
324            sage: P = E(0,8)
325            sage: P
326            (0 : 1 : 1)
327            sage: P in E
328            True
329        """
330        if not isinstance(P, ell_point.EllipticCurvePoint):
331            try:
332                P = self(P)
333            except TypeError:
334                return False
335        if P.curve() == self:
336            return True
337        x, y, a = P[0], P[1], self.ainvs()
338        return y**2 + a[0]*x*y + a[2]*y == x**3 + a[1]*x**2 + a[3]*x + a[4]
339
340    def __call__(self, *args, **kwds):
341        """
342        EXAMPLES:
343            sage: E = EllipticCurve([0, 0, 1, -1, 0])
344
345        The point at infinity, which is the 0 element of the group:
346            sage: E(0)
347            (0 : 1 : 0)
348
349        The origin is a point on our curve:
350            sage: P = E([0,0])
351            sage: P
352            (0 : 0 : 1)
353
354        The curve associated to a point:
355            sage: P.curve()
356            Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field
357
358        Points can be specified by given a 2-tuple or 3-tuple
359            sage: E([0,0])
360            (0 : 0 : 1)
361            sage: E([0,1,0])
362            (0 : 1 : 0)
363
364        Over a field, points are normalized so the right-most nonzero entry is 1:
365            sage: E(105, -69, 125)
366            (21/25 : -69/125 : 1)
367
368        We create points on an elliptic curve over a prime finite field.
369            sage: E = EllipticCurve([GF(7)(0), 1])
370            sage: E([2,3])
371            (2 : 3 : 1)
372            sage: E([0,0])
373            Traceback (most recent call last):
374            ...
375            TypeError: coordinates [0, 0, 1] do not define a point on Elliptic Curve defined by y^2  = x^3 +1 over Finite Field of size 7
376
377        We create a point on an elliptic curve over a number field.
378            sage: x = polygen(RationalField())
379            sage: K = NumberField(x**3 + x + 1, 'a'); a = K.gen()
380            sage: E = EllipticCurve([a,a])
381            sage: E
382            Elliptic Curve defined by y^2  = x^3 + a*x + a over Number Field in a with defining polynomial x^3 + x + 1
383            sage: E = EllipticCurve([K(1),1])
384            sage: E
385            Elliptic Curve defined by y^2  = x^3 + x +1 over Number Field in a with defining polynomial x^3 + x + 1
386            sage: P = E([a,0,1])
387            sage: P
388            (a : 0 : 1)
389            sage: P+P
390            (0 : 1 : 0)
391
392        Another example involving p-adics:
393            sage: E = EllipticCurve('37a1')
394            sage: P = E([0,0]); P
395            (0 : 0 : 1)
396            sage: R = pAdicField(3,20)
397            sage: Ep = E.base_extend(R); Ep
398            Elliptic Curve defined by y^2 + (1+O(3^20))*y = x^3 + (2+2*3+2*3^2+2*3^3+2*3^4+2*3^5+2*3^6+2*3^7+2*3^8+2*3^9+2*3^10+2*3^11+2*3^12+2*3^13+2*3^14+2*3^15+2*3^16+2*3^17+2*3^18+2*3^19+O(3^20))*x over 3-adic Field with capped relative precision 20
399            sage: Ep(P)
400            (0 : 0 : 1 + O(3^20))       
401        """
402        if len(args) == 1 and args[0] == 0:
403            R = self.base_ring()
404            return self.point([R(0),R(1),R(0)], check=False)
405        if isinstance(args[0],
406              (ell_point.EllipticCurvePoint_field, ell_point.EllipticCurvePoint)):
407            args = tuple(args[0])
408        return plane_curve.ProjectiveCurve_generic.__call__(self, *args, **kwds)
409   
410    def lift_x(self, x, all=False):
411        """
412        Given the x-coordinate of a point on the curve, use the defining
413        polynomial to find all affine points on this curve with the
414        given x-coordinate.
415       
416        EXAMPLES:
417            sage: E = EllipticCurve('37a'); E
418            Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field
419            sage: E.lift_x(1)
420            (1 : 0 : 1)
421            sage: E.lift_x(2)
422            (2 : 2 : 1)
423            sage: E.lift_x(1/4, all=True)
424            [(1/4 : -3/8 : 1), (1/4 : -5/8 : 1)]
425           
426        There are no rational points with x-cordinate 3.
427            sage: E.lift_x(3)
428            Traceback (most recent call last):
429            ...
430            ValueError: No point with x-coordinate 3 on Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field
431           
432        However, there are two such points in $E(\R)$:
433            sage: E.change_ring(RR).lift_x(3, all=True)
434            [(3.00000000000000 : 4.42442890089805 : 1.00000000000000), (3.00000000000000 : -5.42442890089805 : 1.00000000000000)]
435           
436        And of course it always works in $E(\C)$:
437            sage: E.change_ring(RR).lift_x(.5, all=True)
438            []
439            sage: E.change_ring(CC).lift_x(.5)
440            (0.500000000000000 : -0.500000000000000 + 0.353553390593274*I : 1.00000000000000)
441           
442        We can perform these operations over finite fields too:
443            sage: E = E.change_ring(GF(17)); E
444            Elliptic Curve defined by y^2 + y = x^3 + 16*x over Finite Field of size 17
445            sage: E.lift_x(7)
446            (7 : 11 : 1)
447            sage: E.lift_x(3)
448            Traceback (most recent call last):
449            ...
450            ValueError: No point with x-coordinate 3 on Elliptic Curve defined by y^2 + y = x^3 + 16*x over Finite Field of size 17
451           
452        Note that there is only one lift with x-coordinate 10 in $E(\F_{17})$.
453            sage: E.lift_x(10, all=True)
454            [(10 : 8 : 1)]
455           
456        We can lift over more exotic rings too.
457            sage: E = EllipticCurve('37a');
458            sage: E.lift_x(pAdicField(17, 5)(6))
459            (6 + O(17^5) : 2 + 16*17 + 16*17^2 + 16*17^3 + 16*17^4 + O(17^5) : 1 + O(17^5))
460            sage: K.<t> = PowerSeriesRing(QQ, 't', 5)
461            sage: E.lift_x(1+t)
462            (1 + t : 2*t - t^2 + 5*t^3 - 21*t^4 + O(t^5) : 1)
463            sage: K.<a> = GF(16)
464            sage: E = E.change_ring(K)
465            sage: E.lift_x(a^3)
466            (a^3 : a^3 + a : 1)
467           
468           
469        AUTHOR: Robert Bradshaw, 2007-04-24
470           
471        TEST:
472            sage: E = EllipticCurve('37a').weierstrass_model().change_ring(GF(17))
473            sage: E.lift_x(3, all=True)
474            []
475            sage: E.lift_x(7, all=True)
476            [(7 : 3 : 1), (7 : 14 : 1)]
477        """
478        a1, a2, a3, a4, a6 = self.ainvs()
479        f = ((x + a2) * x + a4) * x + a6
480        K = self.base_ring()
481        x += K(0)
482        one = x.parent()(1)
483        if a1.is_zero() and a3.is_zero():
484            if f.is_square():
485                if all:
486                    ys = f.sqrt(all=True)
487                    return [self.point([x, y, one], check=False) for y in ys]
488                else:
489                    return self.point([x, f.sqrt(), one], check=False)
490        else:
491            b = (a1*x + a3)
492            D = b*b + 4*f
493            if K.characteristic() == 2:
494                R = PolynomialRing(K, 'y')
495                F = R([-f,b,1])
496                ys = F.roots()
497                if all:
498                    return [self.point([x, y[0], one], check=False) for y in ys]
499                elif len(ys) > 0:
500                    return self.point([x, ys[0][0], one], check=False)
501            elif D.is_square():
502                if all:
503                    return [self.point([x, (-b+d)/2, one], check=False) for d in D.sqrt(all=True)]
504                else:
505                    return self.point([x, (-b+D.sqrt())/2, one], check=False)
506        if all:
507            return []
508        else:
509            raise ValueError, "No point with x-coordinate %s on %s"%(x, self)
510
511    def _homset_class(self, *args, **kwds):
512        return homset.SchemeHomsetModule_abelian_variety_coordinates_field(*args, **kwds)
513
514    def __getitem__(self, n):
515        """
516        """
517        raise NotImplementedError, "not implemented."
518
519    def __is_over_RationalField(self):
520        return isinstance(self.base_ring(), rings.RationalField)
521
522    def change_ring(self, R):
523        """
524        Return the elliptic curve defined by coercing the a-invariants
525        of this elliptic curve into the ring R.
526
527        INPUT:
528            R -- ring
529
530        OUTPUT:
531            an elliptic curve
532
533        EXAMPLES:
534            sage: E = EllipticCurve([0, 0, 1, -1, 0])
535            sage: E.change_ring(GF(3))
536            Elliptic Curve defined by y^2 + y = x^3 + 2*x over Finite Field of size 3
537        """
538        return constructor.EllipticCurve(R, [R(a) for a in self.ainvs()])
539
540    def is_on_curve(self, x, y):
541        """
542        Returns True if the (x,y) is an affine point on this curve.
543        """
544        a = self.ainvs()
545        return y**2 +a[0]*x*y + a[2]*y == x**3 + a[1]*x**2 + a[3]*x + a[4]
546
547    def a_invariants(self):
548        """
549        The a-invariants of this elliptic curve.
550
551        EXAMPLES:
552            sage: E = EllipticCurve([1,2,3,4,5])
553            sage: E.a_invariants()
554            [1, 2, 3, 4, 5]
555            sage: E = EllipticCurve([0,1])
556            sage: E
557            Elliptic Curve defined by y^2  = x^3 +1 over Rational Field
558            sage: E.a_invariants()
559            [0, 0, 0, 0, 1]
560            sage: E = EllipticCurve([GF(7)(3),5])
561            sage: E.a_invariants()
562            [0, 0, 0, 3, 5]
563        """
564        return self.__ainvs
565
566    ainvs = a_invariants
567
568    def b_invariants(self):
569        """
570        The b-invariants of this elliptic curve.
571
572        EXAMPLES:
573            sage: E = EllipticCurve([0, -1, 1, -10, -20])
574            sage: E.b_invariants()
575            (-4, -20, -79, -21)
576            sage: E = EllipticCurve([-4,0])
577            sage: E.b_invariants()
578            (0, -8, 0, -16)
579
580            sage: E = EllipticCurve([1,2,3,4,5])
581            sage: E.b_invariants()
582            (9, 11, 29, 35)
583            sage: E.b2()
584            9
585            sage: E.b4()
586            11
587            sage: E.b6()
588            29
589            sage: E.b8()
590            35
591           
592        ALGORITHM: These are simple functions of the a invariants.
593
594        AUTHOR: William Stein, 2005-04-25
595        """
596        try:
597            return self.__b_invariants
598        except AttributeError:
599            a1,a2,a3,a4,a6 = tuple(self.ainvs())
600            self.__b_invariants = a1*a1 + 4*a2, \
601                                  a1*a3 + 2*a4, \
602                                  a3**2 + 4*a6, \
603                                  a1**2 * a6 + 4*a2*a6 - a1*a3*a4 + a2*a3**2 - a4**2
604            return self.__b_invariants
605
606    def b2(self):
607        """
608        EXAMPLES:
609            sage: E = EllipticCurve([1,2,3,4,5])
610            sage: E.b2()
611            9
612        """
613        try:
614            return self.__b_invariants[0]
615        except AttributeError:
616            return self.b_invariants()[0]
617       
618    def b4(self):
619        """
620        EXAMPLES:
621            sage: E = EllipticCurve([1,2,3,4,5])
622            sage: E.b4()
623            11
624        """
625        try:
626            return self.__b_invariants[1]
627        except AttributeError:
628            return self.b_invariants()[1]
629
630    def b6(self):
631        """
632        EXAMPLES:
633            sage: E = EllipticCurve([1,2,3,4,5])
634            sage: E.b6()
635            29
636        """
637        try:
638            return self.__b_invariants[2]
639        except AttributeError:
640            return self.b_invariants()[2]
641
642    def b8(self):
643        """
644        EXAMPLES:
645            sage: E = EllipticCurve([1,2,3,4,5])
646            sage: E.b8()
647            35
648        """
649        try:
650            return self.__b_invariants[3]
651        except AttributeError:
652            return self.b_invariants()[3]
653
654    def c_invariants(self):
655        """
656        The c-invariants of this elliptic curve.
657
658        EXAMPLES:
659            sage: E = EllipticCurve([0, -1, 1, -10, -20])
660            sage: E.c_invariants()
661            (496, 20008)
662            sage: E = EllipticCurve([-4,0])
663            sage: E.c_invariants()
664            (192, 0)
665
666        ALGORITHM: These are simple functions of the b invariants.
667
668        AUTHOR: William Stein, 2005-04-25
669        """
670        try:
671            return self.__c_invariants
672        except AttributeError:
673            b2,b4,b6,b8 = self.b_invariants()
674            self.__c_invariants = b2**2 - 24*b4,\
675                                  -b2**3 + 36*b2*b4 - 216*b6    # note: c6 is wrong in Silverman, but right in Cremona
676            return self.__c_invariants
677   
678    def c4(self):
679        try:
680            return self.__c_invariants[0]
681        except AttributeError:
682            pass
683        return self.c_invariants()[0]
684
685    def c6(self):
686        try:
687            return self.__c_invariants[1]
688        except AttributeError:
689            pass
690        return self.c_invariants()[1]
691
692
693    def base_extend(self, R):
694        return constructor.EllipticCurve(R, [R(a) for a in self.a_invariants()])
695
696    def base_ring(self):
697        """
698        Returns the base ring of the elliptic curves.
699       
700        EXAMPLES:
701            sage: E = EllipticCurve(GF(49, 'a'), [3,5])
702            sage: E.base_ring()
703            Finite Field in a of size 7^2
704           
705            sage: E = EllipticCurve([1,1])
706            sage: E.base_ring()
707            Rational Field
708
709        """
710       # NOT supported yet
711       #     sage: E = EllipticCurve(Z, [3,5])
712       #     sage: E.base_ring()
713       #     Integer Ring
714        return self.__base_ring
715
716    base_field = base_ring
717
718    def a1(self):
719        """
720        EXAMPLES:
721            sage: E = EllipticCurve([1,2,3,4,6])
722            sage: E.a1()
723            1
724        """
725        return self.__ainvs[0]
726   
727    def a2(self):
728        """
729        EXAMPLES:
730            sage: E = EllipticCurve([1,2,3,4,6])
731            sage: E.a2()
732            2
733        """
734        return self.__ainvs[1]
735   
736    def a3(self):
737        """
738        EXAMPLES:
739            sage: E = EllipticCurve([1,2,3,4,6])
740            sage: E.a3()
741            3
742        """
743        return self.__ainvs[2]
744   
745    def a4(self):
746        """
747        EXAMPLES:
748            sage: E = EllipticCurve([1,2,3,4,6])
749            sage: E.a4()
750            4
751        """
752        return self.__ainvs[3]
753   
754    def a6(self):
755        """
756        EXAMPLES:
757            sage: E = EllipticCurve([1,2,3,4,6])
758            sage: E.a6()
759            6
760        """
761        return self.__ainvs[4]
762
763    def gens(self):
764        raise NotImplementedError
765
766    def gen(self, i):
767        return self.gens()[i]
768       
769    def quadratic_twist(self, D):
770        """
771        Return the quadratic twist of this curve by D.
772
773        EXAMPLES:
774            sage: E = EllipticCurve([GF(1103)(1), 0, 0, 107, 340]); E
775            Elliptic Curve defined by y^2 + x*y  = x^3 + 107*x + 340 over Finite Field of size 1103
776            sage: E.quadratic_twist(-1)
777            Elliptic Curve defined by y^2  = x^3 + 770*x + 437 over Finite Field of size 1103
778        """
779        a = self.ainvs()
780        ap = [0, 0, 0]
781        ap.append(-27*D**2*a[0]**4 - 216*D**2*a[0]**2*a[1] + 648*D**2*a[0]*a[2] - 432*D**2*a[1]**2 + 1296*D**2*a[3])
782        ap.append(54*D**3*a[0]**6 + 648*D**3*a[0]**4*a[1]
783                  - 1944*D**3*a[0]**3*a[2] + 2592*D**3*a[0]**2*a[1]**2
784                  - 3888*D**3*a[0]**2*a[3] - 7776*D**3*a[0]*a[1]*a[2]
785                  + 3456*D**3*a[1]**3 - 15552*D**3*a[1]*a[3]
786                  + 11664*D**3*a[2]**2 + 46656*D**3*a[4])
787        import constructor
788        return constructor.EllipticCurve(self.base_ring(), ap)
789
790    def rst_transform(self, r, s, t):
791        """
792        Transforms the elliptic curve using the unimodular (u=1) transform with standard parameters [r,s,t].
793       
794        Returns the transformed curve.
795        """
796        ##Ported from John Cremona's code implementing Tate's algorithm.
797        (a1, a2, a3, a4, a6) = self.a_invariants()
798        a6 += r*(a4 + r*(a2 + r)) - t*(a3 + r*a1 + t)
799        a4 += -s*a3 + 2*r*a2 - (t + r*s)*a1 + 3*r*r - 2*s*t
800        a3 += r*a1 + 2*t
801        a2 += -s*a1 + 3*r - s*s
802        a1 += 2*s
803        return constructor.EllipticCurve([a1, a2, a3, a4, a6])
804   
805    def scale_curve(self, u):
806        """
807        Transforms the elliptic curve using scale factor $u$,
808        i.e. multiplies $c_i$ by $u^i$.
809       
810        Returns the transformed curve.
811        """
812        ##Ported from John Cremona's code implementing Tate's algorithm.
813        (a1,a2,a3,a4,a6) = self.a_invariants()
814        a1 *= u
815        a2 *= u^2
816        a3 *= u^3
817        a4 *= u^4
818        a6 *= u^6
819        return constructor.EllipticCurve([a1, a2, a3, a4, a6])
820
821    def discriminant(self):
822        """
823        Returns the discriminant of this elliptic curve.
824
825        EXAMPLES:
826            sage: E = EllipticCurve([0,0,1,-1,0])
827            sage: E.discriminant()
828            37
829            sage: E = EllipticCurve([0, -1, 1, -10, -20])
830            sage: E.discriminant()
831            -161051
832
833            sage: E = EllipticCurve([GF(7)(2),1])
834            sage: E.discriminant()
835            1
836        """
837        try:
838            return self.__discriminant
839        except AttributeError:
840            b2, b4, b6, b8 = self.b_invariants()
841            self.__discriminant = -b2**2*b8 - 8*b4**3 - 27*b6**2 + 9*b2*b4*b6
842            return self.__discriminant
843
844    def j_invariant(self):
845        """
846        Returns the j-invariant of this elliptic curve.
847
848        EXAMPLES:
849            sage: E = EllipticCurve([0,0,1,-1,0])
850            sage: E.j_invariant()
851            110592/37
852            sage: E = EllipticCurve([0, -1, 1, -10, -20])
853            sage: E.j_invariant()
854            -122023936/161051
855            sage: E = EllipticCurve([-4,0])
856            sage: E.j_invariant()
857            1728
858
859            sage: E = EllipticCurve([GF(7)(2),1])
860            sage: E.j_invariant()
861            1
862
863           
864        """
865        try:
866            return self.__j_invariant
867        except AttributeError:
868            c4, _ = self.c_invariants()
869            self.__j_invariant = c4**3 / self.discriminant()
870            return self.__j_invariant
871
872
873##     def pseudo_torsion_polynomial(self, n, x=None, cache=None):
874##         r"""
875##         Returns the n-th torsion polynomial (division polynomial), without
876##         the 2-torsion factor if n is even, as a polynomial in $x$.
877
878##         These are the polynomials $g_n$ defined in Mazur/Tate (``The p-adic
879##         sigma function''), but with the sign flipped for even $n$, so that
880##         the leading coefficient is always positive.
881
882##         The full torsion polynomials may be recovered as follows:
883##         \begin{itemize}
884##         \item $\psi_n = g_n$ for odd $n$.
885##         \item $\psi_n = (2y + a_1 x + a_3) g_n$ for even $n$.
886##         \end{itemize}
887
888##         Note that the $g_n$'s are always polynomials in $x$, whereas the
889##         $\psi_n$'s require the appearance of a $y$.
890
891##         SEE ALSO:
892##             -- torsion_polynomial()
893##             -- multiple_x_numerator()
894##             -- multiple_x_denominator()
895
896##         INPUT:
897##             n -- positive integer, or the special values -1 and -2 which
898##                  mean $B_6 = (2y + a_1 x + a_3)^2$ and $B_6^2$ respectively
899##                  (in the notation of Mazur/Tate).
900##             x -- optional ring element to use as the "x" variable. If x
901##                  is None, then a new polynomial ring will be constructed over
902##                  the base ring of the elliptic curve, and its generator will
903##                  be used as x. Note that x does not need to be a generator of
904##                  a polynomial ring; any ring element is ok. This permits fast
905##                  calculation of the torsion polynomial *evaluated* on any
906##                  element of a ring.
907##             cache -- optional dictionary, with integer keys. If the key m
908##                  is in cache, then cache[m] is assumed to be the value of
909##                  pseudo_torsion_polynomial(m) for the supplied x. New entries
910##                  will be added to the cache as they are computed.
911
912##         ALGORITHM:
913##             -- Recursion described in Mazur/Tate. The recursive formulae are
914##             evaluated $O((log n)^2)$ times.
915
916##         TODO:
917##             -- for better unity of code, it might be good to make the regular
918##             torsion_polynomial() function use this as a subroutine.
919
920##         AUTHORS:
921##             -- David Harvey (2006-09-24)
922
923##         EXAMPLES:
924##            sage: E = EllipticCurve("37a")
925##            sage: E.pseudo_torsion_polynomial(1)
926##            1
927##            sage: E.pseudo_torsion_polynomial(2)
928##            1
929##            sage: E.pseudo_torsion_polynomial(3)
930##            3*x^4 - 6*x^2 + 3*x - 1
931##            sage: E.pseudo_torsion_polynomial(4)
932##            2*x^6 - 10*x^4 + 10*x^3 - 10*x^2 + 2*x + 1
933##            sage: E.pseudo_torsion_polynomial(5)
934##            5*x^12 - 62*x^10 + 95*x^9 - 105*x^8 - 60*x^7 + 285*x^6 - 174*x^5 - 5*x^4 - 5*x^3 + 35*x^2 - 15*x + 2
935##            sage: E.pseudo_torsion_polynomial(6)
936##            3*x^16 - 72*x^14 + 168*x^13 - 364*x^12 + 1120*x^10 - 1144*x^9 + 300*x^8 - 540*x^7 + 1120*x^6 - 588*x^5 - 133*x^4 + 252*x^3 - 114*x^2 + 22*x - 1
937##            sage: E.pseudo_torsion_polynomial(7)
938##            7*x^24 - 308*x^22 + 986*x^21 - 2954*x^20 + 28*x^19 + 17171*x^18 - 23142*x^17 + 511*x^16 - 5012*x^15 + 43804*x^14 - 7140*x^13 - 96950*x^12 + 111356*x^11 - 19516*x^10 - 49707*x^9 + 40054*x^8 - 124*x^7 - 18382*x^6 + 13342*x^5 - 4816*x^4 + 1099*x^3 - 210*x^2 + 35*x - 3
939##            sage: E.pseudo_torsion_polynomial(8)
940##            4*x^30 - 292*x^28 + 1252*x^27 - 5436*x^26 + 2340*x^25 + 39834*x^24 - 79560*x^23 + 51432*x^22 - 142896*x^21 + 451596*x^20 - 212040*x^19 - 1005316*x^18 + 1726416*x^17 - 671160*x^16 - 954924*x^15 + 1119552*x^14 + 313308*x^13 - 1502818*x^12 + 1189908*x^11 - 160152*x^10 - 399176*x^9 + 386142*x^8 - 220128*x^7 + 99558*x^6 - 33528*x^5 + 6042*x^4 + 310*x^3 - 406*x^2 + 78*x - 5
941
942##            sage: E.pseudo_torsion_polynomial(18) % E.pseudo_torsion_polynomial(6) == 0
943##            True
944
945##           An example to illustrate the relationship with torsion points.
946##            sage: F = GF(11)
947##            sage: E = EllipticCurve(F, [0, 2]); E
948##            Elliptic Curve defined by y^2  = x^3 + 2 over Finite Field of size 11
949##            sage: f = E.pseudo_torsion_polynomial(5); f
950##            5*x^12 + x^9 + 8*x^6 + 4*x^3 + 7
951##            sage: f.factor()
952##            (5) * (x^2 + 5) * (x^2 + 2*x + 5) * (x^2 + 5*x + 7) * (x^2 + 7*x + 7) * (x^2 + 9*x + 5) * (x^2 + 10*x + 7)
953
954##           This indicates that the x-coordinates of all the 5-torsion points
955##           of $E$ are in $GF(11^2)$, and therefore the y-coordinates are in
956##           $GF(11^4)$.
957         
958##            sage: K = GF(11^4, 'a')
959##            sage: X = E.change_ring(K)
960##            sage: f = X.pseudo_torsion_polynomial(5)
961##            sage: x_coords = [root for (root, _) in f.roots()]; x_coords
962##            [10*a^3 + 4*a^2 + 5*a + 6,
963##             9*a^3 + 8*a^2 + 10*a + 8,
964##             8*a^3 + a^2 + 4*a + 10,
965##             8*a^3 + a^2 + 4*a + 8,
966##             8*a^3 + a^2 + 4*a + 4,
967##             6*a^3 + 9*a^2 + 3*a + 4,
968##             5*a^3 + 2*a^2 + 8*a + 7,
969##             3*a^3 + 10*a^2 + 7*a + 8,
970##             3*a^3 + 10*a^2 + 7*a + 3,
971##             3*a^3 + 10*a^2 + 7*a + 1,
972##             2*a^3 + 3*a^2 + a + 7,
973##             a^3 + 7*a^2 + 6*a]
974
975##           Now we check that these are exactly the x coordinates of the
976##           5-torsion points of E.
977##            sage: for x in x_coords:
978##            ...       y = (x**3 + 2).square_root()
979##            ...       P = X([x, y])
980##            ...       assert P.order(disable_warning=True) == 5
981
982##           todo: need to show an example where the 2-torsion is missing
983
984##         """
985##         if cache is None:
986##             cache = {}
987##         else:
988##             try:
989##                 return cache[n]
990##             except KeyError:
991##                 pass
992
993##         if x is None:
994##             x = rings.PolynomialRing(self.base_ring(), 'x').gen()
995
996##         b2, b4, b6, b8 = self.b_invariants()
997
998##         n = int(n)
999##         if n <= 4:
1000##             if n == -1:
1001##                 answer = 4*x**3 + b2*x**2 + 2*b4*x + b6
1002##             elif n == -2:
1003##                 answer = self.pseudo_torsion_polynomial(-1, x, cache) ** 2
1004##             elif n == 1 or n == 2:
1005##                 answer = x.parent()(1)
1006##             elif n == 3:
1007##                 answer = 3*x**4 + b2*x**3 + 3*b4*x**2 + 3*b6*x + b8
1008##             elif n == 4:
1009##                 answer = -self.pseudo_torsion_polynomial(-2, x, cache) + \
1010##                          (6*x**2 + b2*x + b4) * \
1011##                          self.pseudo_torsion_polynomial(3, x, cache)
1012##             else:
1013##                 raise ValueError, "n must be a positive integer (or -1 or -2)"
1014##         else:
1015##             if n % 2 == 0:
1016##                 m = (n-2) // 2
1017##                 g_mplus3 = self.pseudo_torsion_polynomial(m+3, x, cache)
1018##                 g_mplus2 = self.pseudo_torsion_polynomial(m+2, x, cache)
1019##                 g_mplus1 = self.pseudo_torsion_polynomial(m+1, x, cache)
1020##                 g_m      = self.pseudo_torsion_polynomial(m,   x, cache)
1021##                 g_mless1 = self.pseudo_torsion_polynomial(m-1, x, cache)
1022##                 answer = g_mplus1 * \
1023##                          (g_mplus3 * g_m**2 - g_mless1 * g_mplus2**2)
1024##             else:
1025##                 m = (n-1) // 2
1026##                 g_mplus2 = self.pseudo_torsion_polynomial(m+2, x, cache)
1027##                 g_mplus1 = self.pseudo_torsion_polynomial(m+1, x, cache)
1028##                 g_m      = self.pseudo_torsion_polynomial(m,   x, cache)
1029##                 g_mless1 = self.pseudo_torsion_polynomial(m-1, x, cache)
1030##                 B6_sqr   = self.pseudo_torsion_polynomial(-2, x, cache)
1031##                 if m % 2 == 0:
1032##                     answer = B6_sqr * g_mplus2 * g_m**3 - \
1033##                              g_mless1 * g_mplus1**3
1034##                 else:
1035##                     answer = g_mplus2 * g_m**3 - \
1036##                              B6_sqr * g_mless1 * g_mplus1**3
1037
1038##         cache[n] = answer
1039##         return answer
1040
1041
1042##     def multiple_x_numerator(self, n, x=None, cache=None):
1043##         r"""
1044##         Returns the numerator of the x-coordinate of the nth multiple of
1045##         a point, using torsion polynomials (division polynomials).
1046
1047##         The inputs n, x, cache are as described in pseudo_torsion_polynomial().
1048
1049##         The result is adjusted to be correct for both even and odd n.
1050
1051##         WARNING:
1052##           -- There may of course be cancellation between the numerator and
1053##           the denominator (multiple_x_denominator()). Be careful. For more
1054##           information on how to avoid cancellation, see Christopher Wuthrich's
1055##           thesis.
1056
1057##         SEE ALSO:
1058##           -- multiple_x_denominator()
1059
1060##         AUTHORS:
1061##            -- David Harvey (2006-09-24)
1062
1063##         EXAMPLES:
1064##           sage: E = EllipticCurve("37a")
1065##           sage: P = E.gens()[0]
1066##           sage: x = P[0]
1067
1068##           sage: (35*P)[0]
1069##           -804287518035141565236193151/1063198259901027900600665796
1070##           sage: E.multiple_x_numerator(35, x)
1071##           -804287518035141565236193151
1072##           sage: E.multiple_x_denominator(35, x)
1073##           1063198259901027900600665796
1074
1075##           sage: (36*P)[0]
1076##           54202648602164057575419038802/15402543997324146892198790401
1077##           sage: E.multiple_x_numerator(36, x)
1078##           54202648602164057575419038802
1079##           sage: E.multiple_x_denominator(36, x)
1080##           15402543997324146892198790401
1081
1082##         An example where cancellation occurs:
1083##           sage: E = EllipticCurve("88a1")
1084##           sage: P = E([2,2])   # fixed choice of generator
1085##           sage: n = E.multiple_x_numerator(11, P[0]); n
1086##           442446784738847563128068650529343492278651453440
1087##           sage: d = E.multiple_x_denominator(11, P[0]); d
1088##           1427247692705959881058285969449495136382746624
1089##           sage: n/d
1090##           310
1091##           sage: 11*P
1092##           (310 : -5458 : 1)
1093         
1094##         """
1095##         if cache is None:
1096##             cache = {}
1097
1098##         if x is None:
1099##             x = rings.PolynomialRing(self.base_ring(), 'x').gen()
1100
1101##         n = int(n)
1102##         if n < 2:
1103##             print "n must be at least 2"
1104
1105##         self.pseudo_torsion_polynomial( -2, x, cache)
1106##         self.pseudo_torsion_polynomial(n-1, x, cache)
1107##         self.pseudo_torsion_polynomial(n  , x, cache)
1108##         self.pseudo_torsion_polynomial(n+1, x, cache)
1109
1110##         if n % 2 == 0:
1111##             return x * cache[-1] * cache[n]**2 - cache[n-1] * cache[n+1]
1112##         else:
1113##             return x * cache[n]**2 - cache[-1] * cache[n-1] * cache[n+1]
1114
1115
1116##     def multiple_x_denominator(self, n, x=None, cache=None):
1117##         r"""
1118##         Returns the denominator of the x-coordinate of the nth multiple of
1119##         a point, using torsion polynomials (division polynomials).
1120
1121##         The inputs n, x, cache are as described in pseudo_torsion_polynomial().
1122
1123##         The result is adjusted to be correct for both even and odd n.
1124
1125##         SEE ALSO:
1126##           -- multiple_x_numerator()
1127
1128##         TODO: the numerator and denominator versions share a calculation,
1129##         namely squaring $\psi_n$. Maybe would be good to offer a combined
1130##         version to make this more efficient.
1131
1132##         EXAMPLES:
1133##            -- see multiple_x_numerator()
1134
1135##         AUTHORS:
1136##            -- David Harvey (2006-09-24)
1137
1138##         """
1139##         if cache is None:
1140##             cache = {}
1141
1142##         if x is None:
1143##             x = rings.PolynomialRing(self.base_ring(), 'x').gen()
1144
1145##         n = int(n)
1146##         if n < 2:
1147##             print "n must be at least 2"
1148
1149##         self.pseudo_torsion_polynomial(-2, x, cache)
1150##         self.pseudo_torsion_polynomial(n , x, cache)
1151
1152##         if n % 2 == 0:
1153##             return cache[-1] * cache[n]**2
1154##         else:
1155##             return cache[n]**2
1156
1157
1158    def torsion_polynomial(self, n, var='x', i=0):
1159        """
1160        Returns the n-th torsion polynomial (a.k.a., division polynomial).
1161
1162        INPUT:
1163            n -- non-negative integer
1164            var -- string; the indeterminate of the polynomial (default: x)
1165            i -- integer, either 0 (default) or 1.
1166           
1167        OUTPUT:
1168            Polynomial -- n-th torsion polynomial, which is a polynomial over
1169                          the base field of the elliptic curve.
1170
1171        SEE ALSO:
1172            pseudo_torsion_polynomial()
1173
1174        EXAMPLES:
1175            sage: E = EllipticCurve([0,0,1,-1,0])
1176            sage: E.division_polynomial(1)
1177            1
1178            sage: E.division_polynomial(2)
1179            4*x^3 - 4*x + 1
1180            sage: E.division_polynomial(3, 'z')
1181            3*z^4 - 6*z^2 + 3*z - 1
1182
1183            sage: E = EllipticCurve([0, -1, 1, -10, -20])
1184            sage: E.torsion_polynomial(0)
1185            0
1186            sage: E.torsion_polynomial(1)
1187            1
1188            sage: E.torsion_polynomial(2)
1189            4*x^3 - 4*x^2 - 40*x - 79
1190            sage: E.torsion_polynomial(3)
1191            3*x^4 - 4*x^3 - 60*x^2 - 237*x - 21
1192            sage: E.torsion_polynomial(4)
1193            8*x^9 - 24*x^8 - 464*x^7 - 2758*x^6 + 6636*x^5 + 34356*x^4 + 53510*x^3 + 99714*x^2 + 351024*x + 459859
1194
1195            sage: E = EllipticCurve([-4,0])
1196            sage: E.torsion_polynomial(2)
1197            4*x^3 - 16*x
1198            sage: E.torsion_polynomial(5)
1199            5*x^12 - 248*x^10 - 1680*x^8 + 19200*x^6 - 32000*x^4 + 51200*x^2 + 4096
1200            sage: E.torsion_polynomial(6)
1201            12*x^19 - 1200*x^17 - 18688*x^15 + 422912*x^13 - 2283520*x^11 + 9134080*x^9 - 27066368*x^7 + 19136512*x^5 + 19660800*x^3 - 3145728*x
1202
1203        AUTHOR: David Kohel (kohel@maths.usyd.edu.au), 2005-04-25
1204        """
1205        n = int(n)
1206        try:
1207            return self.__torsion_polynomial[n]
1208        except AttributeError:
1209            self.__torsion_polynomial = {}
1210        except KeyError:
1211            pass
1212        E = self; i=int(i)
1213        if n < 0:
1214            raise ValueError, "n must be a non-negative integer."
1215        if i > 1 : 
1216            raise ValueError, "i must be 0 or 1."
1217       
1218        R = rings.PolynomialRing(E.base_ring(), var)
1219        if i == 1:
1220            if n == 0:
1221                f = E.torsion_polynomial(1)
1222                E.__torsion_polynomial[n] = f
1223                return f
1224            else: 
1225                x = R.gen()
1226                psi2 = E.torsion_polynomial(2)
1227                if n%2 == 0: 
1228                    f = x * psi2 * (E.torsion_polynomial(n)//psi2)**2 - \
1229                        E.torsion_polynomial(n+1) * E.torsion_polynomial(n-1)
1230                    E.__torsion_polynomial[n] = f
1231                    return f
1232                else:
1233                    f = x * E.torsion_polynomial(n)**2 - \
1234                        (E.torsion_polynomial(n+1)//psi2) * E.torsion_polynomial(n-1)
1235                    E.__torsion_polynomial[n] = f
1236                    return f
1237               
1238        else:
1239           
1240            if n == 0: return R(0)
1241            if n == 1: return R(1)
1242            x = R.gen()
1243            b2, b4, b6, b8 = E.b_invariants()
1244            psi2 = 4*x**3 + b2*x**2 + 2*b4*x + b6
1245            if n == 2: 
1246                f = psi2
1247                E.__torsion_polynomial[n] = f; return f
1248            if n == 3:
1249                f = 3*x**4 + b2*x**3 + 3*b4*x**2 + 3*b6*x + b8
1250                E.__torsion_polynomial[n] = f; return f
1251            if n == 4:
1252                f = psi2 * (2*x**6 + b2*x**5 + 5*b4*x**4 + 10*b6*x**3 \
1253                    + 10*b8*x**2  + (b2*b8 - b4*b6)*x + b4*b8 - b6**2)
1254                E.__torsion_polynomial[n] = f; return f
1255            if n%2 == 0:
1256                m = n//2 
1257                if m%2 == 0:
1258                    f = E.torsion_polynomial(m) * ( \
1259                        (E.torsion_polynomial(m+2)//psi2) * E.torsion_polynomial(m-1)**2 - \
1260                        (E.torsion_polynomial(m-2)//psi2) * E.torsion_polynomial(m+1)**2)
1261                    E.__torsion_polynomial[n] = f; return f
1262                else:
1263                    f = psi2 * E.torsion_polynomial(m)*( \
1264                        E.torsion_polynomial(m+2) * (E.torsion_polynomial(m-1)//psi2)**2 - \
1265                        E.torsion_polynomial(m-2) * (E.torsion_polynomial(m+1)//psi2)**2)
1266                    E.__torsion_polynomial[n] = f; return f                   
1267            else: 
1268                m = n//2 
1269                if m%2 == 0:
1270                    f = psi2 * \
1271                        E.torsion_polynomial(m+2) * (E.torsion_polynomial(m)//psi2)**3 - \
1272                        E.torsion_polynomial(m-1) * E.torsion_polynomial(m+1)**3
1273                    E.__torsion_polynomial[n] = f; return f                   
1274                else:
1275                    f = E.torsion_polynomial(m+2) * E.torsion_polynomial(m)**3 - psi2 * \
1276                        E.torsion_polynomial(m-1)*(E.torsion_polynomial(m+1)//psi2)**3
1277                    E.__torsion_polynomial[n] = f; return f
1278
1279    division_polynomial = torsion_polynomial
1280           
1281    def weierstrass_model(self):
1282        """
1283        Return a model of the form $y^2 = x^3 + a*x + b$ for this curve.
1284
1285        EXAMPLES:
1286            sage: E = EllipticCurve('37a')
1287            sage: print E
1288            Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field
1289            sage: F = E.weierstrass_model()
1290            sage: print F
1291            Elliptic Curve defined by y^2  = x^3 - x + 1/4 over Rational Field
1292            sage: print F.minimal_model() == E.minimal_model()
1293            True
1294
1295            sage: E = EllipticCurve([1,2,3,4,5])
1296            sage: F = E.weierstrass_model()
1297            sage: print F
1298            Elliptic Curve defined by y^2  = x^3 + 61/16*x + 127/32 over Rational Field
1299            sage: print F.minimal_model() == E.minimal_model()
1300            True       
1301        """
1302        import constructor
1303        c4, c6 = self.c_invariants()
1304        return constructor.EllipticCurve([-c4/(2**4*3), -c6/(2**5*3**3)])
1305       
1306   
1307
1308    ##############################################################################
1309    # Plotting
1310    ##############################################################################
1311
1312    def plot(self, xmin=None, xmax=None, **args):
1313        """
1314        Draw a graph of this elliptic curve.
1315
1316        INPUT:
1317            xmin, xmax -- points will be computed at least within this
1318                          rings, but possibly farther.  These may be
1319                          left off.
1320            **args -- all other options are passed to the line graphing
1321                      primitive.
1322
1323        EXAMPLES:
1324            sage: E = EllipticCurve([0,-1])
1325            sage: plot(E, rgbcolor=hue(0.7))
1326            Graphics object consisting of 1 graphics primitive
1327        """
1328        RR = rings.RealField()
1329        K = self.base_ring()
1330        try:
1331            RR._coerce_(K(1))
1332        except TypeError:
1333            raise NotImplementedError, "Plotting of curves over %s not implemented yet"%K
1334        a1, a2, a3, a4, a6 = self.ainvs()
1335        R = rings.PolynomialRing(rings.RealField(), 'x')
1336        x = R.gen()
1337        d = 4*x**3 + (a1**2 + 4*a2)*x**2 + (2*a3*a1 + 4*a4)*x + (a3**2 + 4*a6)
1338        def f1(z):
1339            return (-(a1*z + a3) + sqrt(abs(d(z))))/2
1340        def f2(z):
1341            return (-(a1*z + a3) - sqrt(abs(d(z))))/2
1342        r = d.roots(multiplicities=False)
1343        r.sort()
1344        if xmax is None:
1345            xmax = r[-1] + 2
1346        xmax = max(xmax, r[-1]+2)
1347        if xmin is None:
1348            xmin = r[0]  - 2
1349        xmin = min(xmin, r[0]-2)
1350        if len(r) == 1:
1351            # one real root; 1 component
1352            I = [(r[0],xmax)]
1353        else:
1354            # three real roots; 2 components
1355            I = [(r[0],r[1]), (r[2],xmax)]
1356        I = [(max(a,xmin),min(b,xmax)) for a,b in I]
1357
1358        g = plot.Graphics()
1359        try:
1360            plot_points = int(args['plot_points'])
1361            del args['plot_points']
1362        except KeyError:
1363            plot_points = 100
1364           
1365        for j in range(len(I)):
1366            a,b = I[j]
1367            delta = (b-a)/float(plot_points)
1368            v = []
1369            w = []
1370            for i in range(plot_points):
1371                x = a + delta*i
1372                v.append((x, f1(x)))
1373                w.append((x, f2(x)))
1374            v.append((b,f1(b)))
1375            w.append((b,f2(b)))
1376            if len(I) == 2 and j == 0:  # two components -- the oh.
1377                g += plot.line(v + list(reversed(w)) + [v[0]], **args)
1378            else:
1379                g += plot.line(list(reversed(v)) + w, **args)
1380        return g
1381
1382    def formal_group(self):
1383        r"""
1384        The formal group associated to this elliptic curve.
1385
1386        EXAMPLES:
1387            sage: E = EllipticCurve("37a")
1388            sage: E.formal_group()
1389            Formal Group associated to the Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field           
1390        """
1391        try:
1392            return self.__formal_group
1393        except AttributeError:
1394            self.__formal_group = formal_group.EllipticCurveFormalGroup(self)
1395            return self.__formal_group
1396
1397    formal = formal_group
1398
1399
1400    def hyperelliptic_polynomials(self):
1401        K = self.base_ring()
1402        R = PolynomialRing(K, 'x')
1403        x = R.gen(0)
1404        a1, a2, a3, a4, a6 = self.ainvs()
1405        return R([a6, a4, a2, 1]), R([a3, a1])
Note: See TracBrowser for help on using the repository browser.