test_div_1.sage (7050B)
1 from sage.schemes.elliptic_curves.ell_generic import is_EllipticCurve 2 3 # Pre-tests on E and P 4 def suitable( E, P, ell ): 5 if not is_EllipticCurve(E): 6 print "E is not an elliptic curve" 7 return False 8 if E.base_field() != QQ: 9 print "E is not defined over Q" 10 return False 11 if not P in E: 12 print "P is not in E" 13 return False 14 if P.has_finite_order(): 15 print "P has finite order" 16 return False 17 if not is_prime(ell): 18 print "ell is not prime" 19 return False 20 return True 21 22 23 # Stupid auxiliary function. Returns higest power of n dividing m. 24 def val( n, m ): 25 ret = 0 26 while m % (n^(ret+1)) == 0: 27 ret += 1 28 return ret 29 30 # Valuation of l-divisibility of the point P 31 # The parameters h and k describe the ell-part of E(F) (k>=h) 32 def divisibility( P, ell, k, h ): 33 v = 0 34 while v != k-1 and P.is_divisible_by(ell^(v+1)): 35 v += 1 36 return v 37 38 # label is the Cremona label of an elliptic curve over Q of rank >= 1. 39 # ell is a rational prime. 40 def test_label( label, ell ): 41 E = EllipticCurve(label) 42 if E.rank() < 1: 43 print "E has rank 0" 44 return 45 P = E.gens()[0] 46 test( E, P, ell ) 47 48 # E is an elliptic curve over Q and P a point of infinite order on E. 49 # ell is a rational prime. 50 def test( E, P, ell ): 51 if not suitable( E, P, ell ): 52 return 53 54 # Setting up for small primes cases 55 print "Computations for small primes starting..." 56 small_primes = [] 57 vl = [] 58 lpart = [] 59 for p in Primes(): 60 if p > 100: 61 break 62 if p == ell: 63 print "Skipping", p, "because = ell" 64 continue 65 if E.discriminant() % p == 0: 66 print "Skipping", p, "because E has bad reduction" 67 continue 68 if P.reduction(p).order() % ell != 0: 69 print "Skipping", p, "because red of P is infinitely ell-divisible" 70 continue 71 small_primes.append(p) 72 print "Working with prime p =", p 73 74 # "torsion level" 75 F = FiniteField(p) 76 vlj = [] 77 lpj = [] 78 for i in range(3): ### I WOULD LIKE TO CHANGE THIS TO SOMETHING BIGGER 79 print "Torsion level", i, "..." 80 81 # We can build the division fields for increasing powers of l 82 # incrementally. To get a division field, we first compute the 83 # splitting field of the division polynomial. We may be off by 84 # a degree 2 extension. If so, by finite field magic we know 85 # exactly which degree 2 extension we need: the unique one! 86 R.<x> = PolynomialRing(F) 87 F.<a> = E.reduction(p).division_polynomial(ell^i).splitting_field() 88 E_red = E.reduction(p).base_extend(F) 89 # extend if necessary 90 k = val(ell,E_red.gens()[0].order()) 91 h = val(ell,E_red.order()) - k 92 if k < i or h < i: 93 F.<a> = F.extension(2) 94 E_red = E_red.base_extend(F) 95 96 P_red = E_red.point(P.reduction(p)) 97 98 print "[Torsion field computed]" 99 100 # l-part of the abelian group E(F_i) 101 #gg = E_red.abelian_group().gens() 102 #if len( gg ) == 1: 103 # lpj.append( ( val(ell,gg[0].order()), 0 ) ) 104 #else: 105 # lpj.append( (val(ell,gg[0].order()), val(ell,gg[1].order())) ) 106 k = val(ell,E_red.gens()[0].order()) 107 h = val(ell,E_red.order()) - k 108 109 lpj.append( ( k, h ) ) 110 vlj.append(divisibility(P_red,ell,k,h)) 111 112 vl.append(vlj) 113 lpart.append(lpj) 114 115 # Output 116 print "Done!" 117 for i in range(3): 118 print "" 119 print "******************************************************" 120 print "Torsion Level:", i 121 print "" 122 rows = [small_primes, [vl[j][i] for j in range(len(small_primes))], 123 [lpart[j][i] for j in range(len(small_primes))] ] 124 print table(rows) 125 print "" 126 127 # Wrapper for densities(E,P,ell) 128 def densities_label( label, ell ): 129 E = EllipticCurve(label) 130 if E.rank() < 1: 131 print "E has rank 0" 132 return 133 P = E.gens()[0] 134 densities( E, P, ell ) 135 136 def ratio( h, k, d, ell ): 137 x1 = max(k-d,0) 138 y1 = max(h-d,0) 139 x2 = max(k-d-1,0) 140 y2 = max(h-d-1,0) 141 up = ell^(x1+y1)-ell^(x2+y2) 142 down = ell^(h+k) 143 if d == 0: 144 return RDF((up+1)/down) 145 return RDF(up/down) 146 147 148 # Computes the densities dens[n] of primes such that P is ell^n-divisible mod p 149 def densities( E, P, ell ): 150 if not suitable( E, P, ell ): 151 return 152 153 n_primes = 0 154 divisible = [] 155 inf_divisible = 0 156 157 # Array for counting divisibility with specified l-part 158 div_part = [[[0 for i in range(10)] for j in range(10)] for k in range(10)] 159 # total number of elements in the reductions that (do not) form the l-parts 160 tot_l_part = 0 161 tot_non_l_part = 0 162 163 for p in Primes(): 164 if p > 10^2: 165 break 166 if p == ell or E.discriminant() % p == 0: 167 continue 168 169 n_primes +=1 170 #print "Working with prime", p 171 172 E_red = E.reduction(p) 173 N = E_red.order() 174 k = val(ell,E_red.gens()[0].order()) 175 h = val(ell,N) - k 176 temp_non_l_part = (N / (ell^(h+k)))-1 177 tot_non_l_part += temp_non_l_part 178 temp_l_part = N - temp_non_l_part 179 tot_l_part += temp_l_part 180 # Debug 181 # print "Group structure at", p, ":" 182 # print E_red.abelian_group() 183 # print "Our result:", N, (k,h), temp_l_part, temp_non_l_part 184 185 if P.reduction(p).order() % ell != 0: 186 inf_divisible += 1 187 continue 188 189 n = divisibility( P.reduction(p), ell, \ 190 val(ell,E.reduction(p).gens()[0].order()), 0 ) # Wrong parameters but ok 191 while len(divisible) < n+1: 192 divisible.append(0) 193 divisible[n] += 1 194 div_part[k][h][n] += 1 195 196 N = len(divisible) 197 divtotal = [0]*N 198 divtotal[N-1] = divisible[N-1] + inf_divisible 199 for i in range(2,N): 200 divtotal[N-i] = divisible[N-i] + divtotal[N-i+1] 201 202 print "Tested primes:", n_primes 203 for i in range(10): 204 for j in range(10): 205 s = 0 206 for l in range(10): 207 s += div_part[i][j][l] 208 if s != 0: 209 print "l-part (%d,%d):"%(i,j) 210 for l in range(10): 211 if div_part[i][j][l] != 0: 212 print "Exactly %d^%d-divisible: %d, expected %0.3f"%\ 213 (ell,l,div_part[i][j][l],ratio(i,j,l,ell)) 214 215 for n in range(1,N): 216 print "At least %d^%d-divisble: %0.3f density (%d times)"%(ell,n,\ 217 RDF(divtotal[n]/n_primes),divtotal[n]) 218 print "Infinitely-divisble: %0.3f density (%d times), expected %0.3f"%\ 219 (RDF(inf_divisible/n_primes),inf_divisible,\ 220 RDF(tot_non_l_part/(tot_l_part+tot_non_l_part))) 221 222 223