""" 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); Alarm(%s); Attach("bsd.m"); 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("magma") 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); Alarm(%s); Attach("bsd.m"); 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("magma") 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]]) 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]]) 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['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 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); Alarm(%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("magma") 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): """ 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 using up to 2 discs for %s."%self curve = [self[k] for k in ['conductor', 'isoclass', 'a_invs', 'rank', 'number']] analyze_single_curve2(curve, timeout=timeout, force=force) self.check_for_output2() def sha_bound2(self): if self['sha_bound2'] == None: return 0 return self['sha_bound2'] 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 table_line(self): N = self.conductor() iso = self.isoclass() num = self.number() r = self.rank() B = self.sha_bound() 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 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 new == True: if os.path.exists(self.__filename(name)): i = 0 while os.path.exists(self.__filename(name + str(i))): i += 1 name = name + 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 = {} def clear(self): self.__data = {} self.save() def run_computation(self, timeout=DEFAULT_TIMEOUT, force=False): """ 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. """ for x in self.levels(): for E in self[x]: E.run_computation(timeout=timeout, force=force) def run_computation2(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. """ for x in self.levels(): for E in self[x]: print E E.run_computation2(timeout=timeout, force=force) def run_computation_2selmer(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. """ for x in self.levels(): 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): 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 for N in data.keys(): if N%50==0: print N for E in data[N]: F.write("%s\n"%E.table_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_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