 CoCalc Public Files2016-02-23-psi-formula.sagews
Author: Nick Bohall
Views : 58
Compute Environment: Ubuntu 18.04 (Deprecated)
%md










︠%auto

a = 5


a

5   def fig_riemann_spectrum_gaps():
t = stats.TimeSeries(zeta_zeros())
g = t.diffs().plot_histogram(bins=500, normalize=False)
g.show()

fig_riemann_spectrum_gaps() def fig_mini_phihat_even():
G = plot_symbolic_phihat(5, 1, 100, 10000, zeros=True)
G.show()

def fig_phihat_even():
for bound in [5, 10, 20]:
G = plot_symbolic_phihat(bound, 2, 100,
plot_points=10^5)
G.show()

def fig_phihat_even_all():
p = [plot_symbolic_phihat(n, 2,100) for n in [5,20,50,500]]
[a.ymin(0) for a in p]
for a in p:
a.show()
g = graphics_array([[a] for a in p],4,1)
g.show()

def symbolic_phihat(bound):
t = var('t')
f = SR(0)#what is this - all expressions learned in calculus
for pn in prime_powers(bound+1):
if pn == 1: continue
p, e = factor(pn)
f += - log(p)/sqrt(pn) * cos(t*log(pn))
return f

#omit 2 and - signs for p = 3 (mod 4)
#do Fourier trans on sum in phi
def plot_symbolic_phihat(bound, xmin, xmax, plot_points=1000, zeros=True):
f = symbolic_phihat(bound)
P = plot(f, (t,xmin, xmax), plot_points=plot_points)
if not zeros:
return P
ym = P.ymax()
Z = []
for y in zeta_zeros():
if y > xmax: break
Z.append(y)
zeros = sum([arrow((x,ym),(x,0),rgbcolor='red',width=0.5,arrowsize=2)
for i, x in enumerate(Z)])
return P + zeros
fig_mini_phihat_even()
#fig_phihat_even()

Error in lines 1-4 Traceback (most recent call last): File "/projects/sage/sage-6.10/local/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 905, in execute exec compile(block+'\n', '', 'single') in namespace, locals File "", line 1, in <module> File "sage/structure/sage_object.pyx", line 1097, in sage.structure.sage_object.save (/projects/sage/sage-6.10/src/build/cythonized/sage/structure/sage_object.c:12013) if not os.path.splitext(filename) or not hasattr(obj, "save"): File "/projects/sage/sage-6.10/local/lib/python/posixpath.py", line 98, in splitext return genericpath._splitext(p, sep, altsep, extsep) File "/projects/sage/sage-6.10/local/lib/python/genericpath.py", line 99, in _splitext sepIndex = p.rfind(sep) AttributeError: 'NoneType' object has no attribute 'rfind'
@auto
def gauss_primepowers(bound):
"""
Computes all powers of Gaussian primes less than a given norm.
Returns them as a list
"""
def pgen(I):
I.is_principal()
assert len(I._reduced_generators) == 1
return I._reduced_generators

def exponentiate(base, power):
res = base
while power > 1:
res = res * base
power = power - 1
return res

K.<i> = QQ[sqrt(-1)]
eles = K.primes_of_bounded_norm(bound-1)
print "eles: ",
print eles
primes = [pgen(I) for I in eles]
print "primes: ",
print primes
primepowers = []
for p in primes:
primepowers.append(p)
n = 2
p = exponentiate(p,n)
while p.norm() < bound:
primepowers.append(p)
n = n + 1
p = exponentiate(p,n)
return primepowers

def zeros(y):
s = os.popen('lcalc --zeros-interval -x 1 -y %s --stepsize .01 --twist-quadratic -s -4 -f -4'%y).readlines()
return [float(z.split()) for z in s]

