Sharedwww / talks / mazur-explicit-formula / explicit.pyxOpen in CoCalc
Author: William A. Stein
1
#########################################################
2
#
3
# How Explicit is the Explict Formula
4
#
5
# (c) William Stein, 2013
6
#########################################################
7
8
import math, os, sys
9
10
from sage.all import prime_range, line, tmp_dir, parallel, text, cached_method, nth_prime, load, EllipticCurve, walltime
11
12
from scipy.special import expi as Ei
13
14
cdef extern from "math.h":
15
double log(double)
16
double sqrt(double)
17
double cos(double)
18
double sin(double)
19
20
cdef double pi = math.pi
21
22
################################################
23
# The raw data -- means
24
################################################
25
"""
26
sage: dp = DataPlots("160a", 10^6)
27
sage: v = dp.data(num_points=5000)
28
sage: plot_step_function(v['raw']['mean'],color='red',thickness=.5) + plot_step_function(v['raw']['delta'], thickness=.5)
29
"""
30
31
def data(E, B, aplist, num_points=None, log_X_scale=True, verbose=True):
32
"""
33
Return graphs of the raw, medium, and well-done delta's,
34
along with their means. By a graph, we mean a list of
35
pairs (X, f(X)), where the function is sampled at various
36
values X.
37
38
INPUT:
39
40
- E -- an elliptic curve over QQ
41
- B -- positive integer
42
- aplist -- list of Fourier coefficients a_p
43
- log_X_scale -- if true, scale x-axis on a logarithmic scale,
44
so the point (x,y) in the returned list indicates that
45
mean...(exp(x)) = y. It is important to do this scaling
46
in this function when *computing* the mean, rather than taking
47
the returned data from this function and then rescaling, because
48
we have billions of sample points and this function only returns
49
a few thousand.
50
- verbose -- print estimate of remaining time to do the computation
51
"""
52
if verbose:
53
print "Time remaining:",
54
tm = walltime()
55
if B is None:
56
if aplist is not None:
57
B = nth_prime(len(aplist))+1
58
else:
59
B = 1000
60
if aplist is None:
61
aplist = E.aplist(B)
62
primes = prime_range(B)
63
assert len(aplist) == len(primes)
64
M = len(aplist)
65
if num_points is None or num_points > M:
66
num_points = len(aplist)
67
cdef int record_modulus = M//num_points
68
69
delta_raw = []; delta_medium = []; delta_well = []
70
mean_raw = []; mean_medium = []; mean_well = []
71
72
cdef double ap, p, last_p=0, X, last_X = 0, \
73
sum_raw = 0, sum_medium = 0, sum_well = 0, \
74
last_sum_raw = 0, last_sum_medium = 0, last_sum_well = 0, \
75
integral_raw = 0, integral_medium = 0, integral_well = 0, \
76
gamma_p = 0, last_gamma_p = 0, EilogX = 0, last_EilogX = 0, \
77
length, plot_X_position, next_plot_X_position, plot_X_delta
78
79
next_plot_X_position = 0
80
if log_X_scale:
81
plot_X_delta = log(primes[-1])/num_points
82
else:
83
plot_X_delta = primes[-1]/num_points
84
85
cdef int i = -1, cnt = 0
86
for ap in aplist:
87
i += 1
88
p = primes[i]
89
X = p
90
91
# Update the sums:
92
# raw sum
93
if ap < 0:
94
gamma_p = -1
95
elif ap > 0:
96
gamma_p = 1
97
98
last_sum_raw = sum_raw
99
sum_raw += gamma_p
100
101
# used below
102
logX = log(X)
103
sqrtX = sqrt(X)
104
105
# medium sum
106
last_sum_medium = sum_medium
107
sum_medium += ap/sqrtX
108
109
# well-done sum
110
last_sum_well = sum_well
111
sum_well += ap*logX/p
112
113
# Update the integrals that appear in the mean
114
if i > 0:
115
length = p - last_p
116
# raw mean -- an integral of log(X)/(X*sqrt(X)) is -2*log(X)/sqrt(X) - 4/sqrt(X)
117
integral_raw += ((-2*logX/sqrtX - 4/sqrtX) - (-2*last_logX/last_sqrtX - 4/last_sqrtX)) * last_sum_raw
118
119
# medium mean -- an integral of log(X)/(X*sqrt(X)) is -2*log(X)/sqrt(X) - 4/sqrt(X)
120
integral_medium += ((-2*logX/sqrtX - 4/sqrtX) - (-2*last_logX/last_sqrtX - 4/last_sqrtX)) * last_sum_medium
121
122
# well done mean -- an integral of 1/(X*log(X)) is log(log(X))
123
loglogX = log(logX)
124
integral_well += (loglogX - last_loglogX) * last_sum_well
125
126
last_sqrtX = sqrtX
127
last_logX = logX
128
last_loglogX = loglogX
129
last_X = X
130
last_p = p
131
132
# Finally, record next data point, if it is time to do so...
133
if log_X_scale:
134
plot_X_position = log(X)
135
else:
136
plot_X_position = X
137
138
if verbose and record_modulus >= 1000 and (i%(10*record_modulus)==0):
139
per = float(i)/M
140
if per > 0.1:
141
print "%.1f"%(walltime(tm)*(1-per)/per),
142
sys.stdout.flush()
143
144
if plot_X_position >= next_plot_X_position:
145
next_plot_X_position += plot_X_delta
146
147
delta_raw.append((plot_X_position, sum_raw))
148
delta_medium.append((plot_X_position, sum_medium*logX/sqrtX))
149
delta_well.append((plot_X_position, sum_well/logX))
150
151
mean_raw.append((plot_X_position, integral_raw/logX))
152
mean_medium.append((plot_X_position, integral_medium/logX))
153
mean_well.append((plot_X_position, integral_well/logX))
154
155
if verbose:
156
print "\n-- Done computing raw plot data; returning."
157
158
return {'raw' : {'delta' : delta_raw, 'mean': mean_raw},
159
'medium' : {'delta' : delta_medium, 'mean': mean_medium},
160
'well' : {'delta' : delta_well, 'mean': mean_well} }
161
162
def theoretical_mean(r, vanishing_symmetric_powers):
163
"""
164
INPUT:
165
166
"""
167
mean = 2/pi - 16/(3*pi)*r
168
for n, order in vanishing_symmetric_powers:
169
assert n%2 != 0
170
k = (n-1)/2
171
mean += (4/pi) * (-1)**(k+1)*(1/(2*k+1) + 1/(2*k+3))*order
172
return mean
173
174
def draw_plot(E, B, vanishing_symmetric_powers=None):
175
if vanishing_symmetric_powers is None:
176
vanishing_symmetric_powers = []
177
mean = theoretical_mean(E.rank(), vanishing_symmetric_powers)
178
d, running_average = data(E, B)
179
g = line(d)
180
g += line([(0,mean), (d[-1][0],mean)], color='darkred')
181
g += line(running_average, color='green')
182
return g
183
184
185
############################################################
186
# Plots of data suitable for display (so number of recorded
187
# sample points is smaller), which benefit from having a
188
# pre-computed aplist table.
189
############################################################
190
191
class DataPlots(object):
192
def __init__(self, lbl, B, data_path=None):
193
self.data_path = data_path
194
self.B = B
195
self.E = EllipticCurve(lbl)
196
if self.data_path is None:
197
print "Computing aplist"
198
self.aplist = self.E.aplist(self.B)
199
else:
200
print "Loading aplist"
201
self.aplist = load('%s/%s-aplist-%s.sobj'%(data_path,lbl,B))
202
print "done"
203
204
@cached_method
205
def data(self, num_points=1000,log_X_scale=True, verbose=True): # num_points = number of sample points in output plot
206
print "Computing raw data of plots"
207
d = data(self.E, B=self.B, aplist=self.aplist, num_points=num_points,
208
log_X_scale=log_X_scale, verbose=verbose)
209
print "Done computing raw data"
210
return d
211
212
############################################################
213
# Plots of error term got by taking partial sum over zeros
214
############################################################
215
216
class OscillatoryTerm(object):
217
"""
218
The Oscillatory Term is a real-valued function of a real variable
219
220
S_n(X) = (1/log(X)) * sum(X^(i*gamma_j)/(i*gamma_j) for j<=n)
221
222
where gamma are the imaginary parts of zeros on the critical strip.
223
224
Return object that when evaluated at a specific number X,
225
returns the sequence S_1(X), S_2(X), ..., of partial sums.
226
"""
227
def __init__(self, E, zeros):
228
self.E = E
229
if not isinstance(zeros, list):
230
zeros = E.lseries().zeros(zeros)
231
self.zeros = [float(x) for x in zeros if x>0] # only include *complex* zeros.
232
233
def __call__(self, double X):
234
cdef int i
235
cdef double logX = log(X), running_sum = 0
236
result = []
237
for i in range(len(self.zeros)):
238
running_sum += cos(logX*self.zeros[i])/self.zeros[i]/logX
239
result.append(running_sum)
240
return result
241
242
def plot(self, double X, **kwds):
243
v = self(X)
244
return line(enumerate(v), **kwds)
245
246
def animation_svg(self, Xvals, output_path, ncpus=1, **kwds):
247
if not isinstance(Xvals, list):
248
raise TypeError, "Xvals must be a list"
249
if not os.path.exists(output_path):
250
os.makedirs(output_path)
251
252
print("Rendering %s frames to %s using %s cpus"%(len(Xvals), output_path, ncpus))
253
254
@parallel(ncpus)
255
def f(X):
256
return self(X)
257
258
V = list(f(Xvals))
259
ymax = max([max(v[1]) for v in V])
260
ymin = min([min(v[1]) for v in V])
261
fnames = []
262
for X, v in V:
263
frame = ( line(enumerate(v), ymax=ymax, ymin=ymin, **kwds) +
264
text("X = %s"%X[0], (len(v)//6, ymax/2.0), color='black', fontsize=16) )
265
fname = "%s.svg"%X[0]
266
frame.save(os.path.join(output_path, fname))
267
fnames.append(fname)
268
269
def animation_pdf(self, Xvals, output_pdf, ncpus=1, **kwds):
270
if not isinstance(Xvals, list):
271
raise TypeError, "Xvals must be a list"
272
273
print("using %s cpus"%ncpus)
274
@parallel(ncpus)
275
def f(X):
276
return self(X)
277
278
V = list(f(Xvals))
279
ymax = max([max(v[1]) for v in V])
280
ymin = min([min(v[1]) for v in V])
281
path = tmp_dir()
282
fnames = []
283
for X, v in V:
284
frame = ( line(enumerate(v), ymax=ymax, ymin=ymin, **kwds) +
285
text("X = %s"%X[0], (len(v)//6, ymax/2.0), color='black', fontsize=16) )
286
fname = "%s.pdf"%X[0]
287
frame.save(os.path.join(path, fname))
288
fnames.append(fname)
289
cmd = "cd '%s' && gs -dNOPAUSE -sDEVICE=pdfwrite -sOUTPUTFILE='%s' -dBATCH %s"%(
290
path, os.path.abspath(output_pdf), ' '.join(fnames))
291
print(cmd)
292
os.system(cmd)
293
294
295
296
##############################################################
297
# Plots of mean value of error term involving sums over zeros
298
##############################################################
299
300
def mean_zero_sum_plot(list zeros, int n, double Xmax):
301
"""
302
Plot this function on a log scale up to X at n sample points:
303
304
-2/log(X) * sum cos(gamma*log(X))/gamma^2
305
306
where g runs over zeros, which should be the positive imaginary
307
parts of the first few zeros of an elliptic curve L-function.
308
Returns list of pairs of doubles (x,y), which are points on the plot
309
"""
310
cdef list v = []
311
zeros = [float(x) for x in zeros if x>0]
312
313
# start at X=2
314
cdef double s, gamma, logX=log(2), logXmax = log(Xmax)
315
cdef double delta = logXmax / n
316
317
while logX <= logXmax:
318
s = 0
319
for gamma in zeros:
320
s += cos(gamma*logX)/(gamma*gamma)
321
v.append((logX, -2*s/logX))
322
logX += delta
323
324
return v
325
326
def zero_sum_plot(list zeros, int n, double Xmax):
327
"""
328
Plot this function on a log scale up to X at n sample points:
329
330
f(X) = (2/log(X)) * sum sin(gamma*log(X))/gamma,
331
332
where gamma runs over positive imaginary parts of the first few
333
zeros of an elliptic curve L-function.
334
335
OUTPUT: list of pairs (log(X), f(X)), which are points on the plot.
336
"""
337
cdef list v = []
338
zeros = [float(x) for x in zeros if x>0]
339
340
# start at X=2
341
cdef double s, gamma, logX=log(2), logXmax = log(Xmax)
342
cdef double delta = logXmax / n
343
344
while logX <= logXmax:
345
s = 0
346
for gamma in zeros:
347
s += sin(gamma*logX)/gamma
348
v.append((logX, 2*s/logX))
349
logX += delta
350
351
return v
352
353
def zero_sum_no_log_plot(list zeros, int n, double Xmax):
354
"""
355
Plot this function on a log scale up to X at n sample points:
356
357
f(X) = 2 * sum sin(gamma*log(X))/gamma,
358
359
where gamma runs over positive imaginary parts of the first few
360
zeros of an elliptic curve L-function.
361
362
OUTPUT: list of pairs (log(X),f(X)), which are points on the plot.
363
"""
364
cdef list v = []
365
zeros = [float(x) for x in zeros if x>0]
366
367
# start at X=2
368
cdef double s, gamma, logX=log(2), logXmax = log(Xmax)
369
cdef double delta = logXmax / n
370
371
while logX <= logXmax:
372
s = 0
373
for gamma in zeros:
374
s += sin(gamma*logX)/gamma
375
v.append((logX, 2*s))
376
logX += delta
377
378
return v
379