attachment:intpts.sage of days22/cremona/project2

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