load("modular_symbol_map.pyx")
def compute_dist(E, d, m_max):
M = E.modular_symbol_space(sign=1)
assert d%2 == 1
N = M.level()
ms_map = ModularSymbolMap(M)
ms_denom = ZZ(ms_map.denom)
inf_zero = M.rational_period_mapping()([oo,0])[0]
print ms_denom, inf_zero
def f(a,b):
return ms_map._eval1(a,b)[0] / ms_denom
ms = [m for m in prime_range(3,m_max+1) if \
gcd(m, N) == 1 and euler_phi(m) % d == 0]
print 'There are %s primes to use up to %s'%(len(ms), m_max)
def alphas(m, d):
assert d%2 == 1
R = Integers(m)
Npow = R(N)^((d-1)//2)
gen = R(primitive_root(m))
n = euler_phi(m)//d
b = gen
h = gen^d
denom = float(sqrt(euler_phi(m)*log(m)))
alphas = []
for i in range(1, (d-1)//2 + 1):
s = 0
for j in range(n):
period = f((Npow * b^i * h^j).lift(), m) + inf_zero
s += period
alphas.append(s / denom)
return alphas
t0 = walltime()
print "Starting..."
data = []
progress = len(ms) // 10 + 1
for i, m in enumerate(ms):
if i % progress == 0:
print 'X',; sys.stdout.flush()
data += alphas(m, d)
return stats.TimeSeries(data)
def running_average(t):
show(stats.TimeSeries(t[:i].mean() for i in range(5,len(t))).plot())
def plot_histogram(t, bins='auto'):
import matplotlib.pyplot as plt
plt.figure(figsize=(14,6))
plt.hist(t.numpy(), bins=bins)
plt.show()
def report(t):
import scipy.stats
print "kurtosis: ", scipy.stats.kurtosis(t.numpy(), fisher=False)
print "mean: ", t.mean()
print "sd: ", t.standard_deviation()
def drew_data(name):
if not name.endswith('.txt'):
name += '.txt'
filename = os.path.join('data', name)
if not os.path.exists(filename):
raise RuntimeError("'%s' does not exist"%filename)
s = '[' + open(filename).read().replace('\n',',') + ']'
return stats.TimeSeries(eval(s))
def get_drew_data(label, d):
v = os.listdir('data')
v = [name for name in v if label + '_%s_3_'%d in name]
v.sort()
name = v[-1]
return drew_data(name)
def drew_report_all(redo=False):
v = [name for name in os.listdir('data') if name.endswith('.txt')]
def mycmp(name1, name2):
c = cmp(int(name1.split('_')[1][:-1]), int(name2.split('_')[1][:-1]))
if c:
return c
c = cmp(name1.split('_')[1], name2.split('_')[1])
if c:
return c
c = cmp(int(name1.split('_')[2]), int(name2.split('_')[2]))
return c
v.sort(mycmp)
w = []
links = ''
prev = None
for name in v:
y = parse_name(name)
if not prev or prev['label'] != y['label'] or prev['min_m'] != y['min_m'] or prev['max_m'] != y['max_m']:
if links:
w.append('<li>%s</li>'%links)
links = ''
if links == '':
links += '%s, %s≤m≤%s <a href="%s">d=%s</a> '%(y['label'], y['min_m'], y['max_m'], name+'.html', y['d'])
else:
links += '<a href="%s">d=%s</a> '%(name+'.html', y['d'])
prev = y
if links:
w.append('<li>%s</li>'%links)
s = '\n'.join(w)
open("plots/index.html", 'w').write('<body style="word-wrap: break-word;"><ul>%s</ul></body>'%s)
if not redo:
v = [name for name in v if not os.path.exists(os.path.join('plots', name+'.html'))]
print "Generating %s files"%len(v)
if len(v) == 0:
return
for x in drew_report(v):
pass
def parse_name(name):
v = name[:-4].split('_')
label = v[1]
d = ZZ(v[2])
min_m = ZZ(v[3])
max_m = ZZ(v[4])
return {'label':label, 'd':d, 'min_m':min_m, 'max_m':max_m}
@parallel
def drew_report(name):
if not name.endswith('.txt'):
raise RuntimeError('name must end in .txt')
if not os.path.exists('plots'):
os.makedirs('plots')
v = name[:-4].split('_')
label = v[1]
d = ZZ(v[2])
min_m = ZZ(v[3])
max_m = ZZ(v[4])
rank = EllipticCurve(label).rank()
t = drew_data(name)
np = t.numpy()
svg = 'plots/%s.svg'%name
import matplotlib.pyplot as plt
import scipy.stats
import scipy.stats.mstats
test = scipy.stats.mstats.normaltest(np)
if not os.path.exists(svg):
plt.figure(figsize=(12,6))
plt.hist(np, bins='auto')
plt.savefig(svg)
html = open("plots/%s.html"%name,'w')
html.write("""
<h1 style="text-align:center; font-family:courier">EllipticCurve('%s') d=%s %s≤m≤%s</h1>
<hr>
<pre style="text-align:center; font-size:12pt">
Rank=%s Mean=%6.6f SD=%6.6f Kurtosis=%6.6f Samples=%s Min=%s Max=%s k2=%6.6f p=%6.6f
</pre>
<hr>
<img src='%s.svg' style="width:100%%">
"""%(label, d, min_m, max_m, rank, t.mean(), t.standard_deviation(),
scipy.stats.kurtosis(np, fisher=False), len(t), t.min(), t.max(), test.statistic, test.pvalue, name)
)
