""" bsd.py -- coordinate computation in MAGMA of bounds on Sha. (c) William Stein, July 2004. Here's what this program does: For each level N < 1000, and for each optimal quotient of conductor N with analyatic rank either 0 or 1, run the MAGMA code ShaPrime(E) for 5 minutes. If it completes in 5 minutes, record the answer, and if it doesn't record that it doesn't. The output is written to the file bsd_out, and the format is N isogeny_class a_invariants B time where B is an integer such that any prime p that does not divide B does not divide the order of Sha. time is the number of seconds that the computation takes. The cremona files curves.1-8000, curves.8001-12000, and curves.12001-20000, curves.20001-25000 should be present in the subdirectory cremona, since it contains the definition of the curves. """ DEFAULT_TIMEOUT = 120 import cPickle, operator, os, time def cputime(t=0): if t==0: return time.clock() return time.clock() - t def gcd(a, b=0): if a == 0: return abs(b) if b == 0: return abs(a) if a < 0: a = -a if b < 0: b = -b while b != 0: c = a % b a = b; b = c return a def pad_number(i,max=100): x = str(i) n = len(str(max)) while len(x) < n: x = "0" + x return x def next_curve(F): """ Returns 4-tuple (N, iso, [a1,a2,a3,a4,a6], rank) where N is the conductor, iso the isogeny letter, ai the a-invariants, and rank the rank of the next curve in the file F, where F is one of Cremona's lists of curves. """ assert isinstance(F, file) try: N,iso,num,a_invs,r,_,_ = F.readline().split() except: return None return eval(N),iso,eval(a_invs), eval(r), num def get_from_out(out, var, error_value=0): i = out.find(var + "{") j = out.find("}"+var) if i==-1 or j==-1: return error_value else: try: print out[i+len(var + "{"):j] return eval(out[i+len(var + "{"):j]) except: return error_value def sha_bound(a_invs, timeout=DEFAULT_TIMEOUT): """ Returns a tuple (B,D,t) where B is an integer such that any prime p that does not divide B does not divide the order of Sha, D is the discriminant used in the Kolyvagin computation, and t is the number of seconds that the computation takes. """ cmd = """ SetColumns(0); Attach("bsd.m"); print "timeout = %s"; E := EllipticCurve(%s); B, D := ShaPrimes(E); printf "B{%%o}B, D{%%o}D",B, D; quit; """%(timeout, a_invs) cmd = cmd.replace("L","") # stupid long integer format in Python. a,b = os.popen2("ulimit -t %s; magma"%timeout) a.write(cmd) a.close() out = b.read() i = out.find("Total time: ") if i == -1: t = 0 else: j = i+len("Total time: ") t = float(out[j:j+4]) B = get_from_out(out, "B") D = get_from_out(out, "D") return B,D,t, out def sha_bound2(a_invs, timeout=DEFAULT_TIMEOUT): """ Returns a tuple (B1,B2,D1,D2,t) where B1, B2 are integers such that any prime p that does not divide gcd(B1,B2) does not divide the order of Sha, D1 and D2 are the discriminants used in the Kolyvagin computation, and t is the number of seconds that the computation takes. """ cmd = """ SetColumns(0); Attach("bsd.m"); print "timeout = %s"; E := EllipticCurve(%s); B, B1, B2, D1, D2, F1, z1, h1, F2, z2, h2 := ShaPrimes2(E); printf "B{%%o}B, B1{%%o}B1, B2{%%o}B2, D1{%%o}D1, D2{%%o}D2, F1{%%o}F1, z1{%%o}z1, h1{%%o}h1, F2{%%o}F2, z2{%%o}z2, h2{%%o}h2",B, B1, B2, D1, D2, F1, z1, h1, F2, z2, h2; quit; """%(timeout, a_invs) cmd = cmd.replace("L","") # stupid long integer format in Python. a,b = os.popen2("ulimit -t %s; magma"%timeout) a.write(cmd) a.close() out = b.read() i = out.find("Total time: ") if i == -1: t = 0 else: j = i+len("Total time: ") t = float(out[j:j+4]) B = get_from_out(out, "B") B1 = get_from_out(out, "B1") B2 = get_from_out(out, "B2") D1 = get_from_out(out, "D1") D2 = get_from_out(out, "D2") F1 = get_from_out(out, "F1") z1 = get_from_out(out, "z1") h1 = get_from_out(out, "h1") F2 = get_from_out(out, "F2") z2 = get_from_out(out, "z2") h2 = get_from_out(out, "h2") return B, B1,B2,D1,D2, F1, z1, h1, F2, z2, h2, t, out def analyze_single_curve(curve, timeout=DEFAULT_TIMEOUT, force=False): """ Writes files to disk named data/[N][iso] containing a single line with data about the given elliptic curve. Concatenate these for complete data. If the file has already been created, nothing is done. """ N, iso, a_invs, rank, num = curve filename = "data/%s%s"%(pad_number(N,99999),iso) if os.path.exists(filename+".data") and not force: return F = open(filename+".data","w") a_invs = str(a_invs).replace(" ","").replace("L","") B,D,t,out = sha_bound(a_invs, timeout=timeout) data = "\t".join([str(x) for x in \ [N,iso,a_invs,rank, B, D, t]]) data = data.replace("L","") print data print out F.write(data) F.close() open(filename+".out","w").write(out) def analyze_single_curve2(curve, timeout=DEFAULT_TIMEOUT, \ force = False): """ Writes files to disk named data/[N][iso] containing a single line with data about the given elliptic curve. Concatenate these for complete data. If the file has already been created, nothing is done. """ N, iso, a_invs, rank, num = curve filename = "data/%s%s-2"%(pad_number(N,99999),iso) if not os.path.exists("data/"): os.makedirs("data") if os.path.exists(filename+".data") and not force: return F = open(filename+".data","w") B,B1,B2,D1,D2,F1,z1,h1,F2,z2,h2,t,out = sha_bound2(a_invs, timeout=timeout) data = "\t".join([str(x).replace(" ","") for x in \ [N,iso,num,a_invs,rank, B, B1, B2, D1, D2, F1,z1,h1,F2,z2,h2,t]]) data = data.replace("L","") # annoying large python ints! print data F.write(data) F.close() open(filename+".out","w").write(out) def cremona_file(): assert os.path.exists("cremona") if os.path.exists("cremona/curves.1-25000"): return open("cremona/curves.1-25000","r") files = ["curves.1-8000","curves.8001-12000",\ "curves.12001-20000","curves.20001-25000"] for x in files: assert os.path.exists("cremona/"+x), "cremona/"+x+" must exists" if not os.path.exists("cremona/curves.1-25000"): os.system("cd cremona; cat %s > curves.1-25000"%(" ".join(files))) assert os.path.exists("cremona/curves.1-25000") return open("cremona/curves.1-25000","r") def analyze_curves(x, timeout=DEFAULT_TIMEOUT): """ Runs MAGMA and writes files to disk named [level][iso_number] for each curve whose level is in x. x -- list of integers or a single integer """ if isinstance(x,int): x = [x] assert isinstance(x,list) if not os.path.exists("data/"): os.makedirs("data") F = cremona_file() m = max(x) while True: curve = next_curve(F) N = curve[0] if N == None or N > m: return if N in x: try: analyze_single_curve(curve, timeout=timeout) except: pass class EllipticCurve: def __init__(self, data): if isinstance(data, str): # line from cremona file N,iso,num,a_invs,r,t,deg = data.split() data = {'conductor':int(N), 'isoclass':iso, 'number':int(num), \ 'a_invs':eval(a_invs), 'rank':eval(r), 'torsion':eval(t), \ 'modular_degree':eval(deg)} self.__data = data self.check_for_output() self.check_for_output2() self.check_for_2selmer() def __repr__(self): return "%s%s%s"%(self.conductor(), self.isoclass(), self.number()) def __cmp__(self, other): if not isinstance(other, EllipticCurve): return -1 if self.conductor() < other.conductor(): return -1 if self.conductor() > other.conductor(): return 1 if len(self.isoclass()) < len(other.isoclass()): return -1 if len(self.isoclass()) > len(other.isoclass()): return 1 if self.isoclass() < other.isoclass(): return -1 if self.isoclass() > other.isoclass(): return 1 if self.number() < other.number(): return -1 if self.number() > other.number(): return 1 return 0 def check_for_output(self): """ Look to see if data about this curve has been computed. If so, read in the output of the computation. """ N = self.conductor() assert N != None iso = self.isoclass() assert iso != None file = "%s%s"%(N,iso) print "Attempting to load data from computation for " + file for x in os.listdir("."): if x[:4] == "data": for i in range(5): if os.path.exists( x+"/"+ "0"*i + file + ".data"): self['data'] = open(x+"/"+"0"*i + file+".data","r").read() try: self['log'] = open(x+"/"+"0"*i + file+".out","r").read() except: print "Warning -- unable to open log file " + \ "0"*i + file + ".out" self['log'] = None self.__parse_output() return #endif #endfor #endif #endfor print "Data for run of %s not found."%file def __parse_output(self): """ Parse the data and log files output by MAGMA when computing data about this elliptic curve. """ if self['data'] == None: self.check_for_output() if self['data'] == None or self['log'] == None: return try: N,iso,num,a_invs,rank, B, D, t = self['data'].split("\t") except ValueError: return assert eval(N) == self.conductor() assert iso == self.isoclass() assert eval(num) == self.number() or self.conductor() == 990 assert eval(a_invs) == self.a_invs() assert eval(rank) == self.rank() self['sha_bound'] = eval(B) self['disc'] = eval(D) self['time'] = eval(t) self['computed'] = True def check_for_output2(self): """ Look to see if data about this curve has been computed. If so, read in the output of the computation. """ N = self.conductor() assert N != None iso = self.isoclass() assert iso != None file = "%s%s-2"%(N,iso) print "Attempting to load data from computation2 for " + file for x in os.listdir("."): if x[:4] == "data": for i in range(5): if os.path.exists( x+"/"+ "0"*i + file + ".data"): self['data2'] = open(x+"/"+"0"*i + file+".data","r").read() try: self['log2'] = open(x+"/"+"0"*i + file+".out","r").read() except: print "WARNING -- unable to open log file " + \ "0"*i + file + ".out" self['log2'] = None self.__parse_output2() return #endif #endfor #endif #endfor print "Data for run of %s not found."%file def __parse_output2(self): """ Parse the data and log files output by MAGMA when computing data about this elliptic curve. """ if self['data2'] == None: self.check_for_output2() if self['data2'] == None or self['log2'] == None: return try: N,iso,num,a_invs,rank, B, B1, B2, D1, D2, F1, z1, h1, F2, z2, h2, t = \ self['data2'].split("\t") except ValueError: return assert eval(N) == self.conductor() assert iso == self.isoclass() assert eval(num) == self.number() or self.conductor() == 990 assert eval(a_invs) == self.a_invs() assert eval(rank) == self.rank() self['sha_bound2'] = eval(B) self['sha_bound2_each'] = [eval(B1), eval(B2)] self['disc2'] = [eval(D1), eval(D2)] self['time2'] = eval(t) self['computed2'] = True self['F1'] = F1 self['z1'] = z1 self['h1'] = h1 self['F2'] = F2 self['z2'] = z2 self['h2'] = h2 self['indexes'] = None def indexes(self): """ The odd part of the indexes of the Heegner points in the group E(K)/tors. """ if self['indexes'] != None: return self['indexes'] log = self['log2'] self['indexes'] = [] if log != None: # extract index i = log.find("Index= ") j = log[i:].find("\n") if i != -1 and j != -1: self['indexes'].append(oddpart(int(log[i+6:i+j]))) log = log[i+j:] i = log.find("Index= ") j = log[i:].find("\n") if i != -1 and j != -1: self['indexes'].append(oddpart(int(log[i+6:i+j]))) return self['indexes'] def index(self): X = self.indexes() if X == None or X == []: return None return GCD(X) def torsions_over_K(self): if self['torsk'] != None: return self['torsk'] log = self['log2'] if log != None: # extract torsions i = log.find("torK= ") if i != -1: j = log[i:].find("\n") self['torsk'] = eval(log[i+len("torK= "):i+j]) return self['torsk'] def torsion_over_K(self): T = self.torsions_over_K() if T == None: return None return LCM(T) def check_for_2selmer(self): """ Look to see if the 2-selmer data about this curve has been computed. If so, read in the output of the computation. """ N = self.conductor() assert N != None iso = self.isoclass() assert iso != None file = "%s%s-2selmer"%(N,iso) print "Attempting to load data from 2selmer for " + file for x in os.listdir("."): if x[:4] == "data": for i in range(5): if os.path.exists( x+"/"+ "0"*i + file + ".data"): print "Found 2selmer data file" self['data_2selmer'] = open(x+"/"+"0"*i + file+".data","r").read() try: self['log_2selmer'] = open(x+"/"+"0"*i + file+".out","r").read() except: print "WARNING -- unable to open log file " + \ "0"*i + file + ".out" self['log2_2selmer'] = None self.__parse_2selmer() return #endif #endfor #endif #endfor print "Data for run of %s not found."%file def __parse_2selmer(self): """ Parse the data and log files output by MAGMA when computing 2-selmer group of this elliptic curve. """ if self['data_2selmer'] == None: self.check_for_2selmer() if self['data_2selmer'] == None or self['log_2selmer'] == None: return #try: N,iso,num,a_invs,two_rank, t = self['data_2selmer'].split("\t") #except ValueError: # return assert eval(N) == self.conductor() assert iso == self.isoclass() assert eval(num) == self.number() or self.conductor() == 990 assert eval(a_invs) == self.a_invs() self['rank_2selmer'] = eval(two_rank) self['time_2selmer'] = eval(t) def rank_2selmer(self): return self['rank_2selmer'] def time_2selmer(self): return self['time_2selmer'] def data_2selmer(self): return self['data_2selmer'] def log_2selmer(self): return self['log_2selmer'] def __getitem__(self, key): if key in self.__data.keys(): return self.__data[key] return None def __setitem__(self, key, value): self.__data[key] = value def compute_2selmer(self, timeout=DEFAULT_TIMEOUT, force=False): """ Computes rank of 2-selmer group of this elliptic curve, unless another process is already computing it. The result is stored as a data file, and is not loaded into this curve if another process is computing this 2-selmer group already. """ filename = "data/%s%s-2selmer"%(pad_number(self.level(),99999),self.isoclass()) if os.path.exists(filename+".data") and not force: return cmd = """ SetColumns(0); print "timeout = %s"; E := EllipticCurve(%s); G,f := TwoSelmerGroup(E); printf "r{%%o}r",#Invariants(G); quit; """%(timeout, self.a_invs()) cmd = cmd.replace("L","") # stupid long integer format in Python. a,b = os.popen2("ulimit -t %s; magma"%timeout) a.write(cmd) a.close() out = b.read() i = out.find("Total time: ") if i == -1: t = 0 else: j = i+len("Total time: ") t = float(out[j:j+4]) r = get_from_out(out, "r", error_value=-1) F = open(filename+".data","w") a_invs = str(self.a_invs()).replace(" ","").replace("L","") data = "\t".join([str(x) for x in [self.level(), self.isoclass(), \ self.number(), a_invs, r, t]]) print data print out F.write(data) F.close() open(filename+".out","w").write(out) self.check_for_2selmer() def run_computation(self, timeout=DEFAULT_TIMEOUT, force=False): """ Use magma to compute BSD bound data about this curve. timeout -- the number of seconds before MAGMA is killed force -- if true, then computation is redone, even if it has been done before. """ print "Computing Kolyvagin data for 1 disc using %s."%self curve = [self[k] for k in ['conductor', 'isoclass', 'a_invs', 'rank', 'number']] analyze_single_curve(curve, timeout=timeout, force=force) self.check_for_output() def run_computation2(self, timeout=DEFAULT_TIMEOUT, force=False, force_if_timeout=False): """ Use magma to compute BSD bound data about this curve. timeout -- the number of seconds before MAGMA is killed force -- if true, then computation is redone, even if it has been done before. force_if_timeout -- if true, then if this computation was previously attempted but timed out, recompute even if force is False """ print "Computing Kolyvagin data using up to 2 discs for %s."%self curve = [self[k] for k in ['conductor', 'isoclass', 'a_invs', 'rank', 'number']] if not force and self.log_file2() != "": # computation has been attempted already if not self.failed_because_of_timeout2(): return if self.failed_because_of_timeout2() and not force_if_timeout: return # at this point we *always* want to do the computation analyze_single_curve2(curve, timeout=timeout, force=True) self['timeout2'] = timeout self.check_for_output2() def sha_bound2(self): if self['sha_bound2'] == None: return 0 return self['sha_bound2'] def sha_bound2_each(self): if self['sha_bound2_each'] == None: return [0,0] return self['sha_bound2_each'] def disc2(self): return self['disc2'] def time2(self): return self['time2'] def log_file2(self): if self['log2'] == None: self['log2'] = "" return self['log2'] def conductor(self): return self['conductor'] def level(self): return self['conductor'] def a_invs(self): return self['a_invs'] def isoclass(self): return self['isoclass'] def number(self): return self['number'] def rank(self): return self['rank'] def torsion(self): return self['torsion'] def sha_bound(self): if self['sha_bound'] == None: return 0 return self['sha_bound'] def log_file(self): if self['log'] == None: self.check_for_output() if self['log'] == None: self['log'] = "" return self['log'] def data_file(self): if self['data'] == None: self.check_for_output() return self['data'] def data2_file(self): if self['data2'] == None: self.check_for_output2() return self['data2'] def disc(self): return self['disc'] def modular_degree(self): return self['modular_degree'] def warning(self): """ True if the log file contains the string WARNING. """ return self.log_file().find("WARNING") != -1 def warning2(self): """ True if the log file 2 (up to 2 discs) contains the string WARNING. """ return self.log_file2().find("WARNING") != -1 def time(self): return self['time'] def time_bound(self): """ The time bound that was used in the computation of the possible_sha_primes for this curve. """ if self['time_bound'] == None: return 120 return self['time_bound'] def computed(self): """ Whether or not MAGMA was run yet on this particular curve. """ if self['computed'] == None: return False return self['computed'] def failed_because_of_L(self): if self['log'] == None: self.check_for_output() if self['log'] == None: return False return self['log'].find("bad syntax") != -1 def failed_because_of_timeout(self): if self['log'] == None: self.check_for_output() if self['log'] == None: return False return self['log'].find("Total time") == -1 def failed_because_of_timeout2(self): if self['log2'] == None: self.check_for_output() if self['log2'] == None: return False return self['log2'].find("Total time") == -1 def table_line(self): N = self.conductor() iso = self.isoclass() num = self.number() r = self.rank() B = self.sha_bound2_each() D = self.disc() t = self.time() a_invs = str(self.a_invs()).replace(" ","").replace("L","") notes = "" if self.warning(): notes += " warning" if self.failed_because_of_L(): notes += " failL" if self.failed_because_of_timeout(): notes += " timeout" x = [N, iso, num, r, B, D, t, a_invs, notes] x = [str(y) for y in x] return "\t".join(x) def table2_line(self): N = self.conductor() iso = self.isoclass() num = self.number() r = self.rank() sha = self.bound() sel = self.rank_2selmer() Bs = self.sha_bound2_each() D = self.disc2() t = self.time2() a_invs = str(self.a_invs()).replace("L","") notes = [] if self.warning2(): notes.append("warning") if self.failed_because_of_timeout2(): notes.append("timeout") notes = ",".join(notes) x = [pad(N,5), pad(iso,4), pad(num,4), pad(sha,4), pad(r,4), \ pad(sel,4), pad(t,5), nospace(Bs), nospace(D),\ nospace(a_invs), notes] return " ".join(x) def table2_curve_gen_line(self): if self['F1'] != None: x = [self['F1'],self['z1'], self['h1'], \ self['F2'],self['z2'], self['h2']] else: x = ["NONE" for i in range(6)] x = [nospace(y) for y in x] N = self.conductor() iso = self.isoclass() num = self.number() r = self.rank() x = [pad(N,5),pad(iso,4),pad(num,4), pad(r,4)] + x return (" ".join(x)).replace("L","") def bound(self): """ Combines everything known about the curve and returns a square- free product of primes divisible only by primes that divide the order of Sha, or 0 if nothing is known about Sha. """ b = 0 x = self.sha_bound() if x != None: b = x x = self.sha_bound2() if x != None: b = gcd(x,b) if self.rank() != None and self.rank() > 1: return 0 if b % 2 == 0: r_selmer = self.rank_2selmer() r_curve = self.rank() tors = self.torsion() if r_selmer == None or r_selmer < 0 \ or r_curve == None or tors == None: return b if tors%2 != 0: t = 0 else: t = 1 if r_selmer == (r_curve + t) or r_selmer == (r_curve + t + 1): # Comparing to either "r_curve + t" or "r_curve+t+1" is # OK, because Sha is finite for this curve, hence Sha[2] # has even dimension. The reason we do this is because # the t computed above, which is got only from the order # of the torsion subgroup, could be one too small (e.g., # Z/2 x Z/4 versus Z/8). b = b/2 return b class Database: def __init__(self, name="default", new=False): if isinstance(name, list): x = name do_load, a, b = True, min(x), max(x) name = "%s-%s"%(pad_number(a,25000),pad_number(b,25000)) else: do_load = False if new == True: if os.path.exists(self.__filename(name)): i = 0 while os.path.exists(self.__filename(name + str(i))): i += 1 name = name + "-v" + str(i) self.__name = name if os.path.exists(self.__filename()): print "Loading %s from disk."%name self.__data = cPickle.load(open(self.__filename(),"r")) else: self.__data = {} if do_load: self.load(a,b) self.save() def clear(self): self.__data = {} self.save() def run_computation(self, timeout=DEFAULT_TIMEOUT, force=False, min_level = None, max_level=None): """ Use magma to compute BSD bound data about all curves in database. timeout -- the number of seconds before MAGMA is killed force -- if true, then computation is redone, even if it has been done before. """ L = self.levels() if min_level == None: min_level = min(L) if max_level == None: max_level = max(L) for x in [N for N in L if N >= min_level and N <= max_level]: for E in self[x]: E.run_computation(timeout=timeout, force=force) def run_computation2(self, timeout=DEFAULT_TIMEOUT, force=False, force_if_timeout=False, min_level=None, max_level=None): """ Use magma to compute BSD bound data about this curve. timeout -- the number of seconds before MAGMA is killed force -- if true, then computation is redone, even if it has been done before. """ L = self.levels() if min_level == None: min_level = min(L) if max_level == None: max_level = max(L) for x in [N for N in L if N >= min_level and N <= max_level]: for E in self[x]: print E E.run_computation2(timeout=timeout, force=force, \ force_if_timeout=force_if_timeout) def run_computation_2selmer(self, timeout=DEFAULT_TIMEOUT, force=False, \ min_level=None, max_level=None): """ Use magma to compute BSD bound data about this curve. timeout -- the number of seconds before MAGMA is killed force -- if true, then computation is redone, even if it has been done before. """ L = self.levels() if min_level == None: min_level = min(L) if max_level == None: max_level = max(L) for x in [N for N in L if N >= min_level and N <= max_level]: for E in self[x]: print E E.compute_2selmer(timeout=timeout, force=force) def __filename(self, name=None): if name == None: name = self.name() if not os.path.exists("bsd_data"): os.makedirs("bsd_data") return "bsd_data/bsd_data-%s.python"%name def __len__(self): return len(self.__data) def name(self): return self.__name def rename(self, newname): self.__name = newname self.save() def levels(self): x = self.__data.keys() x.sort() return x def curves(self): """ Returns a flat ordered list of all curves in the database. """ return reduce(operator.add, [self[x] for x in self.levels()]) def load(self, min, max): F = cremona_file() while True: try: L = F.readline() N = eval(L.split()[0]) if N > max: break if N < min: continue except: break if len(L) == 0: break self.include(EllipticCurve(L)) #endwhile self.save() def save(self): cPickle.dump(self.__data,open(self.__filename(),"w")) def __repr__(self): return "Database %s of %s levels (%s)."%(self.name(), len(self.__data), self.__filename()) def __getitem__(self, N): if self.__data.has_key(N): return self.__data[N] return [] def __setitem__(self, N, value): if not self.__data.has_key(N): self.__data[N] = {} self.__data[N] = value def include(self, E): assert isinstance(E, EllipticCurve) x = self[E.level()] for i in range(len(x)): if x[i] == E: x[i] = E return x.append(E) x.sort() self[E.level()] = x def table_complete(self): """ Create a file tables/bsd_bounds.min-max that gives the data in self. """ N = self.levels() if not os.path.exists("tables"): os.makedirs("tables") file = "tables/bsd_bounds.%s-%s"%(min(N), max(N)) F = open(file,"w") data = self.__data K = data.keys() K.sort() for N in K: if N%50==0: print N for E in data[N]: F.write("%s\n"%E.table_line()) def table_shabounds(self): """ Create a file tables/name-shabound that gives the data in self. """ N = self.levels() nm = self.name() if nm == "default": nm = "%s-%s"%(min(N),max(N)) if not os.path.exists("tables"): os.makedirs("tables") file = "tables/%s-shabound"%nm F = open(file,"w") data = self.__data K = data.keys() K.sort() for N in K: if N%50==0: print N for E in data[N]: F.write("%s\n"%E.table2_line()) def table_curvegens(self): """ Create a file tables/curvegens.min-max that gives the curves and generators computed by MAGMA that were used in the computation of bounds on the order of Sha. These could be incorrect, and the purpose of this file is to provide input to a program that computes provably correct generators, for comparison. """ N = self.levels() nm = self.name() if nm == "default": nm = "%s-%s"%(min(N),max(N)) if not os.path.exists("tables"): os.makedirs("tables") file = "tables/%s-curvegens"%nm F = open(file,"w") data = self.__data K = data.keys() K.sort() for N in K: if N%50==0: print N for E in data[N]: F.write("%s\n"%E.table2_curve_gen_line()) def curves_with_warning(self): """ Sub-database of curves such that the log for the sha bound computation using a single discriminant divisible by at least two primes contains a warning. """ d = Database("%s-warning"%self.name(), new=True) for N in self.levels(): if N%50==0: print N for E in self[N]: if E.warning(): d.include(E) d.save() return d def curves_with_warning2(self): """ Sub-database of curves such that the log for the sha bound computation using up to two discriminants contains a warning. """ d = Database("%s-warning2-"%self.name(), new=True) for N in self.levels(): if N%50==0: print N for E in self[N]: if E.warning2(): d.include(E) d.save() return d def curves_with_timeout(self): """ Sub-database of curves whose computation failed because they ran out of time (MAGMA was killed after a few minutes). """ d = Database("%s-timeout"%self.name(), new=True) for N in self.levels(): if N%50==0: print N for E in self[N]: if E.failed_because_of_timeout(): print E d.include(E) d.save() return d def curves_with_L(self): """ Sub-database of curves whose computation failed because the coefficients of Weierstrass equations were python long integers, so python put a stupid L after them. """ d = Database("%s-L"%self.name(), new=True) for N in self.levels(): if N%50==0: print N for E in self[N]: if E.failed_because_of_L(): print E d.include(E) d.save() return d def curves_of_rank(self, rank=0): """ Sub-database of curves of given rank. """ d = Database("%s-rank%s"%(self.name(),rank), new=True) for N in self.levels(): if N%50==0: print N for E in self[N]: if E.rank() == rank: print E d.include(E) d.save() return d def curves_with_index_at_least(self, x=3): """ The sub-database of curves such that the gcd of the indexes of the Heegner point is at least x. If the index or indexes are unknown, the curve is not included in this list. """ d = Database("%s-index%s"%(self.name(),x), new=True) for N in self.levels(): if N%50==0: print N for E in self[N]: if E.index() != None and E.index() >= x: print E d.include(E) d.save() return d def curves_with_index_divbyprime_at_least(self, x=3): """ The sub-database of curves such that the gcd of the indexes of the Heegner point is at least x. If the index or indexes are unknown, the curve is not included in this list. Only primes up to 300 are actually used. """ P = [z for z in [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293] if z >= x] d = Database("%s-index_divbyprime_%s"%(self.name(),x), new=True) for N in self.levels(): if N%50==0: print N for E in self[N]: if E.index() != None: ind = E.index() for p in P: if ind%p == 0: print E d.include(E) break d.save() return d def curves_with_no_bound(self): """ Sub-database of curves whose Sha bound using any computation so far is 0. """ d = Database("%s-no_bound"%self.name(), new=True) for N in self.levels(): if N%50==0: print N for E in self[N]: if E.sha_bound() == 0 and E.sha_bound2() == 0: print E d.include(E) d.save() return d def nospace(x): return str(x).replace(" ","") def pad(x,n): s=str(x) s += " "*(n-len(s)) return s def __GCD_list(v): if len(v) == 0: return 1 if len(v) == 1: return v[0] g = v[0] for i in range(1,len(v)): g = GCD(g, v[i]) return g def GCD(a, b=0): """ The greatest commond divisor of a and b. >>> GCD(97,100) 1 >>> GCD(97 * 10**15, 19**20 * 97**2) 97 """ if isinstance(a,list): return __GCD_list(a) if a == 0: return abs(b) if b == 0: return abs(a) if a < 0: a = -a if b < 0: b = -b while b != 0: c = a % b a = b; b = c return a def oddpart(x): while x%2 == 0 and x: x = x // 2 return x def LCM(a, b=None): """The least common multiple of a and b. >>> LCM(97,100) 9700 >>> LCM(0,2) 0 >>> LCM(-3,-5) 15 """ if isinstance(a,list): return __LCM_list(a) return int(long((a*b)/GCD(a,b))) def __LCM_list(v): if len(v) == 0: return 1 x = v[0] for i in range(1,len(v)): x = LCM(x, v[i]) return x