#omit 2 and - signs for p = 3 (mod 4)
#do Fourier trans on sum in phi
def gauss_symbolic_phihat(bound):
t = var('t')
f = SR(0)#what is this - all expressions learned in calculus
for pn in prime_powers(bound+1):
if pn == 1: continue
if pn % 2 != 0:
p, e = factor(pn)
if p % 4 == 3:
f += log(p)/sqrt(pn) * cos(t*log(pn))
else:
f += -log(p)/sqrt(pn) * cos(t*log(pn))
return f

def gauss_plot_symbolic_phihat(bound, xmin, xmax, plot_points=1000, show_zeros=True):
f = gauss_symbolic_phihat(bound)
P = plot(f, (t,xmin, xmax), plot_points=plot_points)
if not show_zeros:
return P
ym = P.ymax()
Z = []
for y in zeros(xmax):
Z.append(y)
zeros_arrows = sum([arrow((x,ym),(x,0),rgbcolor='red',width=0.5,arrowsize=2)
for i, x in enumerate(Z)])
return P + zeros_arrows

def gauss_fig_mini_phihat_even():
G = gauss_plot_symbolic_phihat(5, 1, 100, 10000)
G.show()

plot_symbolic_phihat(5, 2,100) %time
frames = []
for n in [6,7,8]:
print n
p = gauss_plot_symbolic_phihat(n, 2,100)
p.xmin(0); p.xmax(100); p.ymin(-2); p.ymax(2)
frames.append(p)

6 7 8 CPU time: 12.98 s, Wall time: 13.63 s
animate(frames).save('a.gif', delay=50)


for n in [5, 20, 50]:
g = gauss_plot_symbolic_phihat(n, 2,100)
g.show()   ︠

zeros(50)

[6.0209489046976, 10.243770304167, 12.988098012312, 16.342607104587, 18.291993196124, 21.450611343983, 23.27837652046, 25.728756425089, 28.359634343025, 29.656384014593, 32.592186527117, 34.199957509213, 36.142880458303, 38.511923141719, 40.322674066691, 41.807084620005, 44.617891058662, 45.599584396792, 47.741562280939, 49.723129323783]
def lbl(name, f):
open("thesis/graphics/%s.tex"%name, 'w').write("$$f(t)=%s$$"%latex(f))
f = gauss_symbolic_phihat(3)
lbl("intro1", f)
p = f.plot(0,20)
save(p, 'thesis/graphics/intro1.pdf')



f = gauss_symbolic_phihat(5); show(f)
lbl('intro2', f)
p = f.plot(0, 20)
p
save(p, 'thesis/graphics/intro2.pdf')


$\displaystyle -\frac{1}{5} \, \sqrt{5} \cos\left(t \log\left(5\right)\right) \log\left(5\right) + \frac{1}{3} \, \sqrt{3} \cos\left(t \log\left(3\right)\right) \log\left(3\right)$ f = gauss_symbolic_phihat(7)
lbl("intro3", f)
p = f.plot(0, 20)
save(p, 'thesis/graphics/intro3.pdf')

f = gauss_symbolic_phihat(11); show(f)
f.plot(0, 20)

$\displaystyle \frac{1}{11} \, \sqrt{11} \cos\left(t \log\left(11\right)\right) \log\left(11\right) + \frac{1}{7} \, \sqrt{7} \cos\left(t \log\left(7\right)\right) \log\left(7\right) - \frac{1}{5} \, \sqrt{5} \cos\left(t \log\left(5\right)\right) \log\left(5\right) + \frac{1}{3} \, \sqrt{3} \cos\left(t \log\left(3\right)\right) \log\left(3\right) + \frac{1}{3} \, \cos\left(t \log\left(9\right)\right) \log\left(3\right)$ f = gauss_symbolic_phihat(13); show(f)
f.plot(0, 20)


