attachment:intpts.sage of days22/cremona/to-do-list

   1 #        File: intpts.sage
   2 #     Created: Thu Jul 01 04:22 PM 2010 C
   3 # Last Change: Thu Jul 01 07:40 PM 2010
   4 # Original Magma Code: Thotsaphon "Nook" Thongjunthug 
   5 # Sage version: Radoslav Kirov, Jackie Anderson
   6 
   7 from sage.misc.all import verbose
   8 
   9 # This function should be detached and included as part of projective points over number field
  10 def abs_log_height(X, gcd_one = True, precision = None):
  11     r''' Computes the height of a point in a projective space over field K.
  12     It assumes the coordinates have gcd equal to 1.
  13     If not use gcd_one flag.
  14     '''
  15     assert isinstance(X,list)
  16     K = X[0].parent()
  17     if precision is None:
  18         RR = RealField()
  19         CC = ComplexField()
  20     else:
  21         RR = RealField(precision)
  22         CC = ComplexField(precision)
  23     places = set([])
  24     if K == QQ:
  25         embs = K.embeddings(RR)
  26         Xideal = X
  27     else:
  28         embs = K.places(precision)
  29         # skipping zero as it currently over K breaks Sage
  30         Xideal = [K.ideal(x) for x in X if not x.is_zero()]
  31     for x in Xideal:
  32         places = places.union(x.denominator().prime_factors())
  33         if not gcd_one:
  34             places = places.union(x.numerator().prime_factors())
  35     if K == QQ:
  36         non_arch = sum([log(max([RR(P)**(-x.valuation(P)) for x in X])) for P in places])
  37     else:
  38         non_arch = sum([P.residue_class_degree() *
  39                         P.absolute_ramification_index() *
  40                         max([x.abs_non_arch(P, precision) for x in X]).log() for P in places])
  41     arch = 0
  42     r,s = K.signature()
  43     for i,f in enumerate(embs):
  44         if i<r:
  45             arch += max([f(x).abs() for x in X]).log()
  46         else:
  47             arch += max([f(x).abs()**2 for x in X]).log()
  48     return (arch+non_arch)/K.degree()
  49     
  50 def compute_c6(E,emb):
  51     x = var('x')
  52     #f = x**3-27*emb(E.c4())*x-54*emb(E.c6())
  53     f = x**3-emb(E.c4()/48)*x-emb(E.c6()/864)
  54     R = f.roots(multiplicities = False)
  55     m = max([CC(r).abs() for r in R])
  56     return 2*m
  57 
  58 def compute_c8(L):
  59     # Original code transformed the lattice generators.
  60     # Here we assume they are of standard form.
  61     w1, w2 = L
  62     m = max(CC(w1).abs(), CC(w2).abs(), CC(w1+w2).abs())
  63     return m 
  64 
  65 # There is a Sage trak ticket implementing this. In the future the library function can be used and this one removed 
  66 def height_pairing_matrix(points, precision = None):
  67     r = len(points)
  68     if precision is None:
  69         RR = RealField()
  70     else:
  71         RR = RealField(precision)
  72     M = MatrixSpace(RR, r)
  73     mat = M()
  74     for j in range(r):
  75         mat[j,j] = points[j].height(precision = precision)
  76     for j in range(r):
  77         for k in range(j+1,r):
  78             mat[j,k] = ((points[j]+points[k]).height(precision=precision) - mat[j,j] - mat[k,k])/2
  79             mat[k,j] = mat[j,k]
  80     return mat
  81 
  82 def c3(L):
  83     return min(map(abs,height_pairing_matrix(L).eigenvalues()))
  84 
  85 def h_E(E):
  86     K = E.base_field()
  87     j = E.j_invariant()
  88     c4, c6 = E.c_invariants()
  89     g2, g3 = c4/12, c6/216
  90     return max(abs_log_height([K(1), g2, g3]), abs_log_height([K(1), j]))
  91 
  92 def h_m(E, P, ElogEmbedP, D7):
  93     K = E.base_field()
  94     return max([P.height(), h_E(E), D7/K.degree()*abs(log(ElogEmbedP))**2])
  95     
  96 def Extra_h_m(E, Periods, D7):
  97     D = E.base_field().degree()
  98     h = h_E(E)
  99     return map(lambda x: max([0, h, D7/D*abs(x)**2]), Periods)
 100 
 101 def d8(E, L, Elog, Periods, D7):
 102     C = [exp(1) * h_E(E)]
 103     D = E.base_field().degree()
 104     for i in range(len(L)):
 105         C.append(h_m(E, L[i], Elog[i], D7) / D)
 106     C += [t / D for t in Extra_h_m(E, Periods, D7)]
 107     return max(C);
 108 
 109 def d9(E, L, Elog, Periods, D7):
 110     D = E.base_field().degree()
 111     C = []
 112     for i in range(len(L)):
 113         tmp = exp(1) * sqrt(D * h_m(E, L[i], Elog[i], D7) / D7)
 114         C.append( tmp / abs(Elog[i]))
 115     Ehm = Extra_h_m(E, Periods, D7)
 116     C += [exp(1) * sqrt(D*Ehm[i]/D7) / abs(Periods[i]) for i in [0,1]]
 117     return min(C);
 118 
 119 def d10(E, L, Elog, Periods, D7):
 120     D = E.base_field().degree()
 121     n = len(L)+2
 122     scalar = 2 * 10**(8 + 7*n) * (2/exp(1))**(2*n**2)
 123     scalar *= (n+1)**(4*n**2 + 10*n) * D**(2*n + 2)
 124     scalar *= (log(d9(E, L, Elog, Periods, D7)))**(-2*n-1)
 125     for i in range(len(L)):
 126         scalar *= h_m(E, L[i], Elog[i], D7)
 127     scalar *= prod(Extra_h_m(E, Periods, D7))   
 128     return scalar
 129 
 130 def RHS(D, r, C9, C10, D9, D10, h, Q, expTors):
 131     bound = (log(log(Q*r*expTors)) + h + log(D*D9))**(r+3)
 132     bound *= D10*(log(Q*r*expTors) + log(D*D9))
 133     bound += log(C9*expTors)
 134     bound /= C10
 135     return bound
 136 
 137 def InitialQ(D, r, Q0, C9, C10, D8, D9, D10, h, expTors):
 138     minQ = max(Q0, exp(D8))
 139     Q = minQ + 1
 140     x = ceil(log(Q, 10)) # x = log_10(Q)
 141     exp_OK = 0   # the exponent that satisfies P.I.
 142     exp_fail = 0 # the exponent that fails P.I.
 143     while 10**(2*x) < RHS(D, r, C9, C10, D9, D10, h, 10**x, expTors):
 144         exp_OK = x # Principal Inequality satisfied
 145         x *= 2     # double x, and retry
 146     exp_fail = x # x that fails the Principal Inequality
 147     
 148     # So now x = log_10(Q) must lie between exp_OK and exp_fail
 149     # Refine x further using "binary search"
 150     while True:
 151         x = floor((exp_OK + exp_fail)/2)
 152         if 10**(2*x) >= RHS(D, r, C9, C10, D9, D10, h, 10**x, expTors): 
 153             exp_fail = x
 154         else:
 155             exp_OK = x
 156         if (exp_fail - exp_OK) <= 1:
 157             break
 158     return exp_fail # over-estimate
 159 
 160 def Faltings_height(E):
 161     K = E.base_field()
 162     c = log(2)
 163     if E.b2().is_zero():
 164         c = 0
 165     h1 = abs_log_height([K(E.discriminant()), K(1)])/6
 166     h2 = K(E.j_invariant()).global_height_arch()/6
 167     h3 = K(E.b2()/12).global_height_arch()
 168     return n(h1 + h2/2 + h3/2 + c)
 169 
 170 def Silverman_height_bounds(E):
 171     K = E.base_field()
 172     mu = Faltings_height(E)
 173     lb = -mu-2.14
 174     ub = abs_log_height([K(E.j_invariant()), K(1)])/12 + mu + 1.922
 175     return lb, ub
 176 
 177 def ReducedQ(E, f, LGen, Elog, C9, C10, Periods, expTors, initQ):
 178     w1, w2 = Periods 
 179     r = len(LGen)
 180     newQ = initQ
 181     # Repeat LLL reduction until no further reduction is possible
 182     while True: 
 183         Q = newQ
 184         S = r*(Q**2)*(expTors**2)
 185         T = 3*r*Q*expTors/sqrt(2)
 186         # Create the basis matrix
 187         C = 1
 188         while True: 
 189             C *= Q**ceil((r+2)/2)    
 190             precbits = int(log(C,2)+10)
 191             L = copy(zero_matrix(ZZ, r+2))
 192             # Elliptic logarithm should have precision "suitable to" C
 193             # e.g. If C = 10**100, then Re(Elog[i]) should be computed
 194             # correctly to at least 100 decimal places
 195             if precbits > Elog[0].prec(): 
 196                 verbose( "Need higher precision, recompute elliptic logarithm ...")
 197                 # Re-compute elliptic logarithm to the right precision
 198                 verbose( "precision in bits %i" % precbits)
 199                 Elog = [ P.elliptic_logarithm(f, precision = precbits) for P in LGen]
 200                 verbose( "Elliptic logarithm recomputed")
 201                 w1, w2 = E.period_lattice(f).normalised_basis( prec = precbits)
 202             # Assign all non-zero entries
 203             for i in range(r): 
 204                 L[i, i] = 1
 205                 L[r, i]   = (C*Elog[i].real_part()).trunc()
 206                 L[r+1, i] = (C*Elog[i].imag_part()).trunc()
 207             L[r , r ]   = (C*w1.real_part()).trunc()
 208             L[r , r+1 ] = (C*w2.real_part()).trunc()
 209             L[r+1, r]   = (C*w1.imag_part()).trunc()
 210             L[r+1, r+1] = (C*w2.imag_part()).trunc()
 211             # LLL reduction and constants
 212             L = L.transpose()
 213             L = L.LLL()
 214             b1 = L[0] # 1st row of reduced basis
 215             # Norm(b1) = square of Euclidean norm
 216             normb1 = sum([i**2 for i in b1])
 217             lhs = 2**(-r-1)*normb1 - S 
 218             if (lhs >= 0) and (sqrt(lhs) > T):
 219                 break
 220         newQ = ( log(C*C9*expTors) - log(sqrt(lhs) - T) ) / C10
 221         newQ = floor(sqrt(newQ))
 222         verbose( "After reduction, Q <= %f" % newQ)
 223         if ((Q - newQ) <= 1) or (newQ <= 1):
 224             break
 225     return newQ
 226 
 227 
 228 #// Search for all integral points on elliptic curves over number fields
 229 #// within a given bound
 230 #// Input:    E = elliptic curve over a number field K
 231 #//           L = a sequence of points in the Mordell-Weil basis for E(K)
 232 #//           Q = a maximum for the absolute bound on all coefficients
 233 #//               in the linear combination of points in L
 234 #// Output:  S1 = sequence of all integral points on E(K) modulo [-1]
 235 #//          S2 = sequence of tuples representing the points as a
 236 #//               linear combination of points in L
 237 #// Option: tol = tolerance used for checking integrality of points.
 238 #//               (Default = 0 - only exact arithmetic will be performed)
 239 
 240 # cyclic group iterator  
 241 # Returns all elements of the cyclic group, remembering all intermediate steps 
 242 # id - identity element
 243 # gens - generators of the group
 244 # mult - orders of the generators
 245 # both_signs - whether to return all elements or one per each {element, inverse} pair
 246 def cyc_iter(id, gens, mult, both_signs = False):
 247     if len(gens) == 0:
 248         yield id,[] 
 249     else:
 250         P = gens[0]
 251         cur = id
 252         if both_signs:
 253             ran = xrange(mult[0])
 254         else:
 255             ran = xrange(mult[0]/2 + 1)
 256         for k in ran: 
 257             for rest, coefs in cyc_iter(id, gens[1:], mult[1:], both_signs or k != 0):
 258                 yield cur + rest, [k] + coefs
 259             cur += P
 260 
 261 #generates elements of form a_1G_1 + ... + a_nG_n
 262 #where |a_i| <= bound and the first non-zero coefficient is positive
 263 def L_points_iter(id, gens, bound, all_zero = True):
 264     if len(gens) == 0:
 265         yield id, []
 266     else:
 267         P = gens[0]
 268         cur = id
 269         for k in xrange(bound+1):
 270             for rest, coefs in L_points_iter(id, gens[1:], bound, all_zero = (all_zero and k == 0)):
 271                 yield cur + rest, [k] + coefs  
 272                 if k!=0 and not all_zero:
 273                    yield -cur + rest, [-k] + coefs 
 274             cur += P
 275 
 276 def integral_points_with_Q(E, L, Q, tol = 0):
 277     r'''Given an elliptic curve E over a number field, its Mordell-Weil basis L, and a positive integer Q, return the sequence of all integral points modulo [-1] of the form P = q_1*L[1] + ... + q_r*L[r] + T with some torsion point T and |q_i| <= Q, followed by a sequence of tuple sequences representing the points as a linear combination of points. An optional tolerance may be given to speed up the computation when checking integrality of points.
 278     '''
 279     assert tol >= 0 
 280     Tors = E.torsion_subgroup()
 281     expTors = Tors.exponent()
 282     OrdG = Tors.invariants() 
 283     Tgens = Tors.gens()
 284     id = E([0,1,0])
 285     total_gens = L + list(Tgens)
 286     verbose( "Generators = %s" % L)
 287     PtsList = []
 288     if tol == 0: 
 289         verbose( "Exact arithmetic")
 290         for P, Pcoeff in L_points_iter(id, L, Q):
 291             for T, Tcoeff in cyc_iter(id, Tgens, OrdG, both_signs = P!=id):
 292                 R = P + T
 293                 if R[0].is_integral() and R[1].is_integral() and R != id:
 294                     Rcoeff = Pcoeff + Tcoeff
 295                     verbose( "%s ---> %s" %(R, Rcoeff))
 296                     PtsList.append(R)
 297         verbose( "*"*45)
 298         return PtsList
 299     # Suggested by John Cremona
 300     # Point search. This is done via arithmetic on complex points on each
 301     # embedding of E. Exact arithmetic will be carried if the resulting
 302     # complex points are "close" to being integral, subject to some tolerance
 303     K = E.base_ring()
 304     r, s = K.signature()
 305     # Set tolerance - This should be larger than 10**(-current precision) to
 306     # avoid missing any integral points. Too large tolerance will slow the
 307     # computation, too small tolerance may lead to missing some integral points.
 308     verbose( "Tolerance = %f" % tol )
 309     if K == QQ:
 310         A = matrix([1])
 311     else:
 312         A = matrix([a.complex_embeddings() for a in K.integral_basis()])
 313         # Note that A is always invertible, so we can take its inverse
 314     B = A.inverse()
 315     point_dict = {}
 316     for emb in K.embeddings(CC):
 317         Ec = EllipticCurve(CC, map(emb,E.a_invariants()))
 318         # need to turn off checking, otherwise the program crashes
 319         Lc = [ Ec.point(map(emb,[P[0],P[1],P[2]]), check = False) for P in L]
 320         Tgensc = [ Ec.point(map(emb,[P[0],P[1],P[2]]), check = False) for P in Tgens]
 321         # Compute P by complex arithmetic for every embedding
 322         id = Ec([0,1,0]) 
 323         for P, Pcoeff in L_points_iter(id, Lc, Q):
 324             for T, Tcoeff in cyc_iter(id, Tgensc, OrdG, both_signs = P!=id):
 325                 R = P + T
 326                 if R == id:
 327                     continue
 328                 key = tuple(Pcoeff + Tcoeff)
 329                 if key in point_dict:
 330                     point_dict[key] += [R]
 331                 else:
 332                     point_dict[key] = [R]
 333     integral_pts = []
 334     false = 0
 335     for Pcoeff in point_dict:
 336         # Check if the x-coordinate of P is "close to" being integral
 337         # If so, compute P exactly and check if it is integral skip P otherwise
 338         XofP = vector([Pemb[1] for Pemb in point_dict[Pcoeff]])
 339         # Write x(P) w.r.t. the integral basis of (K)
 340         # Due to the floating arithmetic, some entries in LX may have very tiny
 341         # imaginary part, which can be thought as zero
 342         LX = XofP * B
 343         try:
 344             LX = [abs( i.real_part() - i.real_part().round() ) for i in LX[0]]
 345         except AttributeError:
 346             LX = [abs(i - i.round()) for i in LX[0]]
 347         if len([1 for dx in LX if dx >= tol]) == 0 :
 348         # x-coordinate of P is not integral, skip P
 349             P = sum([Pcoeff[i]*Pt for i, Pt in enumerate(total_gens)])
 350             if P[0].is_integral() and P[1].is_integral(): 
 351                 integral_pts.append(P)
 352                 verbose( "%s ---> %s" %(P, Pcoeff))
 353             else:
 354                 false += 1
 355     verbose( 'false positives %i'% false)
 356     return integral_pts 
 357 # Compute all integral points on elliptic curve over a number field
 358 # This function simply computes a suitable bound Q, and return
 359 # IntegralPoints(E, L, Q    tol = ...). 
 360 # Input        E = elliptic curve over a number field K
 361 #                     L = a sequence of points in the Mordell-Weil basis for E(K)
 362 # Output    S1 = sequence of all integral points on E(K) modulo [-1]
 363 #                    S2 = sequence of tuples representing the points as a
 364 #                             linear combination of points in L
 365 # Option tol = tolerance used for checking integrality of points.
 366 #                             (Default = 0 - only exact arithmetic will be performed) 
 367 
 368 def calculate_Q(E,L):
 369     if len(L) == 0:
 370         return 0
 371     K = E.base_ring()
 372     expTors = E.torsion_subgroup().exponent()
 373     r, s = K.signature()
 374     pi = RR(math.pi) 
 375     b2 = E.b_invariants()[0] 
 376     # Global constants (independent to the embedding of E)
 377     gl_con = {}
 378     gl_con['c2'] = C2 = - Silverman_height_bounds(E)[0] 
 379     gl_con['c3'] = C3 = c3(L)
 380     gl_con['h_E'] = h = h_E(E)
 381     verbose( "Global constants")
 382     for name, val in gl_con.iteritems():
 383         verbose( '%s = %s'%(name, val))
 384     verbose( "-"*45)
 385     Q = []
 386     # Find the most reduced bound on each embedding of E
 387     # Sage bug, QQ.places() is not implemented
 388     # and K.embeddings gives each complex places twice
 389     if K is QQ:
 390         infplaces = K.embeddings(CC)  
 391     else:
 392         infplaces = K.places()
 393     for i, f in enumerate(infplaces):
 394         # Bug in P.complex_logarithm(QQ.embeddings(CC)[0]) 
 395         if K is QQ:
 396             emb = None
 397         else:
 398             emb = f
 399         if i < r: 
 400             nv = 1
 401             verbose( "Real embedding #%i" % i )
 402         else:
 403             nv = 2
 404             verbose( "Complex embedding #%i" % (i-r))
 405         # Create complex embedding of E
 406         # Embedding of all points in Mordell-Weil basis
 407         # Find complex elliptic logarithm of each embedded point
 408         precbits = floor(100*log(10,2)) 
 409         Elog = [P.elliptic_logarithm(emb, precision = precbits) for P in L]
 410         Periods = E.period_lattice(emb).normalised_basis();
 411         # Local constants (depending on embedding)
 412         # C9, C10 are used for the upper bound on the linear form in logarithm
 413         loc_con = {}
 414         loc_con['c4'] = C4 = C3 * K.degree() / (nv*(r+s))
 415         loc_con['c5'] = C5 = C2 * K.degree() / (nv*(r+s))
 416         loc_con['c6'] = C6 = compute_c6(E,f)
 417         loc_con['delta'] = delta = 1 + (nv-1)*pi
 418         loc_con['c8'] = C8 = compute_c8(Periods)
 419         loc_con['c7'] = C7 = 8*(delta**2) + (C8**2)*abs(f(b2))/12
 420         loc_con['c9'] = C9 = C7*exp(C5/2)
 421         loc_con['c10'] = C10 = C4/2
 422         loc_con['Q0'] = Q0 = sqrt( ( log(C6+abs(f(b2))/12) + C5) / C4 )
 423 
 424         # Constants for David's lower bound on the linear form in logarithm
 425         # Magma and Sage output periods in different order, need to swap w1 and w2
 426         w2, w1 = Periods # N.B. periods are already in "standard" form
 427         loc_con['d7'] = D7 = 3*pi / ((abs(w2)**2) * (w2/w1).imag_part())
 428         loc_con['d8'] = D8 = d8(E, L, Elog, Periods, D7)
 429         loc_con['d9'] = D9 = d9(E, L, Elog, Periods, D7)
 430         loc_con['d10'] = D10 = d10(E, L, Elog, Periods, D7)
 431         for name, val in loc_con.iteritems():
 432             verbose( "{0} = {1}".format(name, val))
 433         # Find the reduced bound for the coefficients in the linear logarithmic form
 434         loginitQ = InitialQ(K.degree(), len(L), Q0, C9, C10, D8, D9, D10, h, expTors)
 435         verbose( "Initial Q <= 10^%f" % loginitQ)
 436         initQ = 10**loginitQ
 437         Q.append(ReducedQ(E, emb, L, Elog, C9, C10, Periods, expTors, initQ))
 438         verbose( "-"*45)
 439     Q = max(Q)
 440     verbose( "Maximum absolute bound on coefficients = %i"% Q)
 441     return Q
 442 
 443 def integral_points(E, L, tol = 0, both_signs = False): 
 444     r'''Given an elliptic curve over a number field and its Mordell-Weil basis, return the sequence of all integral points modulo [-1], followed by a sequence of tuple sequences representing the points as a linear combination of points. An optional tolerance may be given to speed up the computation when checking integrality of points. (This function simply computes a suitable Q and call
 445 IntegralPoints(E, L, Q tol = ...)
 446 '''
 447     assert tol >= 0
 448     id = E([0,1,0]) 
 449     Q = calculate_Q(E,L)
 450     int_points = integral_points_with_Q(E, L, Q, tol = tol)
 451     if both_signs:
 452         int_points += [-P for P in int_points]
 453     int_points.sort()
 454     return int_points