CoCalc Shared Filescode / drew-scratch.sagews
%load 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 [10] [2.0486181837244093] 13 [-40] [-2.8839648458376375] 19 [20] [2.197771062856737] 31 [0] [0.9852354646140221] 37 [-90] [-2.8943688907286744] 43 [-50] [-0.87519498326369] 61 [0] [1.2734658255021074] 67 [0] [1.3206373835571303] 73 [-30] [0.5120641372641447] 79 [-10] [1.1375191405871572] 97 [-170] [-2.5290558632256266] 11 [1.00000000000000, 0.0266115104590982, 3.52512683066233, -3.53963545770662, 20.9232254851203, -37.1832518540298, 147.698802179072, -328.021994193367, 1114.46611503492, -2766.60791251367, 8702.78902798087]
%time
a=compute_dist(EllipticCurve("11a1"),2001,10000)
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 1 primes to use up to 10000 Starting... X
1000: -0.000137 0.010476 -0.000002 0.000332 0.000001 0.000016 0.000000 0.000001 0.000000 0.000000 CPU time: 0.25 s, Wall time: 0.26 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 :
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︠