$\displaystyle -\frac{1}{13} \, \sqrt{13} \cos\left(t \log\left(13\right)\right) \log\left(13\right) + \frac{1}{11} \, \sqrt{11} \cos\left(t \log\left(11\right)\right) \log\left(11\right) + \frac{1}{7} \, \sqrt{7} \cos\left(t \log\left(7\right)\right) \log\left(7\right) - \frac{1}{5} \, \sqrt{5} \cos\left(t \log\left(5\right)\right) \log\left(5\right) + \frac{1}{3} \, \sqrt{3} \cos\left(t \log\left(3\right)\right) \log\left(3\right) + \frac{1}{3} \, \cos\left(t \log\left(9\right)\right) \log\left(3\right)$ f = gauss_symbolic_phihat(20); show(f)
f.plot(0, 20)

$\displaystyle \frac{1}{19} \, \sqrt{19} \cos\left(t \log\left(19\right)\right) \log\left(19\right) - \frac{1}{17} \, \sqrt{17} \cos\left(t \log\left(17\right)\right) \log\left(17\right) - \frac{1}{13} \, \sqrt{13} \cos\left(t \log\left(13\right)\right) \log\left(13\right) + \frac{1}{11} \, \sqrt{11} \cos\left(t \log\left(11\right)\right) \log\left(11\right) + \frac{1}{7} \, \sqrt{7} \cos\left(t \log\left(7\right)\right) \log\left(7\right) - \frac{1}{5} \, \sqrt{5} \cos\left(t \log\left(5\right)\right) \log\left(5\right) + \frac{1}{3} \, \sqrt{3} \cos\left(t \log\left(3\right)\right) \log\left(3\right) + \frac{1}{3} \, \cos\left(t \log\left(9\right)\right) \log\left(3\right)$ %time
f = gauss_symbolic_phihat(200); show(f)
f.plot(0, 20, plot_points=100)

