CoCalc Public Fileswww / talks / mazur-explicit-formula / explicit.pyx
Author: William A. Stein
1#########################################################
2#
3# How Explicit is the Explict Formula
4#
5#   (c) William Stein, 2013
6#########################################################
7
8import math, os, sys
9
10from sage.all import prime_range, line, tmp_dir, parallel, text, cached_method, nth_prime, load, EllipticCurve, walltime
11
12from scipy.special import expi as Ei
13
14cdef extern from "math.h":
15    double log(double)
16    double sqrt(double)
17    double cos(double)
18    double sin(double)
19
20cdef double pi = math.pi
21
22################################################
23# The raw data -- means
24################################################
25"""
26sage: dp = DataPlots("160a", 10^6)
27sage: v = dp.data(num_points=5000)
28sage: plot_step_function(v['raw']['mean'],color='red',thickness=.5) + plot_step_function(v['raw']['delta'], thickness=.5)
29"""
30
31def 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
162def 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
174def 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
191class 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:
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
216class 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
300def 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
326def 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
353def 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