Changeset 7681:04ba6e30689b

Timestamp:

11/14/07 07:13:48 (6 years ago)

Author:

Robert Bradshaw <robertwb@…>

Branch:

default

Message:

Integrate new simon two descent code

Location:

sage

Files:

• categories/morphism.pyx (modified) (1 diff)
• schemes/elliptic_curves/ell_generic.py (modified) (2 diffs)
• schemes/elliptic_curves/ell_number_field.py (modified) (1 diff)
• schemes/elliptic_curves/ell_rational_field.py (modified) (1 diff)
• schemes/elliptic_curves/gp_simon.py (modified) (3 diffs)

sage/categories/morphism.pyx

@@ -213 +213 @@
     def __pow__(self, n, dummy):
         return self
+        
+    def __invert__(self):
+        return self
             
 

sage/schemes/elliptic_curves/ell_generic.py

@@ -44 +44 @@
 import sage.modular.modform as modform
 import sage.functions.transcendental as transcendental
+
+from sage.categories.morphism import IdentityMorphism
+from sage.categories.homset import Hom
+from sage.rings.arith import lcm
 
 # Schemes
@@ -1304 +1308 @@
         return constructor.EllipticCurve([-c4/(2**4*3), -c6/(2**5*3**3)])
         
+    def integral_model(self):
+        denom = lcm([a.denominator() for a in self.ainvs()])
+        if denom != 1:
+            raise NotImplementedError, "model must be integral for now"
+        else:
+            parent = self(0).parent()
+            return self, IdentityMorphism(Hom(parent, parent))
+
     
 

sage/schemes/elliptic_curves/ell_number_field.py

@@ -132 +132 @@
         #if self.torsion_order() % 2 == 0:
         #    raise ArithmeticError, "curve must not have rational 2-torsion\nThe *only* reason for this is that I haven't finished implementing the wrapper\nin this case.  It wouldn't be too difficult.\nPerhaps you could do it?!  Email me (wstein@gmail.com)."
-        F = self.integral_weierstrass_model()
-        a1,a2,a3,a4,a6 = F.a_invariants()
+        # F = self.integral_weierstrass_model()
+        # a1,a2,a3,a4,a6 = F.a_invariants()
         x = PolynomialRing(self.base_ring(), 'x').gen(0)
-        t = simon_two_descent(a2,a4,a6, 
+        t = simon_two_descent(self,
                               verbose=verbose, lim1=lim1, lim3=lim3, limtriv=limtriv,
                               maxprob=maxprob, limbigprime=limbigprime)
         prob_rank = Integer(t[0])
         two_selmer_rank = Integer(t[1])
-        prob_gens = [F(P) for P in t[2]]
+        prob_gens = [self(P) for P in t[2]]
         return prob_rank, two_selmer_rank, prob_gens
                 

sage/schemes/elliptic_curves/ell_rational_field.py

@@ -858 +858 @@
             (8, 8)
         """
-        if self.torsion_order() % 2 == 0:
-            raise ArithmeticError, "curve must not have rational 2-torsion\nThe *only* reason for this is that I haven't finished implementing the wrapper\nin this case.  It wouldn't be too difficult.\nPerhaps you could do it?!  Email me (wstein@gmail.com)."
-        F = self.integral_weierstrass_model()
-        a1,a2,a3,a4,a6 = F.a_invariants()
-        t = simon_two_descent(a2,a4,a6, verbose=verbose, lim1=lim1, lim3=lim3, limtriv=limtriv,
+        t = simon_two_descent(self, verbose=verbose, lim1=lim1, lim3=lim3, limtriv=limtriv,
                               maxprob=maxprob, limbigprime=limbigprime)
         prob_rank = rings.Integer(t[0])
         two_selmer_rank = rings.Integer(t[1])
-        prob_gens = [F(P) for P in t[2]]
+        prob_gens = [self(P) for P in t[2]]
         return prob_rank, two_selmer_rank, prob_gens
 

sage/schemes/elliptic_curves/gp_simon.py

@@ -28 +28 @@
 
 gp = None
-def init(K, two_torsion):
+def init():
     global gp
     if gp is None:
         gp = Gp(script_subdirectory='simon')
         gp.read("ell.gp")
-        if K == QQ and two_torsion == 1:
-            gp.read("ellQ.gp")
+        gp.read("ellQ.gp")
         gp.read("qfsolve.gp")
         gp.read("resultant3.gp")
 
 
-def simon_two_descent(A, B, C, verbose=0, lim1=5, lim3=50, limtriv=10, maxprob=20, limbigprime=30):
+def simon_two_descent(E, verbose=0, lim1=5, lim3=50, limtriv=10, maxprob=20, limbigprime=30):
 
-    K = A.parent()
-    x = PolynomialRing(K, 'x').gen(0)
-    f = x**3 + A*x**2 + B*x + C
-    factors = f.factor()
-    two_torsion = 1 << (len(factors) - 1)
+    init()
+    K = E.base_ring()
+    F, transform = E.integral_model()
     
-    init(K, two_torsion)
-    shift = 0
-    
-    if K == QQ:
-        gp.eval("K = bnfinit(y-1)") # TODO: Is there a better way to get QQ?
-    else:
+    if K != QQ:
         # Simon's program requires that this name be y. 
         with localvars(K.polynomial().parent(), 'y'):
@@ -62 +54 @@
             print "%s = Mod(y,K.pol);" % K.gen()
     
-    if two_torsion == 1:
-        if K == QQ:
-            cmd = 'main(%s,%s,%s);'%(A,B,C)
-        else:
-            # Simon's program requires that the coefficients do NOT all lie in a subfield of K
-            if A.minpoly().degree() != K.degree():
-                shift = K.gen()
-                shifted = f(x+shift)
-                C, B, A, _ = shifted
-            cmd = 'main(K, %s,%s,%s);'%(A,B,C)
+    if K == QQ:
+        cmd = 'ellrank([%s,%s,%s,%s,%s]);' % tuple(F.ainvs())
+    else:
+        cmd = 'bnfellrank(K, [%s,%s,%s,%s,%s]);' % tuple(F.ainvs())
             
-    elif two_torsion == 2:
-        if C != 0:
-            # If E[2] = Z/2Z, it must be of the form y^2 = x^3 + Ax^2 + Bx
-            shift = -factors[0][0][0]
-            shifted = f(x+shift)
-            C, B, A, _ = shifted
-        cmd = 'main2(K, [0, %s, 0, %s, 0])'%(A,B)
-
-    else:  # two_torsion == 4:
-        cmd = 'complete(K, %s,%s,%s)'%(factors[0][0][0], factors[1][0][0], factors[2][0][0])
-
     gp('DEBUGLEVEL=%s; LIM1=%s; LIM3=%s; LIMTRIV=%s; MAXPROB=%s; LIMBIGPRIME=%s;'%(
         verbose, lim1, lim3, limtriv, maxprob, limbigprime))
+
     if verbose >= 2:
         print cmd
@@ -101 +77 @@
         return args[0]
     ans = sage_eval(v, {'Mod': gp_mod, 'y': K.gen(0)})
-    if two_torsion != 1:
-        # selmer group rank and generating points not
-        # returned for non-trivial 2-torsion
-        ans = [ans, -1, []]
-    if shift:
-        # Undo the x -> x+a translation above
-        ans[2] = [[P[0]+shift, P[1]] for P in ans[2]]
+    inv_transform = ~transform
+    ans[2] = [inv_transform(F(P)) for P in ans[2]]
     return ans