$\displaystyle \frac{1}{199} \, \sqrt{199} \cos\left(t \log\left(199\right)\right) \log\left(199\right) - \frac{1}{197} \, \sqrt{197} \cos\left(t \log\left(197\right)\right) \log\left(197\right) - \frac{1}{193} \, \sqrt{193} \cos\left(t \log\left(193\right)\right) \log\left(193\right) + \frac{1}{191} \, \sqrt{191} \cos\left(t \log\left(191\right)\right) \log\left(191\right) - \frac{1}{181} \, \sqrt{181} \cos\left(t \log\left(181\right)\right) \log\left(181\right) + \frac{1}{179} \, \sqrt{179} \cos\left(t \log\left(179\right)\right) \log\left(179\right) - \frac{1}{173} \, \sqrt{173} \cos\left(t \log\left(173\right)\right) \log\left(173\right) + \frac{1}{167} \, \sqrt{167} \cos\left(t \log\left(167\right)\right) \log\left(167\right) + \frac{1}{163} \, \sqrt{163} \cos\left(t \log\left(163\right)\right) \log\left(163\right) - \frac{1}{157} \, \sqrt{157} \cos\left(t \log\left(157\right)\right) \log\left(157\right) + \frac{1}{151} \, \sqrt{151} \cos\left(t \log\left(151\right)\right) \log\left(151\right) - \frac{1}{149} \, \sqrt{149} \cos\left(t \log\left(149\right)\right) \log\left(149\right) + \frac{1}{139} \, \sqrt{139} \cos\left(t \log\left(139\right)\right) \log\left(139\right) - \frac{1}{137} \, \sqrt{137} \cos\left(t \log\left(137\right)\right) \log\left(137\right) + \frac{1}{131} \, \sqrt{131} \cos\left(t \log\left(131\right)\right) \log\left(131\right) + \frac{1}{127} \, \sqrt{127} \cos\left(t \log\left(127\right)\right) \log\left(127\right) - \frac{1}{113} \, \sqrt{113} \cos\left(t \log\left(113\right)\right) \log\left(113\right) - \frac{1}{109} \, \sqrt{109} \cos\left(t \log\left(109\right)\right) \log\left(109\right) + \frac{1}{107} \, \sqrt{107} \cos\left(t \log\left(107\right)\right) \log\left(107\right) + \frac{1}{103} \, \sqrt{103} \cos\left(t \log\left(103\right)\right) \log\left(103\right) - \frac{1}{101} \, \sqrt{101} \cos\left(t \log\left(101\right)\right) \log\left(101\right) - \frac{1}{97} \, \sqrt{97} \cos\left(t \log\left(97\right)\right) \log\left(97\right) - \frac{1}{89} \, \sqrt{89} \cos\left(t \log\left(89\right)\right) \log\left(89\right) + \frac{1}{83} \, \sqrt{83} \cos\left(t \log\left(83\right)\right) \log\left(83\right) + \frac{1}{79} \, \sqrt{79} \cos\left(t \log\left(79\right)\right) \log\left(79\right) - \frac{1}{73} \, \sqrt{73} \cos\left(t \log\left(73\right)\right) \log\left(73\right) + \frac{1}{71} \, \sqrt{71} \cos\left(t \log\left(71\right)\right) \log\left(71\right) + \frac{1}{67} \, \sqrt{67} \cos\left(t \log\left(67\right)\right) \log\left(67\right) - \frac{1}{61} \, \sqrt{61} \cos\left(t \log\left(61\right)\right) \log\left(61\right) + \frac{1}{59} \, \sqrt{59} \cos\left(t \log\left(59\right)\right) \log\left(59\right) - \frac{1}{53} \, \sqrt{53} \cos\left(t \log\left(53\right)\right) \log\left(53\right) + \frac{1}{47} \, \sqrt{47} \cos\left(t \log\left(47\right)\right) \log\left(47\right) + \frac{1}{43} \, \sqrt{43} \cos\left(t \log\left(43\right)\right) \log\left(43\right) - \frac{1}{41} \, \sqrt{41} \cos\left(t \log\left(41\right)\right) \log\left(41\right) - \frac{1}{37} \, \sqrt{37} \cos\left(t \log\left(37\right)\right) \log\left(37\right) + \frac{1}{31} \, \sqrt{31} \cos\left(t \log\left(31\right)\right) \log\left(31\right) - \frac{1}{29} \, \sqrt{29} \cos\left(t \log\left(29\right)\right) \log\left(29\right) + \frac{1}{23} \, \sqrt{23} \cos\left(t \log\left(23\right)\right) \log\left(23\right) + \frac{1}{19} \, \sqrt{19} \cos\left(t \log\left(19\right)\right) \log\left(19\right) - \frac{1}{17} \, \sqrt{17} \cos\left(t \log\left(17\right)\right) \log\left(17\right) - \frac{1}{13} \, \sqrt{13} \cos\left(t \log\left(13\right)\right) \log\left(13\right) + \frac{1}{11} \, \sqrt{11} \cos\left(t \log\left(11\right)\right) \log\left(11\right) + \frac{1}{7} \, \sqrt{7} \cos\left(t \log\left(7\right)\right) \log\left(7\right) - \frac{1}{25} \, \sqrt{5} \cos\left(t \log\left(125\right)\right) \log\left(5\right) - \frac{1}{5} \, \sqrt{5} \cos\left(t \log\left(5\right)\right) \log\left(5\right) + \frac{1}{9} \, \sqrt{3} \cos\left(t \log\left(27\right)\right) \log\left(3\right) + \frac{1}{3} \, \sqrt{3} \cos\left(t \log\left(3\right)\right) \log\left(3\right) - \frac{1}{13} \, \cos\left(t \log\left(169\right)\right) \log\left(13\right) + \frac{1}{11} \, \cos\left(t \log\left(121\right)\right) \log\left(11\right) + \frac{1}{7} \, \cos\left(t \log\left(49\right)\right) \log\left(7\right) - \frac{1}{5} \, \cos\left(t \log\left(25\right)\right) \log\left(5\right) + \frac{1}{9} \, \cos\left(t \log\left(81\right)\right) \log\left(3\right) + \frac{1}{3} \, \cos\left(t \log\left(9\right)\right) \log\left(3\right)$ CPU time: 33.98 s, Wall time: 34.76 s
zeros(100)

