CoCalc Shared Filescode / compute_lambda.sage
1load("modular_symbol_map.pyx")
2
3def 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
48def running_average(t):
49    show(stats.TimeSeries(t[:i].mean() for i in range(5,len(t))).plot())
50
51def 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
57def 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
63def 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
72def 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
79def 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 = []
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']:
102
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
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
123def 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
132def 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">
159Rank=%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)