CoCalc Shared Filescode / compute_lambda.sageOpen in CoCalc with one click!
Authors: Jennifer Balakrishnan, Maarten Derickx, JB Nam, James Rickards, William A. Stein, Andrew Sutherland, Nicholas Triantafillou
1
load("modular_symbol_map.pyx")
2
3
def compute_dist(E, d, m_max):
4
M = E.modular_symbol_space(sign=1)
5
assert d%2 == 1
6
N = M.level()
7
ms_map = ModularSymbolMap(M)
8
ms_denom = ZZ(ms_map.denom)
9
inf_zero = M.rational_period_mapping()([oo,0])[0]
10
print ms_denom, inf_zero
11
def f(a,b):
12
return ms_map._eval1(a,b)[0] / ms_denom
13
14
# much more ms, since this code is massively faster...
15
ms = [m for m in prime_range(3,m_max+1) if \
16
gcd(m, N) == 1 and euler_phi(m) % d == 0]
17
18
print 'There are %s primes to use up to %s'%(len(ms), m_max)
19
20
def alphas(m, d):
21
assert d%2 == 1
22
R = Integers(m)
23
Npow = R(N)^((d-1)//2)
24
gen = R(primitive_root(m))
25
n = euler_phi(m)//d
26
b = gen
27
h = gen^d
28
denom = float(sqrt(euler_phi(m)*log(m)))
29
alphas = []
30
for i in range(1, (d-1)//2 + 1):
31
s = 0
32
for j in range(n):
33
period = f((Npow * b^i * h^j).lift(), m) + inf_zero
34
s += period
35
alphas.append(s / denom)
36
return alphas
37
38
t0 = walltime()
39
print "Starting..."
40
data = []
41
progress = len(ms) // 10 + 1
42
for i, m in enumerate(ms):
43
if i % progress == 0:
44
print 'X',; sys.stdout.flush()
45
data += alphas(m, d)
46
return stats.TimeSeries(data)
47
48
def running_average(t):
49
show(stats.TimeSeries(t[:i].mean() for i in range(5,len(t))).plot())
50
51
def plot_histogram(t, bins='auto'):
52
import matplotlib.pyplot as plt
53
plt.figure(figsize=(14,6))
54
plt.hist(t.numpy(), bins=bins)
55
plt.show()
56
57
def report(t):
58
import scipy.stats
59
print "kurtosis: ", scipy.stats.kurtosis(t.numpy(), fisher=False)
60
print "mean: ", t.mean()
61
print "sd: ", t.standard_deviation()
62
63
def drew_data(name):
64
if not name.endswith('.txt'):
65
name += '.txt'
66
filename = os.path.join('data', name)
67
if not os.path.exists(filename):
68
raise RuntimeError("'%s' does not exist"%filename)
69
s = '[' + open(filename).read().replace('\n',',') + ']'
70
return stats.TimeSeries(eval(s))
71
72
def get_drew_data(label, d):
73
v = os.listdir('data')
74
v = [name for name in v if label + '_%s_3_'%d in name]
75
v.sort()
76
name = v[-1]
77
return drew_data(name)
78
79
def drew_report_all(redo=False):
80
v = [name for name in os.listdir('data') if name.endswith('.txt')]
81
def mycmp(name1, name2):
82
c = cmp(int(name1.split('_')[1][:-1]), int(name2.split('_')[1][:-1])) # cond
83
if c:
84
return c
85
c = cmp(name1.split('_')[1], name2.split('_')[1]) # isog class
86
if c:
87
return c
88
c = cmp(int(name1.split('_')[2]), int(name2.split('_')[2])) # d
89
return c
90
91
v.sort(mycmp)
92
w = []
93
links = ''
94
prev = None
95
for name in v:
96
y = parse_name(name)
97
98
if not prev or prev['label'] != y['label'] or prev['min_m'] != y['min_m'] or prev['max_m'] != y['max_m']:
99
if links:
100
w.append('<li>%s</li>'%links)
101
links = ''
102
103
if links == '':
104
links += '%s, %s&le;m&le;%s&nbsp;&nbsp;&nbsp;<a href="%s">d=%s</a>&nbsp;&nbsp;&nbsp;'%(y['label'], y['min_m'], y['max_m'], name+'.html', y['d'])
105
else:
106
links += '<a href="%s">d=%s</a>&nbsp;&nbsp;&nbsp;'%(name+'.html', y['d'])
107
108
prev = y
109
110
if links:
111
w.append('<li>%s</li>'%links)
112
113
s = '\n'.join(w)
114
open("plots/index.html", 'w').write('<body style="word-wrap: break-word;"><ul>%s</ul></body>'%s)
115
if not redo:
116
v = [name for name in v if not os.path.exists(os.path.join('plots', name+'.html'))]
117
print "Generating %s files"%len(v)
118
if len(v) == 0:
119
return
120
for x in drew_report(v):
121
pass
122
123
def parse_name(name):
124
v = name[:-4].split('_') # alpha_63a_13_3_1000000
125
label = v[1]
126
d = ZZ(v[2])
127
min_m = ZZ(v[3])
128
max_m = ZZ(v[4])
129
return {'label':label, 'd':d, 'min_m':min_m, 'max_m':max_m}
130
131
@parallel
132
def drew_report(name):
133
if not name.endswith('.txt'):
134
raise RuntimeError('name must end in .txt')
135
if not os.path.exists('plots'):
136
os.makedirs('plots')
137
v = name[:-4].split('_') # alpha_63a_13_3_1000000
138
label = v[1]
139
d = ZZ(v[2])
140
min_m = ZZ(v[3])
141
max_m = ZZ(v[4])
142
rank = EllipticCurve(label).rank()
143
t = drew_data(name)
144
np = t.numpy()
145
svg = 'plots/%s.svg'%name
146
import matplotlib.pyplot as plt
147
import scipy.stats
148
import scipy.stats.mstats
149
test = scipy.stats.mstats.normaltest(np) # see https://docs.scipy.org/doc/scipy-0.15.1/reference/generated/scipy.stats.mstats.normaltest.html
150
if not os.path.exists(svg):
151
plt.figure(figsize=(12,6))
152
plt.hist(np, bins='auto')
153
plt.savefig(svg)
154
html = open("plots/%s.html"%name,'w')
155
html.write("""
156
<h1 style="text-align:center; font-family:courier">EllipticCurve('%s') &nbsp;&nbsp;&nbsp;&nbsp; d=%s &nbsp;&nbsp;&nbsp;&nbsp;&nbsp; %s&le;m&le;%s</h1>
157
<hr>
158
<pre style="text-align:center; font-size:12pt">
159
Rank=%s Mean=%6.6f SD=%6.6f Kurtosis=%6.6f Samples=%s Min=%s Max=%s k2=%6.6f p=%6.6f
160
</pre>
161
<hr>
162
<img src='%s.svg' style="width:100%%">
163
"""%(label, d, min_m, max_m, rank, t.mean(), t.standard_deviation(),
164
scipy.stats.kurtosis(np, fisher=False), len(t), t.min(), t.max(), test.statistic, test.pvalue, name)
165
)