attachment:sha_fast.sage of days22/wuthrich

   1 def minPolyGen(conductor,degree):
   2 	"""
   3 	Give an integer m for which the multiplicative group of
   4 	of ZZ/mZZ is cyclic then for each divisor d of euler_phi(m), there
   5 	will be a unique subfield of Q(zeta_m) of degree d. This returns
   6 	this polynomial which generates such an extension. 
   7 	
   8 	EXAMPLE:
   9 	  m=7, d=3
  10 	  K.<a> = NumberField(minPolyGen(7,3))
  11 	"""
  12 	n = conductor
  13 	d = degree
  14 
  15 	# check that the Z/nZ has cyclic multiplicative group
  16 	if not n % 2 == 0 and not n.is_prime_power():
  17 		raise ValueError, 'Invalid input because (ZZ/%sZZ)* is not cyclic' % n
  18 	if n % 2 == 0:
  19 		nprime = Integer(n/2)
  20 		if not nprime.is_prime_power() or nprime % 4 == 0:
  21 			raise ValueError, 'Invalid input because (ZZ/%sZZ)* is not cyclic' % n
  22  
  23 	# check that there will be such a field of degree d in side QQ(zeta_n)
  24 	if euler_phi(n) % d != 0:
  25 		raise ValueError, 'No field exists because %s does not divide %s=phi(%s)' % (d,euler_phi(n),n)
  26 
  27 	f = euler_phi(n)/d
  28         R = IntegerModRing(n)
  29         g = R.unit_gens()[0]
  30         zetap = CC.zeta(n)
  31 	
  32 	# create a list alpha of all the Galois conjugates
  33         alpha = []
  34         for i in range(d):
  35                 alpha.append(0)
  36                 for j in range(f):
  37                         alpha[i] += zetap^(Integer(g^(d*j+i)))
  38 
  39         S.<x> = ZZ[]
  40         the_poly = prod(x - a for a in alpha)
  41 	coeff = [CC(x).real_part().round() for x in the_poly.coefficients()]
  42 	new_poly = S(0)
  43         for i in range(len(coeff)):
  44                 new_poly += coeff[i]*S.gen()^i
  45 	return new_poly
  46 
  47 def shaOrderFast(E,K,mod_symb,m,precision=10^(-10)):
  48 	"""
  49 	E = EllipticCurve/Q we want #Sha(E/K) for
  50 	K = Field to check over
  51 	mod_symb = modular symbols of E
  52 	m = conductor of K (should be 1 mod d)
  53 	precision = a precision to which we will consider something an integer
  54 
  55 	Specifically, the below algorithm will return a complex number, which should be an
  56 	integer, that represents the order of Sha. If the CC-number is more than <precision>
  57 	away then it will raise an exception	
  58 	"""
  59 	print '\t checking conductor %s on curve %s' % (m,E.cremona_label())
  60 	if E.conductor() % m == 0:
  61 		raise ValueError, 'field conductor m=%s was not coprime to E.conductor()=%s' % (m,E.conductor()) 
  62 	d = K.degree()
  63 	EK = E.change_ring(K)
  64 	tor_order = EK.torsion_subgroup().order()^2
  65 	tamagawa_factor = EK.tamagawa_product_bsd()
  66 	product_result = 1
  67 	R = IntegerModRing(m)
  68 	g = R.unit_gens()[0]
  69 	z = CC.zeta(m-1)
  70 	for j in range(d):
  71 		# creat ell_chi
  72 		if j == 0:
  73 			product_result *= mod_symb(0)
  74 		else:
  75 			product_result *= sum(z^((m-1)*j*k/d)*mod_symb(Integer(g^k)/m) for k in range(m-1))
  76 	real_result = product_result*tor_order/(tamagawa_factor*E.real_components()^d)
  77 	int_result = CC(product_result*tor_order/tamagawa_factor).real_part().round()/E.real_components()^d
  78 	if (real_result - int_result).abs() > precision:
  79 		raise ValueError, 'the result %s was too far from its nearest integer %s' % (real_result, int_result)
  80 	return int_result
  81 
  82 def sha_fast(p,field_conductor_bound=100,curve_conductor_bound=20,curve_conductor_lower_bound=11,spacing=30,filename='sha_fast_data.txt', precision=10^(-10)):
  83 	"""
  84 	p = degree of fields K to consider sha(E/K)
  85 	field_conductor_bound = bound on the conductor of the number field
  86 	curve_conductor_bound = UPPER bound on conductor of elliptic_curves to consider
  87 	curve_conductor_lower_bound = LOWER bound on conductor of elliptic_curves to consider
  88 	spacing = formatting
  89 	filename = filename if you want to specify one
  90 	"""
  91 	if filename == 'sha_fast_data.txt':
  92 		filename = 'sha_data_%s_%s_%s_%s.txt' % (p, field_conductor_bound, curve_conductor_lower_bound, curve_conductor_bound)
  93 	candidates = [q for q in prime_range(field_conductor_bound) if q % p == 1]
  94         print 'Candidates field conductors initialized...'
  95 	fields=[NumberField(minPolyGen(q,p),'a') for q in candidates]
  96 	print 'Fields initialized...'
  97 	file = open(filename, 'a')
  98 	print 'Writing to file %s' % filename
  99 	file.write('Data for fields of degree %s of conductor < %s with curves having conductor between %s and %s\n' % (p, field_conductor_bound, curve_conductor_lower_bound, curve_conductor_bound))
 100 	file.write('%s %s %s %s\n' % ('Curve label'.ljust(spacing), '#Sha(E/K)'.ljust(spacing), 'Field conductor'.ljust(spacing), 'Field degree'.ljust(spacing)))
 101 	for E in CremonaDatabase().iter_optimal([curve_conductor_lower_bound..curve_conductor_bound]):
 102 		print 'Beginning curve %s' % E.cremona_label()
 103 		try:
 104 			M = E.modular_symbol()
 105 			for q in candidates:
 106                         	to_write = ''                        
 107 				try:
 108                                 	shaOrder = shaOrderFast(E,fields[candidates.index(q)],M,q, precision)
 109                                 	to_write = '%s %s %s %s\n' % (str(E.cremona_label()).ljust(spacing), is_square(shaOrder) and str(shaOrder).ljust(spacing) or (str(shaOrder)+'***').ljust(spacing), str(q).ljust(spacing), str(fields[candidates.index(q)].degree()))
 110                         	except ValueError as detail:
 111                                 	to_write = '^^Curve %s threw exception:%s\n' % (E.cremona_label(), detail)
 112                         	except:
 113                                 	to_write = '^^Curve %s threw exception: unrecorded\n' % E.cremona_label()
 114 				finally:
 115                                 	file.write(to_write)
 116 
 117 		except:
 118 			file.write('Curve %s did not compute its modular symbols correctly\n' % E.cremona_label())	
 119 	file.close()
 120 	print 'Finished'