source: sage/schemes/elliptic_curves/ell_generic.py @ 6528:0c4c7907d97f

Revision 6528:0c4c7907d97f, 47.8 KB checked in by William Stein <wstein@…>, 6 years ago (diff)

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