[6.0209489046976, 10.243770304167, 12.988098012312, 16.342607104587, 18.291993196124, 21.450611343983, 23.27837652046, 25.728756425089, 28.359634343025, 29.656384014593, 32.592186527117, 34.199957509213, 36.142880458303, 38.511923141719, 40.322674066691, 41.807084620005, 44.617891058662, 45.599584396792, 47.741562280939, 49.723129323783, 51.686093452871, 52.768820767805, 55.267543584699, 56.934374055202, 58.116707110674, 60.421713949008, 62.008632285768, 63.714641118785, 64.976170573096, 67.636920863546, 68.365884503834, 70.185879908802, 72.155484974382, 73.767635521486, 75.143121647433, 76.69630320343, 78.809998314321, 80.210131238367, 81.213951626883, 83.666656014471, 84.731740363782, 86.57766016839, 87.629718119588, 89.801131616696, 91.349703814698, 92.237499910454, 94.16661958596, 96.136011161781, 96.961741579417, 98.755300415755]
g = f.derivative()

f.find_local_maximum(5.5,6.5)

(4.625184318694296, 6.0932418136308746)
%time f = gauss_symbolic_phihat(20000)
f.find_local_maximum(5.5,6.5)

CPU time: 1.40 s, Wall time: 1.42 s (9.609716603100983, 6.0190156444824714)
%time f = gauss_symbolic_phihat(100000)
f.find_local_maximum(5.5,6.5)

CPU time: 9.16 s, Wall time: 9.49 s (10.306059947208393, 6.0248752427482559)
g = f._fast_float_()

/projects/sage/sage-6.10/local/lib/python2.7/site-packages/sage/symbolic/expression_conversions.py:1387: DeprecationWarning: Substitution using function-call syntax and unnamed arguments is deprecated and will be removed from a future release of Sage; you can use named arguments instead, like EXPR(x=..., y=...) See http://trac.sagemath.org/5930 for details. return FastFloatConverter(ex, *vars)()
%time plot(g, 5, 100) CPU time: 8.78 s, Wall time: 8.88 s
p = 1e9
N(log(p)/sqrt(p))

0.000655327206019062

@cached_function
def zeros(y, ring):
s = os.popen('lcalc --zeros-interval -x 1 -y %s --stepsize .01 --twist-quadratic -s %s -f %s'%(y, ring.discriminant(),
return [float(z.split()) for z in s]

#omit 2 and - signs for p = 3 (mod 4)
#do Fourier trans on sum in phi
def arb_symbolic_phihat(bound,delta):
t = var('t')
f = SR(0)#what is this - all expressions learned in calculus
for pn in prime_powers(bound+1):
if pn == 1: continue
if pn % 2 != 0:
p, e = factor(pn)
#print "Kronecker for prime %s of delta %s is %s" % (p, delta, kronecker(delta, p))
f += -kronecker(delta,p) * log(p)/sqrt(pn) * cos(t*log(pn))
return f

def arb_plot_symbolic_phihat(bound, xmin, xmax, delta = -1, plot_points=1000, show_zeros=True):
f = arb_symbolic_phihat(bound, delta)
K.<i> = QQ[sqrt(delta)]
P = plot(f, (t,xmin, xmax), plot_points=plot_points)
if not show_zeros:
return P
ym = P.ymax()
Z = []
for y in zeros(xmax, K):
Z.append(y)
zeros_arrows = sum([arrow((x,ym),(x,0),rgbcolor='red',width=0.5,arrowsize=2)
for i, x in enumerate(Z)])
return P + zeros_arrows

p = arb_plot_symbolic_phihat(100, 2, 30, -5)
p.show() g = gauss_plot_symbolic_phihat(100, 2, 30)
g.show() K.<i> = QQ[sqrt(-1)]
K.discriminant??