Ticket #3674: sage-trac3674new-extra2.patch

File sage-trac3674new-extra2.patch, 7.0 KB (added by cremona, 9 years ago)

Apply this after the preceding two to 3.0.6 or 3.1.rc0

  • sage/schemes/elliptic_curves/ell_point.py

    # HG changeset patch
    # User John Cremona <john.cremona@gmail.com>
    # Date 1218213448 -3600
    # Node ID b4bdd819f08911aec9ba47ab462c7a9bcc39b1c5
    # Parent  d7f85e1920ee445b0e9643a14aa8974aa1d21dde
    #3674: some bug fixes + more efficient version of integral_points_with_bounded_mw_ceoffs()
    
    diff -r d7f85e1920ee -r b4bdd819f089 sage/schemes/elliptic_curves/ell_point.py
    a b class EllipticCurvePoint_field(AdditiveG 
    818818        if disc < 0: #Connected Case
    819819            # Here we use formulae equivalent to those in Cohen, but better
    820820            # behaved when roots are close together
     821
     822            # For some curves (e.g. 3314b3) the default precision is not enough!
     823            while len(real_roots)!=1:
     824                precision*=2
     825                RR=rings.RealField(precision)
     826                real_roots = pol.roots(RR,multiplicities=False)
    821827            try:
    822828                assert len(real_roots) == 1
    823829            except:
    824830                raise ValueError, ' no or more than one real root despite disc < 0'
    825831            e1 = real_roots[0]
    826             roots = pol.roots(CC,multiplicities=False)
     832            roots = pol.roots(rings.ComplexField(precision),multiplicities=False)
    827833            roots.remove(e1)
    828834            e2,e3 = roots
    829835            zz = (e1-e2).sqrt() # complex
    class EllipticCurvePoint_field(AdditiveG 
    843849            return z
    844850
    845851        else:                    #Disconnected Case, disc > 0
     852            # For some curves (e.g. 2370i5) the default precision is not enough!
     853            while len(real_roots)!=3:
     854                precision*=2
     855                RR=rings.RealField(precision)
     856                real_roots = pol.roots(RR,multiplicities=False)
     857
    846858            real_roots.sort() # increasing order
    847859            real_roots.reverse() # decreasing order e1>e2>e3
    848860            try:
  • sage/schemes/elliptic_curves/ell_rational_field.py

    diff -r d7f85e1920ee -r b4bdd819f089 sage/schemes/elliptic_curves/ell_rational_field.py
    a b class EllipticCurve_rational_field(Ellip 
    15251525            points = self.gens()
    15261526        else:
    15271527            for P in points:
    1528                 assert P.curve() is self
     1528                assert P.curve() == self
    15291529
    15301530        r = len(points)
    15311531        if precision is None:
    class EllipticCurve_rational_field(Ellip 
    38863886            bounded by $H_q$.  The bound $H_q$ will be computed at
    38873887            runtime.
    38883888            """
     3889            from sage.misc.all import cartesian_product_iterator
     3890            from sage.groups.generic import multiples
    38893891            xs=set()
    3890             for i in range(ceil(((2*H_q+1)**r)/2)):
    3891                 koeffs = Z(i).digits(base=2*H_q+1)
    3892                 koeffs = [0]*(r-len(koeffs)) + koeffs # pad with 0s
    3893                 P = sum([(koeffs[j]-H_q)*mw_base[j] for j in range(r)],self(0))
    3894                 for Q in tors_points: # P + torsion points (includes 0)
    3895                     tmp = P + Q
    3896                     if not tmp.is_zero():
    3897                         x = tmp[0]
    3898                         if x.is_integral():
    3899                             xs.add(x)
     3892            N=H_q
     3893           
     3894            def use(P):
     3895                """
     3896                Helper function to record x-coord of a point if integral.
     3897                """
     3898                if not P.is_zero():
     3899                    xP = P[0]
     3900                    if xP.is_integral():
     3901                        xs.add(xP)
     3902                       
     3903            def use_t(R):
     3904                """
     3905                Helper function to record x-coords of a point +torsion if integral.
     3906                """
     3907                for T in tors_points:
     3908                    use(R+T)
     3909                       
     3910            # We use a naive method when the number of possibilities is small:
     3911
     3912            if r==1 or N<=10:
     3913                for P in multiples(mw_base[0],N+1):
     3914                    use_t(P)
     3915                return xs
     3916           
     3917            # Otherwise is is very very much faster to first compute
     3918            # the linear combinations over RR, and only compte them as
     3919            # rational points if they are approximately integral
     3920
     3921            def is_approx_integral(P):
     3922                return (abs(P[0]-P[0].round()))<0.001 and (abs(P[1]-P[1].round()))<0.001
     3923           
     3924            RR = RealField() #(100)
     3925            ER = self.change_ring(RR)
     3926            Rgens = [ER.lift_x(P[0]) for P in mw_base]
     3927            for i in range(r):
     3928                if abs(Rgens[i][1]-mw_base[i][1])>abs((-Rgens[i])[1]-mw_base[i][1]):
     3929                    Rgens[i] = -Rgens[i]
     3930            for ni in cartesian_product_iterator([range(-N,N+1) for i in range(r-1)]+[range(N+1)]):
     3931                RP=sum([ni[i]*Rgens[i] for i in range(r)],ER(0))
     3932                for T in tors_points:
     3933                    if is_approx_integral(RP+ER(T)):
     3934                        P = sum([ni[i]*mw_base[i] for i in range(r)],T)
     3935                        use(P)
    39003936            return xs
     3937           
     3938            # Below this point we keep earlier experimental code which
     3939            # is slower when either N or r is large
     3940           
     3941#           if r==2:
     3942#               for P in multiples(mw_base[0],N+1):
     3943#                   for Q in multiples(mw_base[1],N+1):
     3944#                       for R in set([Q+P,Q-P]):
     3945#                           use_t(R)
     3946#               return xs
     3947#           
     3948#           if r==3:
     3949#               for P in multiples(mw_base[0],N+1):
     3950#                   for Q in multiples(mw_base[1],N+1):
     3951#                       PQ = [P+Q,P-Q]
     3952#                       PQ = set(PQ + [-R for R in PQ]) # {\pm P\pm Q}
     3953#                       for R in multiples(mw_base[2],N+1):
     3954#                           for S in PQ:
     3955#                               use_t(R+S)
     3956#               return xs
     3957#           
     3958#           # general rank
     3959#           mults=[list(multiples(mw_base[i],N+1)) for i in range(r)]
     3960#           for i in range(r-1):
     3961#               mults[i] = [-P for P in mults[i] if not P.is_zero()] + mults[i]
     3962#           for Pi in cartesian_product_iterator(mults):
     3963#               use_t(sum(Pi,self(0)))
     3964#           return xs
     3965#
     3966#  older, even slower  code:
     3967#
     3968#           for i in range(ceil(((2*H_q+1)**r)/2)):
     3969#               koeffs = Z(i).digits(base=2*H_q+1)
     3970#               koeffs = [0]*(r-len(koeffs)) + koeffs # pad with 0s
     3971#               P = sum([(koeffs[j]-H_q)*mw_base[j] for j in range(r)],self(0))
     3972#               for Q in tors_points: # P + torsion points (includes 0)
     3973#                   tmp = P + Q
     3974#                   if not tmp.is_zero():
     3975#                       x = tmp[0]
     3976#                       if x.is_integral():
     3977#                           xs.add(x)
     3978#           return xs
    39013979        ##############################  end  #################################
    39023980   
    39033981        # END Internal functions #############################################