CoCalc Shared Filescode / maarten-scratch.sagewsOpen in CoCalc with one click!
Authors: Alex Best, Maarten Derickx
%load modular_symbol_map.pyx
Compiling ./modular_symbol_map.pyx...
sage_server.MAX_OUTPUT = 1000000 for mf in ["431a", "431b", "433a"]: fp=open(mf,"w") M = EllipticCurve(mf+"1").modular_symbol_space(sign=1) MSinfz = M.rational_period_mapping()([oo,0])[0] d = M.dimension() assert d == 1 f = ModularSymbolMap(M) X,MSden = f.C._clear_denom() X=[a[0] for a in X] N = M.level() P1 = f.P1 MStab = [[X[P1.index(a,b)] for b in Integers(N)] for a in Integers(N)] fp.write("int MStab_%s[] = {\n"%mf) for i in range(N): fp.write("%d,"%MStab[i][0]) for j in range(1,N): fp.write("%d,"%MStab[i][j]) fp.write("\n") fp.write("};\n\n") fp.write("int MSden_%s = %d;\n"%(mf,MSden)) fp.write("int MSinfz_%s = %d;\n"%(mf,MSinfz)) fp.write("\n") fp.close()
def contfrac_q(a, b): """ Given coprime integers `a`, `b`, compute the sequence `q_i` of denominators in the continued fraction expansion of `a/b` """ qi = [1] if a == 0: return qi if b == 0: raise ZeroDivisionError # one iteration isn't needed, since we are only computing the q_i. a,b = b, a%b i=1 while b: if i >= 2: qi.append(qi[-1]*(a//b) + qi[-2]) else: qi.append(a//b) a,b = b, a%b i += 1 return qi def eval(a,b,tab): q = contfrac_q(a,b) v = 0 sign = 1 for i in range(1,len(q)): v += tab[(sign*q[i])%N][q[i-1]%N] sign *= -1 return v def alphas(m,d,N,tab): assert is_prime(m) and d%2 == 1 R = Integers(m) Npow = R(N)^((d-1)//2) gen = R(primitive_root(m)) n = (m-1)//d b = gen h = gen^d return [sum([eval((Npow*b^i*h^j).lift(),m,tab) for j in range(n)]) for i in range(1,(d-1)//2+1)] def normalize_alphas(alist,m,MSden,MSinfz): # should be normalized sqrt(euler-phi(m)/d * C_E + D_E), but computing C_E and D_E is a pain, and d is fixed return [(a/MSden+((m-1)//d)*MSinfz)/float(sqrt((m-1)*log(m))) for a in alist] MOMENTS=10 def moments(alist): M=[0.0 for i in range(MOMENTS+1)] for a in alist: apow = 1.0 for i in range(MOMENTS+1): M[i]+=apow apow *= a return [M[i]/len(alist) for i in range(MOMENTS+1)]
d = 3 ms = [m for m in prime_range(3,100) if gcd(m,N)==1 and m % d == 1] a = [] for m in ms: alist = alphas(m,d,N,MStab) print m, alist, normalize_alphas(alist,m,MSden,MSinfz) a += normalize_alphas(alist,m,MSden,MSinfz) print len(a), moments(a)
7 [2] [0.5853194810641169] 13 [-2] [-0.3604956057297047] 19 [2] [0.27472138285709213] 31 [4] [0.39409418584560885] 37 [-8] [-0.7016651856311937] 43 [2] [0.1591263605933982] 61 [4] [0.25469316510042145] 67 [-2] [-0.12005794395973912] 73 [-2] [-0.11379203050314328] 79 [6] [0.32500546873918773] 97 [20] [0.9543607031040101] 11 [1.00000000000000, 0.150119089225459, 0.211787029485880, 0.0737345846554369, 0.113856308766378, 0.0636422712804702, 0.0839059400262076, 0.0601951443657269, 0.0692500441872164, 0.0567079792107614, 0.0600532480700175]
%time load("compute_lambda.sage") a=compute_dist(EllipticCurve("37a1"),29,1000000) plot_histogram(a) print "" print "%d:"%len(a), M = moments(a) for i in range(1,MOMENTS+1): print "%.6f "%M[i],
2 1 There are 2816 primes to use up to 1000000 Starting... X X X X X X X X X X X
39424: 0.000072 0.636646 0.000728 1.645757 -0.010270 9.505303 1.751655 111.743050 109.795762 2240.739724 CPU time: 2827.92 s, Wall time: 2828.03 s
capture?
File: /projects/sage/sage-7.5/local/lib/python2.7/site-packages/smc_sagews/sage_salvus.py Signature : capture(self, code=None, stdout=None, stderr=None, append=False, echo=False) Docstring : Capture or ignore the output from evaluating the given code. (SALVUS only). Use capture as a block decorator by placing either %capture or %capture(optional args) at the beginning of a cell or at the beginning of a line. If you use just plain %capture then stdout and stderr are completely ignored. If you use %capture(args) you can redirect or echo stdout and stderr to variables or files. For example if you start a cell with this line: %capture(stdout='output', stderr=open('error','w'), append=True, echo=True) then stdout is appended (because append=True) to the global variable output, stderr is written to the file 'error', and the output is still displayed in the output portion of the cell (echo=True). INPUT: * stdout -- string (or object with write method) to send stdout output to (string=name of variable) * stderr -- string (or object with write method) to send stderr output to (string=name of variable) * append -- (default: False) if stdout/stderr are a string, append to corresponding variable * echo -- (default: False) if True, also echo stdout/stderr to the output cell.
File: /projects/sage/sage-7.5/local/lib/python2.7/site-packages/smc_sagews/sage_server.py Docstring : Cell execution state object and wrapper for access to special SageMathCloud functionality. An instance of this object is created each time you execute a cell. It has various methods for sending different types of output messages, links to files, etc. Type 'help(smc)' for more details. OUTPUT LIMITATIONS -- There is an absolute limit on the number of messages output for a given cell, and also the size of the output message for each cell. You can access or change those limits dynamically in a worksheet as follows by viewing or changing any of the following variables: sage_server.MAX_STDOUT_SIZE # max length of each stdout output message sage_server.MAX_STDERR_SIZE # max length of each stderr output message sage_server.MAX_MD_SIZE # max length of each md (markdown) output message sage_server.MAX_HTML_SIZE # max length of each html output message sage_server.MAX_TEX_SIZE # max length of tex output message sage_server.MAX_OUTPUT_MESSAGES # max number of messages output for a cell. And: sage_server.MAX_OUTPUT # max total character output for a single cell; computation # terminated/truncated if sum of above exceeds this.
"%s/431a.txt"%os.environ['HOME']
'/projects/68c8b2b8-03ba-44d4-a0d1-5d771c8cb465/431a.txt'
fp=open("foo.txt","w") fp.write("hi") fp.close()
d47c1e60-002b-477f-b5c9-42a39c7b7711