# williamstein/mazur-explicit-formula

Added a bunch of functions for talk.

1 parent 939e8c0 commit 2d2b510d53c6e72d1ad18762f09afdee75a9aff3 committed Mar 2, 2013
Showing with 124 additions and 23 deletions.
1. +23 −1 explicit.pyx
2. +78 −4 go2.sage
 @@ -7,7 +7,7 @@ import math, os, sys -from sage.all import prime_range, line, tmp_dir, parallel, text, cached_method, nth_prime, load, EllipticCurve, walltime, TimeSeries +from sage.all import prime_range, line, tmp_dir, parallel, text, cached_method, nth_prime, load, EllipticCurve, walltime, TimeSeries, animate from scipy.special import expi as Ei @@ -243,6 +243,28 @@ class OscillatoryTerm(object): v = self(X) return line(enumerate(v), **kwds) + def animation(self, Xvals, output_path=None, ncpus=1, **kwds): + if not isinstance(Xvals, list): + raise TypeError, "Xvals must be a list" + if output_path is not None and not os.path.exists(output_path): + os.makedirs(output_path) + + print("Rendering %s frames to %s using %s cpus"%(len(Xvals), output_path, ncpus)) + + @parallel(ncpus) + def f(X): + return self(X) + + V = list(f(Xvals)) + ymax = max([max(v[1]) for v in V]) + ymin = min([min(v[1]) for v in V]) + frames = [] + for X, v in V: + frame = ( line(enumerate(v), ymax=ymax, ymin=ymin, **kwds) + + text("X = %s"%X[0], (len(v)//6, ymax/2.0), color='black', fontsize=16) ) + frames.append(frame) + return animate(frames) + def animation_svg(self, Xvals, output_path, ncpus=1, **kwds): if not isinstance(Xvals, list): raise TypeError, "Xvals must be a list"
 @@ -378,7 +378,7 @@ def sato_tate_animation(lbl, frames=50, B=10^9, **kwds): ##################################### # 2. Chebeshev Bias ##################################### -def chebeshev_distribution_animation(E, frames=50, B=10^9, **kwds): +def chebeshev_distribution_animation(frames=50, B=10^9, **kwds): v = chebeshev_data(B) return animated_histogram(v, frames, **kwds) @@ -419,10 +419,84 @@ class Madison: target = 'sato-tate-animation-%s.gif'%lbl if self.done(target): return - sato_tate_animation(lbl, frames=10, B=10^8, bins=10000).save(self.path(target)) + sato_tate_animation(lbl, frames=100, B=10^9, bins=10000, ymax=1).save(self.path(target)) + + def primes_mod_4(self): + target = "primes-mod-4.gif" + if self.done(target): + return + chebeshev_distribution_animation(frames=100, B=10^9, bins=10).save(self.path(target)) + + def prime_race_mod4(self): + target = "primes-race-mod4.svg" + if self.done(target): + return + v = chebeshev_bias_data(10^9) + v.plot(thickness=0.3, plot_points=10^6).save(self.path(target), figsize=[8,3]) + + def riemann_oscillatory(self): + target = "riemann-oscillatory.svg" + if self.done(target): return + g = plot(lambda x: zeta_oscillatory(x, 0,10000), (3, 1000), plot_points=10000) + g.save(self.path(target)) + + def riemann_epsilon(self): + target = "riemann-epsilon.svg" + if self.done(target): return + v = zeta_eps(1000, 10000) + line(v, thickness=.5).save(self.path(target),figsize=[8,3]) + + def data_plots(self, lbl, B=10^9, force=False): + targets = ["raw-data-%s.svg"%lbl, "raw-mean-%s.svg"%lbl, + "medium-data-%s.svg"%lbl, "medium-mean-%s.svg"%lbl, + "well-data-%s.svg"%lbl, "well-mean-%s.svg"%lbl] + if not force and all(self.done(t) for t in targets): + return + dp = DataPlots(lbl, B, data_path='data') + v = dp.data(num_points=5000, log_X_scale=True) + for w in ['raw','medium','well']: + g = plot_step_function(v[w]['delta'],thickness=1,fontsize=18) + g.save('madison/%s-data-%s.svg'%(w,lbl), gridlines=True, figsize=[10,4]) + g = plot_step_function(v[w]['mean'],thickness=1,fontsize=18) + g.save('madison/%s-mean-%s.svg'%(w,lbl), gridlines=True, figsize=[10,4]) + + def oscillatory(self, lbl, B=1e30): + target = "oscillatory-%s.gif"%lbl + if self.done(target): return + z = zeros(lbl) + v = [line(zero_sum_no_log_plot(z[:k], 10000, B), thickness=.4, figsize=[10,4]) for k in [5,50,..,500]] + ymax = max([g.ymax() for g in v]) + ymin = min([g.ymin() for g in v]) + a = animate(v, ymin=ymin, ymax=ymax) + a.save(self.path(target)) + + def bayesian(self, lbl, B=1e30): + target = "bayesian-%s.svg"%lbl + if self.done(target): return + z = zeros(lbl) + t = TimeSeries(zero_sum_distribution1(z, 100000, math.log(B))) + g = t.plot_histogram(bins=1000) + mean = t.mean(); sd = t.standard_deviation() + pdf = lambda x: 1/(sd*sqrt(2*pi)) * exp(-(x-mean)^2/(2*sd^2)) + key = "%s: sd=%.2f, mean=%.2f"%(lbl, sd, mean) + g += text(key, (-4*sd,1), color='black') + g += plot(pdf, (x,mean-4*sd,mean+4*sd), color='red', thickness=2) + g.save(self.path(target), figsize=[10,6]) + + def madison(): m = Madison() m.sato_tate('11a') - - + m.sato_tate('5077a') + m.primes_mod_4() + m.prime_race_mod4() + m.riemann_oscillatory() + m.riemann_epsilon() + m.data_plots('11a') + m.data_plots('5077a') + m.data_plots('2379b') + m.data_plots('128b') + m.oscillatory('11a') + m.oscillatory('128b') + m.oscillatory('5077a')
 @@ -32,7 +32,7 @@

1. Birch and Swinnerton-Dyer's Hunch and Sato-Tate

Sato-Tate for 11a (rank 0) for $p<10^9$

-

Sato-Tate for 5077a (rank 3) for $p^10^9$

+

Sato-Tate for 5077a (rank 3) for $p<10^9$

@@ -43,14 +43,11 @@

Distribution of primes mod 4 for $p < 10^9$

Prime race mod 4:   $\sum_{2 < p\leq X} \left(\frac{-1}{p}\right)$ for $p < 10^9$

- +

3. Riemann's Explicit Formula

-

Von Mangoldt for $p<10^9$

- -

Oscillatory Term

@@ -86,7 +83,6 @@

5. Sarnak's Letter - Conjectures

Raw Mean on a Log Scale for 11a (rank 0) for $p<10^9$

-

Raw Mean on a Log Scale for 5077a (rank 3) for $p<10^9$

@@ -101,33 +97,42 @@

Well Done Mean on a Log Scale for 11a (rank 0) for $p<10^9$

Well Done Mean on a Log Scale for 5077a (rank 3) for $p<10^9$

+
-

6. A Minimalist Conjecture

+

6. Minimalist Conjecture?

+ + + +

Raw Mean on a Log Scale for 128b (rank 0, but Sym^{7} extra vanishing) for $p<10^9$

+ + +

Raw Mean on a Log Scale for 2379b (rank 0, but Sym^{3} extra vanishing) for $p<10^9$

+ + -

Raw Mean on a Log Scale for 11a (rank 0) for $p<10^9$ with Minimalist Line

- -

Raw Mean on a Log Scale for 128b (rank 0) for $p<10^9$ with Minimalist Line

-(but here there is a nontrivial higher zero!) -
- -

7. Oscillatory Term -- the Bayesian

-

Plot of Oscillatory Term S(X,T) for $X<10^{100}$ (log scale) and with $T$ changing from small to $10000$ (animated) for 11a (rank 0)

+

Plot of Oscillatory Term S(X,T)/log(X) for $X<10^{30}$ (log scale) and with $T$ changing from small to $500$ (animated) for 11a (rank 0)

-

Plot of Distribution of Values of above for $T=10000$ for 11a (rank 0)

+

Plot of Oscillatory Term S(X,T)/log(X) for $X<10^{30}$ (log scale) and with $T$ changing from small to $500$ (animated) for 128b (rank 0, sym pow)

+ + +

Plot of Oscillatory Term S(X,T)/log(X) for $X<10^{30}$ (log scale) and with $T$ changing from small to $500$ (animated) for 5077a (rank 3)

+ + +

Plot of Distribution of Values of S(X,10000) for 11a (rank 0)

-

Plot of Distribution of Values of above for $T=10000$ for 128b (rank 0)

+

Plot of Distribution of Values of S(X,10000) for 128b (rank 0)

-

Plot of Distribution of Values of above for $T=10000$ for 5077a (rank 3)

+

Plot of Distribution of Values of S(X,10000) for 5077a (rank 3)