CoCalc Public Filesbook / rh / code / code.sageOpen with one click!
Author: William A. Stein
Views : 65
Compute Environment: Ubuntu 18.04 (Deprecated)
1
2
# Options
3
4
##############################################################
5
# Drawing figures (one or all)
6
##############################################################
7
8
def draw(fig=None, dir='illustrations/',ext='pdf'):
9
if isinstance(fig, str):
10
print "Drawing %s... "%fig,
11
sys.stdout.flush()
12
t = walltime()
13
eval('fig_%s(dir,ext)'%fig)
14
print " (time = %s seconds)"%walltime(t)
15
return
16
17
if fig is None:
18
figs = ['_'.join(x.split('_')[1:]) for x in globals() if x.startswith('fig_')]
19
elif isinstance(fig, list):
20
figs = fig
21
for fig in figs:
22
draw(fig)
23
return
24
25
##############################################################
26
# Factorization trees
27
##############################################################
28
def fig_factor_tree(dir, ext):
29
g = FactorTree(6).plot(labels={'fontsize':60},sort=True)
30
g.save(dir + '/factor_tree_6.%s'%ext, axes=False, axes_pad=0.1, typeset='latex')
31
32
g = FactorTree(12).plot(labels={'fontsize':50},sort=True)
33
g.save(dir + '/factor_tree_12.%s'%ext, axes=False, axes_pad=0.1, typeset='latex')
34
35
set_random_seed(3)
36
g = FactorTree(12).plot(labels={'fontsize':50},sort=False)
37
g.save(dir + '/factor_tree_12b.%s'%ext, axes=False,axes_pad=0.1, typeset='latex')
38
39
set_random_seed(0)
40
for w in ['a', 'b']:
41
g = FactorTree(300).plot(labels={'fontsize':40},sort=False)
42
g.save(dir + '/factor_tree_300_%s.%s'%(w,ext), axes=False, axes_pad=0.1, typeset='latex')
43
44
set_random_seed(0)
45
g = FactorTree(6469693230).plot(labels={'fontsize':14},sort=False)
46
g.save(dir + '/factor_tree_big.%s'%ext, axes=False, axes_pad=0.1, typeset='latex')
47
48
49
class FactorTree:
50
"""
51
A factorization tree.
52
53
EXAMPLES::
54
55
sage: FactorTree(100)
56
Factorization trees of 100
57
sage: R.<x> = QQ[]
58
sage: FactorTree(x^3-1)
59
Factorization trees of x^3 - 1
60
"""
61
def __init__(self, n):
62
"""
63
INPUT:
64
65
- `n` -- number of element of polynomial ring
66
"""
67
self.__n = n
68
69
def value(self):
70
"""
71
Return the n such that self is FactorTree(n).
72
73
EXAMPLES::
74
75
sage: FactorTree(100).value()
76
100
77
"""
78
return self.__n
79
80
def __repr__(self):
81
"""
82
EXAMPLES::
83
84
sage: FactorTree(100).__repr__()
85
'Factorization trees of 100'
86
"""
87
return "Factorization trees of %s"%self.__n
88
89
def plot(self, lines=None, labels=None, sort=False):
90
"""
91
INPUT:
92
93
- ``lines`` -- optional dictionary of options passed
94
to the line command
95
96
- ``labels`` -- optional dictionary of options passed to
97
the text command for the divisor labels
98
99
- ``sort`` -- bool (default: False); if True, the primes
100
divisors are found in order from smallest to largest;
101
otherwise, the factor tree is draw at random
102
103
EXAMPLES::
104
105
sage: FactorTree(2009).plot(labels={'fontsize':30},sort=True)
106
107
We can make factor trees of polynomials in addition to integers::
108
109
sage: R.<x> = QQ[]
110
sage: F = FactorTree((x^2-1)*(x^3+2)*(x-5)); F
111
Factorization trees of x^6 - 5*x^5 - x^4 + 7*x^3 - 10*x^2 - 2*x + 10
112
sage: F.plot(labels={'fontsize':15},sort=True)
113
"""
114
if lines is None:
115
lines = {'rgbcolor':(.5,.5,1)}
116
if labels is None:
117
labels = {'fontsize':16}
118
n = self.__n
119
rows = []
120
v = [(n,None,0)]
121
self._ftree(rows, v, sort=sort)
122
return self._plot_ftree(rows, lines, labels)
123
124
def _ftree(self, rows, v, sort):
125
"""
126
A function that is used recurssively internally by the plot function.
127
128
INPUT:
129
130
- ``rows`` -- list of lists of integers
131
132
- ``v`` -- list of triples of integers
133
134
- ``sort`` -- bool (default: False); if True, the primes
135
divisors are found in order from smallest to largest;
136
otherwise, the factor tree is draw at random
137
138
EXAMPLES::
139
140
sage: F = FactorTree(100)
141
sage: rows = []; v = [(100,None,0)]; F._ftree(rows, v, True)
142
sage: rows
143
[[(100, None, 0)],
144
[(2, 100, 0), (50, 100, 0)],
145
[(None, None, None), (2, 50, 1), (25, 50, 1)],
146
[(None, None, None), (None, None, None), (5, 25, 2), (5, 25, 2)]]
147
sage: v
148
[(100, None, 0)]
149
"""
150
if len(v) > 0:
151
# add a row to g at the ith level.
152
rows.append(v)
153
w = []
154
for i in range(len(v)):
155
k, _,_ = v[i]
156
if k is None or k.is_irreducible():
157
w.append((None,None,None))
158
else:
159
div = divisors(k)
160
if sort:
161
d = div[1]
162
else:
163
z = divisors(k)[1:-1]
164
d = z[randint(0,len(z)-1)]
165
w.append((d,k,i))
166
e = k//d
167
if e == 1:
168
w.append((None,None))
169
else:
170
w.append((e,k,i))
171
if len(w) > len(v):
172
self._ftree(rows, w, sort=sort)
173
174
def _plot_ftree(self, rows, lines, labels):
175
"""
176
Used internally by the plot method.
177
178
INPUT:
179
180
- ``rows`` -- list of lists of triples
181
182
- ``lines`` -- dictionary of options to pass to lines commands
183
184
- ``labels`` -- dictionary of options to pass to text command to label factors
185
186
EXAMPLES::
187
188
sage: F = FactorTree(100)
189
sage: rows = []; v = [(100,None,0)]; F._ftree(rows, v, True)
190
sage: F._plot_ftree(rows, {}, {})
191
"""
192
g = Graphics()
193
for i in range(len(rows)):
194
cur = rows[i]
195
for j in range(len(cur)):
196
e, f, k = cur[j]
197
if not e is None:
198
if e.is_irreducible():
199
c = (1r,0r,0r)
200
else:
201
c = (0r,0r,.4r)
202
g += text("$%s$"%latex(e), (j*2-len(cur),-i), rgbcolor=c, **labels)
203
if not k is None and not f is None:
204
g += line([(j*2-len(cur),-i), (k*2-len(rows[i-1]),-i+1)], axes=False, **lines)
205
return g
206
207
208
##############################################################
209
# Bag of primes
210
##############################################################
211
212
def bag_of_primes(steps):
213
"""
214
This works up to 9. It took a day using specialized factoring
215
(via GMP-ECM) to get step 10.
216
217
EXAMPLES::
218
219
sage: bag_of_primes(5)
220
[2, 3, 7, 43, 13, 139, 3263443]
221
"""
222
bag = [2]
223
for i in range(steps):
224
for p in prime_divisors(prod(bag)+1):
225
bag.append(p)
226
print bag
227
228
229
##############################################################
230
# Questions about numbers
231
##############################################################
232
def fig_questions(dir, ext):
233
g = questions(100,17,20)
234
g.save(dir + '/questions.%s'%ext, axes=False, typeset='latex')
235
236
def questions(n=100,k=17,fs=20):
237
set_random_seed(k)
238
g = text("?",(5,5),rgbcolor='grey', fontsize=200)
239
g += sum(text("$%s$"%p,(random()*10,random()*10),rgbcolor=(p/(2*n),p/(2*n),p/(2*n)),fontsize=fs)
240
for p in primes(n))
241
return g
242
243
244
##############################################################
245
# Sieve of Eratosthenes
246
##############################################################
247
248
def fig_erat(dir,ext):
249
# sieving out 2,3,5,7
250
for p in [2,3,5,7]:
251
sieve_step(p,100).save(dir+'/sieve100-%s.%s'%(p,ext), typeset='latex')
252
253
sieve_step(13,200).save(dir+'/sieve200.%s'%ext, typeset='latex')
254
255
256
def number_grid(c, box_args=None, font_args=None, offset=0):
257
"""
258
INPUT:
259
c -- list of length n^2, where n is an integer.
260
The entries of c are RGB colors.
261
box_args -- additional arguments passed to box.
262
font_args -- all additional arguments are passed
263
to the text function, e.g., fontsize.
264
offset -- use to fine tune text placement in the squares
265
266
OUTPUT:
267
Graphics -- a plot of a grid that illustrates
268
those n^2 numbers colored according to c.
269
"""
270
if box_args is None: box_args = {}
271
if font_args is None: font_args = {}
272
273
try:
274
n = sqrt(ZZ(len(c)))
275
except ArithmeticError:
276
raise ArithmeticError, "c must have square length"
277
G = Graphics()
278
k = 0
279
for j in reversed(range(n)):
280
for i in range(n):
281
col = c[int(k)]
282
R = line([(i,j),(i+1,j),(i+1,j+1),(i,j+1),(i,j)],
283
thickness=.2, **box_args)
284
d = dict(box_args)
285
if 'rgbcolor' in d.keys():
286
del d['rgbcolor']
287
P = polygon([(i,j),(i+1,j),(i+1,j+1),(i,j+1),(i,j)],
288
rgbcolor=col, **d)
289
G += P + R
290
if col != (1,1,1):
291
G += text(str(k+1), (i+.5+offset,j+.5), **font_args)
292
k += 1
293
G.axes(False)
294
G.xmin(0); G.xmax(n); G.ymin(0); G.ymax(n)
295
G.set_aspect_ratio('automatic')
296
return G
297
298
def sieve_step(p, n, gone=(1,1,1), prime=(1,0,0), \
299
multiple=(.6,.6,.6), remaining=(.9,.9,.9),
300
fontsize=11,offset=0):
301
"""
302
Return graphics that illustrates sieving out multiples of p.
303
Numbers that are a nontrivial multiple of primes < p are shown in
304
the gone color. Numbers that are a multiple of p are shown in a
305
different color. The number p is shown in yet a third color.
306
307
INPUT:
308
p -- a prime (or 1, in which case all points are colored "remaining")
309
n -- a SQUARE integer
310
gone -- rgb color for integers that have been sieved out
311
prime -- rgb color for p
312
multiple -- rgb color for multiples of p
313
remaining -- rgb color for integers that have not been sieved out yet
314
and are not multiples of p.
315
"""
316
if p == 1:
317
c = [remaining]*n
318
else:
319
exclude = prod(prime_range(p)) # exclude multiples of primes < p
320
c = []
321
for k in range(1,n+1):
322
if k <= p and is_prime(k):
323
c.append(prime)
324
elif k == 1 or (gcd(k,exclude) != 1 and not is_prime(k)):
325
c.append(gone)
326
elif k%p == 0:
327
c.append(multiple)
328
else:
329
c.append(remaining)
330
# end for
331
# end if
332
return number_grid(c,{'rgbcolor':(0.2,0.2,0.2)},{'fontsize':fontsize, 'rgbcolor':(0,0,0)},offset=offset)
333
334
335
##############################################################
336
# Similar rates of growth
337
##############################################################
338
339
def fig_simrates(dir,ext):
340
# similar rates
341
G = similar_rates()
342
G.save(dir + "/similar_rates.%s"%ext,figsize=[8,3], typeset='latex')
343
344
def similar_rates():
345
"""
346
Draw figure fig:simrates illustrating similar rates.
347
348
EXAMPLES::
349
350
sage: similar_rates()
351
"""
352
var('X')
353
A = 2*X^2 + 3*X - 5
354
B = 3*X^2 - 2*X + 1
355
G = plot(A/B, (1,100))
356
G += text("$A(X)/B(X)$", (70,.58), rgbcolor='black', fontsize=14)
357
H = plot(A, (X,1,100), rgbcolor='red') + plot(B, (X,1,100))
358
H += text("$A(X)$", (85,8000), rgbcolor='black',fontsize=14)
359
H += text("$B(X)$", (60,18000), rgbcolor='black',fontsize=14)
360
a = graphics_array([[H,G]])
361
return a
362
363
364
365
366
##############################################################
367
# Proportion of Primes to to X
368
##############################################################
369
370
def fig_log(dir, ext):
371
g = plot(log, 1/3, 100, thickness=2)
372
g.save(dir + '/log.%s'%ext, figsize=[8,3], gridlines=True, fontsize=18, typeset='latex')
373
374
def fig_proportion_primes(dir,ext):
375
for bound in [100,1000,10000]:
376
g = proportion_of_primes(bound)
377
g.save(dir + '/proportion_primes_%s.%s'%(bound,ext), typeset='latex')
378
379
def plot_step_function(v, vertical_lines=True, **args):
380
r"""
381
Return the line that gives the plot of the step function f defined
382
by the list v of pairs (a,b). Here if (a,b) is in v, then f(a) = b.
383
384
INPUT:
385
386
- `v` -- list of pairs (a,b)
387
388
EXAMPLES::
389
390
sage: plot_step_function([(i,sin(i)) for i in range(5,20)])
391
"""
392
# make sorted copy of v (don't change in place, since that would be rude).
393
v = list(sorted(v))
394
if len(v) <= 1:
395
return line([]) # empty line
396
if vertical_lines:
397
w = []
398
for i in range(len(v)):
399
w.append(v[i])
400
if i+1 < len(v):
401
w.append((v[i+1][0],v[i][1]))
402
return line(w, **args)
403
else:
404
return sum(line([v[i],(v[i+1][0],v[i][1])], **args) for i in range(len(v)-1))
405
406
407
def proportion_of_primes(bound, **args):
408
"""
409
Return a graph of the step function that assigns to X the
410
proportion of numbers between 1 and X of primes.
411
412
INPUTS:
413
414
- `bound` -- positive integer
415
416
- additional arguments are passed to the line function.
417
418
EXAMPLES::
419
420
sage: proportion_of_primes(100)
421
"""
422
v = []
423
k = 0
424
for n in range(1,bound+1):
425
if is_prime(n):
426
k += 1
427
v.append((n,k/n))
428
return plot_step_function(v, **args)
429
430
431
##############################################################
432
# Prime Gaps
433
##############################################################
434
435
@cached_function
436
def prime_gaps(maxp):
437
"""
438
Return the sequence of prime gaps obtained using primes up to maxp.
439
440
EXAMPLES::
441
442
sage: prime_gaps(100)
443
[1, 2, 2, 4, 2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8]
444
"""
445
P = prime_range(maxp+1)
446
return [P[i+1] - P[i] for i in range(len(P)-1)]
447
448
@cached_function
449
def prime_gap_distribution(maxp):
450
"""
451
Return list v such that v[i] is how many times i is a prime gap
452
among the primes up to maxp.
453
454
EXAMPLES::
455
456
sage: prime_gap_distribution(100)
457
[0, 1, 8, 0, 7, 0, 7, 0, 1]
458
sage: prime_gap_distribution(1000)
459
[0, 1, 35, 0, 40, 0, 44, 0, 15, 0, 16, 0, 7, 0, 7, 0, 0, 0, 1, 0, 1]
460
"""
461
h = prime_gaps(maxp)
462
v = [0]*(max(h)+1)
463
for gap in h: v[gap] += 1
464
return v
465
466
def fig_primegapdist(dir,ext):
467
v = prime_gap_distribution(10^7)[:50]
468
b = bar_chart(v)
469
b.save(dir+"/primegapdist.%s"%ext, fontsize=18, figsize=[9,3], ticks=[[2,4,6,8]+[10,20,..,40]+[48],[20000,60000,90000]], tick_formatter=['latex','latex'], typeset='latex')
470
471
"""
472
# The table in the book...
473
474
def f(B):
475
v = prime_gap_distribution(10^B)
476
z = [v[i] if i<len(v) else 0 for i in [2,4,6,8,100,252]]
477
print "$10^{%s}$ & "%B + " & ".join([str(a) for a in z]) + r"\\\hline"
478
479
for B in [1..8]:
480
f(B)
481
"""
482
483
484
def prime_gap_plots(maxp, gap_sizes):
485
"""
486
Return a list of graphs of the functions Gap_k(X) for 0<=X<=maxp,
487
for each k in gap_sizes. The graphs are lists of pairs (X,Y) of
488
integers.
489
490
INPUT:
491
- maxp -- positive integer
492
- gap_sizes -- list of integers
493
"""
494
P = prime_range(maxp+1)
495
v = [[(0,0)] for i in gap_sizes]
496
k = dict([(g,i) for i, g in enumerate(gap_sizes)])
497
for i in range(len(P)-1):
498
g = P[i+1] - P[i]
499
if g in k:
500
w = v[k[g]]
501
w.append( (P[i+1],w[-1][1]) )
502
w.append( (P[i+1],w[-1][1]+1) )
503
return v
504
505
506
def fig_primegap_race(dir, ext):
507
"""
508
Draw plots showing the race for gaps of size 2, 4, 6, and 8.
509
"""
510
X = 7000
511
gap_sizes = [2,4,6,8]
512
#X = 100000
513
#gap_sizes = [i for i in range(2,50) if i%2==0]
514
v = prime_gap_plots(X, gap_sizes)
515
516
P = sum(line(x) for x in v)
517
P += sum( text( "Gap %s"%gap_sizes[i], (v[i][-1][0]*1.04, v[i][-1][1]), color='black', fontsize=8)
518
for i in range(len(v)))
519
520
P.save(dir + '/primegap_race.%s'%ext, figsize=[9,3], gridlines=True, typeset='latex')
521
return P
522
523
524
525
##############################################################
526
# Multiplicatively even and odd table
527
##############################################################
528
def mult_even_odd_count(bound):
529
"""
530
Return list v of length bound such that v[n] is a pair (a,b) where
531
a is the number of multiplicatively even positive numbers <= n and b is the
532
number of multiplicatively odd positive numbers <= n.
533
534
INPUT:
535
536
- ``bound`` -- a positive integer
537
538
EXAMPLES::
539
540
We make the table in the paper::
541
542
sage: mult_even_odd_count(17)
543
[(0, 0), (1, 0), (1, 1), (1, 2), (2, 2), (2, 3), (3, 3), (3, 4), (3, 5), (4, 5), (5, 5), (5, 6), (5, 7), (5, 8), (6, 8), (7, 8), (8, 8)]
544
"""
545
par = mult_parities_python(bound)
546
a = 0; b = 0
547
v = [(a,b)]
548
for n in range(1,bound):
549
if par[n] == 0: # mult even
550
a += 1
551
else:
552
b += 1
553
v.append((a,b))
554
return v
555
556
def fig_multpar(dir,ext):
557
for n in [10^k for k in reversed([1,2,3,4,5,6])] + [17]:
558
file = dir + '/liouville-%s.%s'%(n,ext)
559
if n >= 1000:
560
time_series = True
561
else:
562
time_series = False
563
g = race_mult_parity(n, time_series=time_series)
564
g.save(file, frame=True, fontsize=20, typeset='latex')
565
566
def race_mult_parity(bound, time_series=False, **kwds):
567
"""
568
Return a plot that shows race between multiplicatively even and
569
odd numbers. More precisely it shows the function f(X) that
570
equals number of even numbers >= 2 and <=X minus the number of odd
571
numbers >= 2 and <= X.
572
573
EXAMPLES::
574
575
sage: race_mult_parity(10^5,time_series=True)
576
sage: race_mult_parity(10^5)
577
"""
578
par = mult_parities_python(bound)[2:]
579
if not time_series:
580
v = [(2,0)]
581
for n in range(bound-2):
582
if par[n] == 0:
583
b = v[-1][1]+1
584
else:
585
b = v[-1][1]-1
586
v.append((v[-1][0], b))
587
v.append((v[-1][0]+1, b))
588
return line(v, **kwds)
589
else:
590
v = [0,0,0]
591
for n in range(bound-2):
592
if par[n] == 0:
593
v.append(v[-1]+1)
594
else:
595
v.append(v[-1]-1)
596
return finance.TimeSeries(v).plot()
597
598
def mult_parities_python(bound, verbose=False):
599
"""
600
Return list v of length bound such that v[n] is 0 if n is
601
multiplicative even, and v[n] is 1 if n is multiplicatively odd.
602
Also v[0] is None.
603
604
This goes up to bound=`10^7` in about 30 seconds.
605
"""
606
v = [None]*bound
607
v[0] = None
608
v[1] = int(0)
609
P = [int(p) for p in prime_range(bound)]
610
for p in P:
611
v[p] = int(1)
612
last = P
613
last_parity = int(1)
614
loops = floor(log(bound,2))+1
615
bound = int(bound)
616
for k in range(loops):
617
cur = []
618
cur_parity = (last_parity+int(1))%int(2)
619
if verbose:
620
print "loop %s (of %s); last = %s"%(k,loops, len(last))
621
for n in last:
622
for p in P:
623
m = n * p
624
if m >= bound:
625
break
626
if v[m] is None:
627
v[m] = cur_parity
628
cur.append(m)
629
last_parity = cur_parity
630
last = cur
631
return v
632
633
##############################################################
634
# LogX over X in "Probabilistic first guess" chapter
635
##############################################################
636
def fig_logXoverX(dir, ext):
637
file = dir + '/logXoverX.%s'%ext
638
x = var('x')
639
xmax = 250
640
G = plot(x/(log(x)-1), 4, xmax, color='blue')
641
G += prime_pi.plot(4, xmax, color='red')
642
G.save(file, figsize=[7,3], typeset='latex')
643
644
##############################################################
645
# Prime pi plots
646
##############################################################
647
648
def fig_prime_pi_aspect1(dir,ext):
649
for n in [25, 100]:
650
p = plot(lambda x: prime_pi(floor(x)), 1,n,
651
plot_points=10000,rgbcolor='red',
652
fillcolor=(.9,.9,.9),fill=True)
653
file = dir + '/prime_pi_%s_aspect1.%s'%(n, ext)
654
p.save(file, aspect_ratio=1, fontsize=16, typeset='latex')
655
656
def fig_prime_pi(dir,ext):
657
for n in [1000, 10000, 100000]:
658
p = plot(lambda x: prime_pi(floor(x)), 1,n,
659
plot_points=10000,rgbcolor='red',
660
fillcolor=(.9,.9,.9),fill=True)
661
file = dir + '/prime_pi_%s.%s'%(n, ext)
662
if n <100000:
663
p.save(file, fontsize=16, typeset='latex')
664
else:
665
p.save(file, fontsize=16, ticks=[[n/2, n], None], tick_formatter=['latex', None], typeset='latex')
666
667
def fig_prime_pi_nofill(dir,ext):
668
for n in [25,38,100,1000,10000]:
669
g = plot_prime_pi(n, rgbcolor='red', thickness=2, fontsize=20)
670
g.save(dir + '/PN_%s.%s'%(n,ext), typeset='latex')
671
n = 100000
672
g = plot_prime_pi(n, rgbcolor='red', thickness=1, fontsize=16, ticks=[[20000,60000,100000],[2000,5000,8000]], tick_formatter=['latex','latex'])
673
g.save(dir + '/PN_%s.%s'%(n,ext), typeset='latex')
674
675
def plot_prime_pi(n = 25, **args):
676
v = [(0,0)]
677
k = 0
678
for p in prime_range(n+1):
679
k += 1
680
v.append((p,k))
681
v.append((n,k))
682
return plot_step_function(v, **args)
683
684
def fig_pi_Li(dir, ext):
685
g = plot_prime_pi(n = 250, rgbcolor='red', thickness=1, fontsize=16)
686
h = plot(Li, 1, 250)
687
(g+h).save(dir + 'pi_Li.%s'%ext, typeset='latex')
688
689
##############################################################
690
# Sieving
691
##############################################################
692
693
def fig_sieves(dir,ext):
694
plot_three_sieves(100, shade=False).save(dir + '/sieve_2_100.%s'%ext, figsize=[9,3], fontsize=18, typeset='latex')
695
plot_all_sieves(1000, shade=True).save(dir + '/sieve1000.%s'%ext,figsize=[9,3], fontsize=18, typeset='latex')
696
697
m=100
698
for z in [3,7]:
699
save(plot_multiple_sieves(m,k=[z]) ,dir+'/sieves%s_100.%s'%(z,ext), xmax=m, figsize=[9,3], fontsize=18)
700
701
def plot_sieve(n, x, poly={}, lin={}, label=True, shade=True):
702
"""
703
Return a plot of the number of integer up to x that are coprime to n.
704
These are the integers that are sieved out using the primes <= n.
705
706
In n is 0 draw a graph of all primes.
707
"""
708
v = range(x+1) # integers 0, 1, ..., x
709
if n == 0:
710
v = prime_range(x)
711
else:
712
for p in prime_divisors(n):
713
v = [k for k in v if k%p != 0 or k==p]
714
# eliminate non-prime multiples of p
715
v = set(v)
716
j = 0
717
w = [(0,j)]
718
for i in range(1,x+1):
719
w.append((i,j))
720
if i in v:
721
j += 1
722
w.append((i,j))
723
w.append((i,0))
724
w.append((0,0))
725
if n == 0:
726
t = "Primes"
727
pos = x,.7*j
728
elif n == 1:
729
t = "All Numbers"
730
pos = x, 1.03*j
731
else:
732
P = prime_divisors(n)
733
if len(P) == 1:
734
t = "Sieve by %s"%P[0]
735
else:
736
t = "Sieve by %s"%(', '.join([str(_) for _ in P]))
737
pos = x,1.05*j
738
F = line(w[:-2], **lin)
739
if shade:
740
F += polygon(w, **poly)
741
if label:
742
F += text(t, pos, horizontal_alignment="right", rgbcolor='black', fontsize=18)
743
F.set_aspect_ratio('automatic')
744
return F
745
746
def plot_three_sieves(m, shade=True):
747
s1 = plot_sieve(1, m, poly={'rgbcolor':(.85,.9,.7)},
748
lin={'rgbcolor':(0,0,0), 'thickness':1}, shade=shade)
749
s2 = plot_sieve(2, m, poly={'rgbcolor':(.75,.8,.6)},
750
lin={'rgbcolor':(0,0,0),'thickness':1}, shade=shade)
751
s3 = plot_sieve(0, m, poly={'rgbcolor':(1,.7,.5)},
752
lin={'rgbcolor':(1,0,0), 'thickness':1}, shade=shade)
753
return s1+s2+s3
754
755
def plot_multiple_sieves(m=100, k = [2,3,5], shade=True):
756
g = Graphics()
757
n = len(k)
758
for i in range(n):
759
c = (1-float(i+1)/n)*0.666
760
if k[i] == 0:
761
z = 0
762
else:
763
z = prod(prime_range(k[i]+1))
764
r = float(i)/n
765
clr = (.85 + 0.15*r,0.9 -0.2*r, 0.9 -0.4*r)
766
if z == 0:
767
clrlin=(1,0,0)
768
else:
769
clrlin=(0,0,0)
770
s = plot_sieve(z, m,
771
poly={'rgbcolor':clr},
772
lin={'rgbcolor':clrlin, 'thickness':1},
773
label=(i==0 or i==n-1), shade=shade)
774
g += s
775
return g
776
777
def plot_all_sieves(x, shade=True):
778
P = [1] + prime_range(int(sqrt(x))+1) + [0]
779
G = plot_multiple_sieves(x, P, shade=shade)
780
return G
781
782
783
##############################################################
784
# Area under plot of 1/log(x)
785
##############################################################
786
787
#auto
788
def area_under_inverse_log(m, **args):
789
r"""
790
This function returns a graphical illustration of `Li(x)` for `x
791
\leq m` viewed as the integral of `1/\log(t)` from 2 to `t`. We
792
also plot primes on the `x`-axis and display the area as text.
793
794
EXAMPLES::
795
796
797
"""
798
f = plot(lambda x: 1/(log(x)),2,m)
799
P = polygon([(2,0)]+list(f[0])+[(m,0)], hue=0.1,alpha=0.4)
800
if False:
801
T = sum([text(str(p),(p,-0.08),vertical_alignment="top",\
802
horizontal_alignment="center", fontsize=6, rgbcolor=(.6,0,0)) \
803
for p in prime_range(m+1)])
804
else:
805
T = Graphics()
806
pr = sum([point((p,0), rgbcolor=(.6,0,0), pointsize=100/log(p)) for p in prime_range(m+1)])
807
L = 'quad_qag(1/log(x), x, 2,%s, 0)'%m
808
fs = 36
809
area = text('Area ~ %f'%(float(maxima(L)[0])),(.5*m,.75),fontsize=fs,rgbcolor='black')
810
primes = text('%s Primes'%len(prime_range(m+1)), (.5*m,-0.3),fontsize=fs,rgbcolor='black')
811
fun = text('1/log(x)',(m/8,1.4),fontsize=fs,rgbcolor='black', horizontal_alignment='left')
812
G = pr + f+P+area+T+primes +fun
813
G.xmax(m+1)
814
G.set_aspect_ratio('automatic')
815
return G
816
817
def fig_inverse_of_log(dir,ext):
818
for m in [30, 100, 1000]:
819
area_under_inverse_log(m).save(dir+'/area_under_log_graph_%s.%s'%(m,ext), figsize=[7,7], typeset='latex')
820
821
822
##############################################################
823
# Comparing Li, pi, and x/log(x)
824
##############################################################
825
def fig_li_pi_loginv(dir,ext):
826
plot_li_pi_loginv(xmax=200).save(dir+'/three_plots.%s'%ext,figsize=[8,3], typeset='latex')
827
828
def plot_li_pi_loginv(xmax=200):
829
var('x')
830
P = plot(x/log(x), (2, xmax))
831
P+= plot(Li, (2, xmax), rgbcolor='black')
832
P+= plot(prime_pi, 2, xmax, rgbcolor='red')
833
return P
834
835
836
##############################################################
837
# Perspective
838
##############################################################
839
840
def fig_primes_line(dir,ext):
841
xmin=1; xmax=38; pointsize=90
842
g = points([(p,0) for p in prime_range(xmax+1)], pointsize=pointsize, rgbcolor='red')
843
g += line([(xmin,0), (xmax,0)], rgbcolor='black')
844
eps = 1/2
845
for n in [xmin..xmax]:
846
g += line([(n,eps), (n,-eps)], rgbcolor='black', thickness=0.5)
847
g += text("$%s$"%n, (n,-6), rgbcolor='black')
848
g.save(dir + '/primes_line.%s'%ext, axes=False,figsize=[9,.7], ymin=-10, ymax=2, typeset='latex')
849
850
##############################################################
851
# Plots of Psi function
852
##############################################################
853
854
def psi_data(xmax):
855
from math import log, pi
856
v = [(0,0), (1,0), (1,log(2*pi))]
857
y = v[-1][1]
858
for pn in prime_powers(2,xmax+1):
859
y += log(factor(pn)[0][0])
860
v.append( (pn,y) )
861
v.append((xmax,y))
862
return v
863
864
def plot_psi(xmax, **kwds):
865
v = psi_data(xmax)
866
return plot_step_function(v, **kwds)
867
868
def fig_psi(dir,ext):
869
for m in [9,38,100,200]:
870
g = plot_psi(m, thickness=2)
871
g.save(dir+'/psi_%s.%s'%(m,ext), aspect_ratio=1, gridlines=True,
872
fontsize=16, figsize=4, typeset='latex')
873
g = plot(lambda x:x,1,1000,rgbcolor='red', thickness=.7) + plot_psi(1000,alpha=0.8, thickness=.7)
874
g.save(dir+'/psi_diag_%s.%s'%(1000,ext),aspect_ratio=1, gridlines=True, fontsize=12, figsize=4, typeset='latex')
875
876
def plot_Psi(xmax, **kwds):
877
v = psi_data(xmax)
878
v = [(log(a),b) for a, b in v if a]
879
return plot_step_function(v, **kwds)
880
881
def fig_Psi(dir, ext):
882
m = 38
883
g = plot_Psi(m, thickness=2)
884
g.save(dir+'/bigPsi_%s.%s'%(m,ext), gridlines=True,
885
fontsize=20, figsize=[6.1,6.1], typeset='latex')
886
887
def fig_Psiprime(dir, ext):
888
g = line([(0,0),(0,100)], rgbcolor='black')
889
xmax = 20
890
for n in [1..xmax]:
891
if is_prime_power(n):
892
if n == 1:
893
h = log(2*pi)
894
else:
895
h = log(factor(n)[0][0])
896
h = float(h)*50
897
g += arrow((log(n),-1),(log(n),h), width=2)
898
g += line([(log(n),-2), (log(n),150)], rgbcolor='black')
899
if n <= 3 or n in [5, 8, 13, 19]:
900
g += text("log(%s)"%n, (log(n),-8), rgbcolor='black', fontsize=12)
901
g += line([(log(n),-2), (log(n),0)], rgbcolor='black')
902
g += line([(-1/2,0), (xmax+1,0)], thickness=2)
903
g.save(dir+'/bigPsi_prime.%s'%ext,
904
xmin=-1/2, xmax=log(xmax),
905
axes=False, gridlines=True, figsize=[8,3], typeset='latex')
906
907
def fig_Phi(dir=0, ext=0):
908
g = line([(0,0),(0,100)], rgbcolor='black')
909
xmax = 20
910
ymax = 50
911
for n in [1..xmax]:
912
if is_prime_power(n):
913
if n == 1:
914
h = log(2*pi)
915
else:
916
h = log(factor(n)[0][0])
917
h /= sqrt(n)
918
h = float(h)*50
919
g += arrow((log(n),0),(log(n),h), width=1)
920
g += arrow((-log(n),0),(-log(n),h), width=1)
921
g += line([(log(n),0),(log(n),100)], color='black', thickness=.3)
922
g += line([(-log(n),0),(-log(n),100)], color='black', thickness=.3)
923
if n in [2, 5, 16]:
924
g += text("log(%s)"%n, (log(n),-5), rgbcolor='black', fontsize=12)
925
g += line([(log(n),-2), (log(n),0)], rgbcolor='black')
926
g += text("log(%s)"%n, (-log(n),-5), rgbcolor='black', fontsize=12)
927
g += line([(-log(n),-2), (-log(n),0)], rgbcolor='black')
928
g += line([(-log(xmax)-1,0), (log(xmax)+1,0)], thickness=2)
929
930
g.save(dir+'/bigPhi.%s'%ext,
931
xmin=-log(xmax), xmax=log(xmax), ymax=ymax,
932
axes=False, figsize=[8,3], typeset='latex')
933
934
935
##############################################################
936
# Sin, etc. waves
937
##############################################################
938
939
def fig_waves(dir,ext):
940
g = plot(sin, -2.1*pi, 2.1*pi, thickness=2)
941
g.save(dir+'/sin.%s'%ext, fontsize=20, typeset='latex')
942
943
x = var('x')
944
c = 5
945
# See for why this is right http://www.phy.mtu.edu/~suits/notefreqs.html
946
g = plot(sin(x), 0, c*pi) + plot(sin(329.0/261*x), 0, c*pi, color='red')
947
g.save(dir+'/sin-twofreq.%s'%ext, fontsize=20, typeset='latex')
948
949
g = plot(sin(x) + sin(329.0/261*x), 0, c*pi)
950
g.save(dir+'/sin-twofreq-sum.%s'%ext, fontsize=20, typeset='latex')
951
952
c=5
953
g = plot(sin(x), -2, c*pi) + plot(sin(x + 1.5), -2, c*pi, color='red')
954
g += text("Phase", (-2.5,.5), fontsize=14, rgbcolor='black')
955
g += arrow((-2.5,.4), (-1.5,0), width=1, rgbcolor='black')
956
g += arrow((-2,.4), (0,0), width=1, rgbcolor='black')
957
g.save(dir+'/sin-twofreq-phase.%s'%ext, fontsize=20, typeset='latex')
958
959
g = plot(sin(x) + sin(329.0/261*x + 0.4), 0, c*pi)
960
g.save(dir+'/sin-twofreq-phase-sum.%s'%ext, fontsize=20, typeset='latex')
961
962
f(x) = sin(x) + sin(329.0/261*x + 0.4)
963
g = points([(i,f(i)) for i in [0,0.1,..,5*pi]])
964
g.save(dir+'/sin-twofreq-phase-sum-points.%s'%ext, fontsize=20, typeset='latex')
965
966
g = points([(i,f(i)) for i in [0,0.1,..,5*pi]])
967
g += plot(f, (0, 5*pi), rgbcolor='black')
968
g.save(dir+'/sin-twofreq-phase-sum-fill.%s'%ext, fontsize=20, typeset='latex')
969
970
f(x) = 0.7*sin(x) + sin(329.0/261*x + 0.4)
971
g = plot(f, (0, 5*pi))
972
g.save(dir+'/sound-ce-general_sum.%s'%ext, fontsize=20, typeset='latex')
973
974
B = bar_chart([0,0,0,0.7, 0, 1])
975
B += text("C", (3.2,-0.05), rgbcolor='black', fontsize=18)
976
B += text("D", (4.2,-0.05), rgbcolor='black', fontsize=18)
977
B += text("E", (5.2,-0.05), rgbcolor='black', fontsize=18)
978
B.save(dir+'/sound-ce-general_sum-blips.%s'%ext, axes=False, xmin=0, fontsize=20, typeset='latex')
979
980
f(x) = 0.7*sin(x) + sin(329.0/261*x + 0.4) + 0.5*sin(300.0/261*x + 0.7) + 0.3*sin(1.5*x + 0.2) + 1.1*sin(4*x+0.1)
981
g = plot(f, (0, 5*pi))
982
g.save(dir + '/complicated-wave.%s'%ext, fontsize=20, typeset='latex')
983
984
##############################################################
985
# Sawtooth
986
##############################################################
987
988
def fig_sawtooth(dir,ext):
989
g = plot_sawtooth(3)
990
g.save(dir+'/sawtooth.%s'%ext, figsize=[8,3], fontsize=20, typeset='latex')
991
992
g = plot_sawtooth_spectrum(18)
993
g.save(dir+'/sawtooth-spectrum.%s'%ext, figsize=[8,3], fontsize=20, typeset='latex')
994
995
def plot_sawtooth(xmax):
996
v = []
997
for x in [0..xmax]:
998
v += [(x,0), (x+1,1), (x+1,0)]
999
return line(v)
1000
1001
def plot_sawtooth_spectrum(xmax):
1002
# the spectrum is a spike of height 1/k at k
1003
return sum([line([(k,0),(k,1/k)],thickness=3) for k in [1..xmax]])
1004
1005
1006
##############################################################
1007
# Pure tone
1008
##############################################################
1009
def fig_pure_tone(dir,ext):
1010
t = var('t')
1011
g = plot(2 * cos(1+t/2), -15, 15, thickness=2)
1012
g.save(dir+'/pure_tone.%s'%ext, figsize=[8,3], fontsize=20, typeset='latex')
1013
1014
1015
1016
##############################################################
1017
# Mixed tone
1018
##############################################################
1019
def fig_mixed_tone(dir,ext):
1020
t = var('t')
1021
f = 5 *cos(-t - 2) + 2 * cos(t/2 + 1) + 3 *cos(2 *t + 4)
1022
g = plot(f, -15, 15, thickness=2)
1023
g.save(dir+'/mixed_tone3.%s'%ext, figsize=[8,3], fontsize=20, typeset='latex')
1024
1025
##############################################################
1026
# Fourier Transforms: second visit
1027
##############################################################
1028
def fig_even_function(dir, ext):
1029
x = var('x')
1030
f = cos(x) + sin(x^2) + sqrt(x)
1031
def g(z):
1032
return f(x=abs(z))
1033
h = plot(g,-4,4)
1034
h.save(dir + '/even_function.%s'%ext, figsize=[9,3], fontsize=18, typeset='latex')
1035
1036
def fig_even_pi(dir, ext):
1037
g1 = prime_pi.plot(0,50, rgbcolor='red')
1038
g2 = prime_pi.plot(0,50, rgbcolor='red')
1039
g2[0].xdata = [-a for a in g2[0].xdata]
1040
g = g1 + g2
1041
g.save(dir + '/even_pi.%s'%ext, figsize=[10,3],xmin=-49,xmax=49, fontsize=18, typeset='latex')
1042
1043
def fig_oo_integral(dir, ext):
1044
t = var('t')
1045
f(t) = 1/(t^2+1)*cos(t)
1046
a = f.find_root(1,2.5)
1047
b = f.find_root(4,6)
1048
c = f.find_root(7,8)
1049
g = plot(f,(t,-13,13), fill='axis', fillcolor='yellow', fillalpha=1, thickness=2)
1050
g += plot(f,(t,-a,a), fill='axis', fillcolor='grey', fillalpha=1, thickness=2)
1051
g += plot(f,(t,-c,-b), fill='axis', fillcolor='grey', fillalpha=1, thickness=2)
1052
g += plot(f,(t,b,c), fill='axis', fillcolor='grey', fillalpha=1, thickness=2)
1053
g += text(r"$\int_{-\infty}^{\,\infty} f(x) dx$", (-7,0.5), rgbcolor='black', fontsize=30)
1054
#g.show(figsize=[9,3], xmin=-10,xmax=10)
1055
g.save(dir+'/oo_integral.%s'%ext, figsize=[9,3], xmin=-10,xmax=10, typeset='latex')
1056
1057
def fig_fourier_machine(dir, ext):
1058
g = text("$f(t)$", (-1/2,1/2), fontsize=30, rgbcolor='black')
1059
g += text(r"$\hat{f}(\theta)$", (3/2,1/2), fontsize=30, rgbcolor='black')
1060
g += line([(0,0),(0,1),(1,1),(1,0),(0,0)],rgbcolor='black',thickness=3)
1061
g += arrow((-1/2+1/9,1/2), (-1/16,1/2), rgbcolor='black')
1062
g += arrow((1+1/16,1/2), (1+1/2-1/9,1/2), rgbcolor='black')
1063
t=var('t')
1064
g += plot((1/2)*t*cos(14*t)+1/2,(t,0,1), fill='axis', thickness=0.8)
1065
g.save(dir+'/fourier_machine.%s'%ext, axes=False, axes_pad=0.1, typeset='latex')
1066
1067
1068
##############################################################
1069
# Distribution section
1070
##############################################################
1071
1072
def fig_simple_staircase(dir,ext):
1073
v = [(-1,0), (0,1), (1,3), (2,3)]
1074
g = plot_step_function(v, thickness=3, vertical_lines=True)
1075
g.save(dir+'/simple_staircase.%s'%ext, fontsize=20, typeset='latex')
1076
1077
1078
##############################################################
1079
# Plots of Fourier transform Phihat_even.
1080
##############################################################
1081
1082
def fig_mini_phihat_even(dir,ext):
1083
G = plot_symbolic_phihat(5, 1, 100, 10000, zeros=False)
1084
G.save(dir + "/mini_phihat_even.%s"%ext, figsize=[9,3], ymin=0, fontsize=18, typeset='latex')
1085
1086
def fig_phihat_even(dir,ext):
1087
for bound in [5, 20, 50, 500]:
1088
G = plot_symbolic_phihat(bound, 2, 100,
1089
plot_points=10^5)
1090
G.save(dir+'/phihat_even-%s.%s'%(bound,ext), figsize=[9,3], ymin=0, fontsize=18, typeset='latex')
1091
1092
def fig_phihat_even_all(dir, ext):
1093
p = [plot_symbolic_phihat(n, 2,100) for n in [5,20,50,500]]
1094
[a.ymin(0) for a in p]
1095
g = graphics_array([[a] for a in p],4,1)
1096
g.save(dir+'/phihat_even_all.%s'%ext, typeset='latex')
1097
1098
def symbolic_phihat(bound):
1099
t = var('t')
1100
f = SR(0)
1101
for pn in prime_powers(bound+1):
1102
if pn == 1: continue
1103
p, e = factor(pn)[0]
1104
f += - log(p)/sqrt(pn) * cos(t*log(pn))
1105
return f
1106
1107
def plot_symbolic_phihat(bound, xmin, xmax, plot_points=1000, zeros=True):
1108
f = symbolic_phihat(bound)
1109
P = plot(f, (t,xmin, xmax), plot_points=plot_points)
1110
if not zeros:
1111
return P
1112
ym = P.ymax()
1113
Z = []
1114
for y in zeta_zeros():
1115
if y > xmax: break
1116
Z.append(y)
1117
zeros = sum([arrow((x,ym),(x,0),rgbcolor='red',width=0.5,arrowsize=2)
1118
for i, x in enumerate(Z)])
1119
return P + zeros
1120
1121
##############################################################
1122
# Calculus pictures
1123
##############################################################
1124
1125
def fig_aplusone(dir,ext):
1126
a = var('a')
1127
g = plot(a+1, -5,8, thickness=3)
1128
g.save(dir + '/graph_aplusone.%s'%ext, gridlines=True, frame=True, fontsize=20, typeset='latex')
1129
1130
def fig_calculus(dir,ext):
1131
x = var('x')
1132
t = 8; f = log(x); fprime = f.diff()
1133
fontsize = 16
1134
g = plot(f, (0.5,t), thickness=2)
1135
g += plot(x*fprime(x=4)+(f(x=4)-4*fprime(x=4)), (-.5,t), rgbcolor='green', thickness=2, zorder=10)
1136
g += point((4,f(x=4)), pointsize=50, rgbcolor='black', zorder=20)
1137
g += plot(fprime, (0.5,t), rgbcolor='red', thickness=2, zorder=15)
1138
g += text("What is the slope of the tangent line?", (3.5,2.2),
1139
fontsize=fontsize, rgbcolor='black')
1140
g += text("Here it is!",(5,.8), fontsize=fontsize, rgbcolor='black')
1141
g += arrow((4.7,.64), (4.05, fprime(x=4.05)+.08), rgbcolor='black')
1142
g += point((4,fprime(x=4)),rgbcolor='black', pointsize=50, zorder=20)
1143
g += text("How to compute the slope? This is Calculus.", (4.3, -0.3),
1144
fontsize=fontsize, rgbcolor='black')
1145
g.save(dir + '/graph_slope_deriv.%s'%ext, gridlines=True, frame=True, fontsize=16, typeset='latex')
1146
1147
def fig_jump(dir,ext):
1148
# straight jump
1149
v = line( [(0,1), (3,1), (3,2), (6,2)], thickness=2)
1150
v.ymin(0)
1151
v.save(dir + '/jump.%s'%ext, fontsize=20, typeset='latex')
1152
1153
# smooth approximation to a jump
1154
e = .7
1155
v = line( [(0,1), (3-e,1)], thickness=2) + line([(3+e,2), (6,2)], thickness=2)
1156
v.ymin(0)
1157
S = spline( [(3-e,1), (3-e+e/20, 1), (3,1.5), (3+e-e/20, 2), (3+e,2)] )
1158
v += plot(S, (3-e, 3+e), thickness=2)
1159
v.save(dir + '/jump-smooth.%s'%ext, fontsize=20, typeset='latex')
1160
1161
# derivatives of smooth jumps
1162
for e in ['7', '2', '05', '01']:
1163
g = smoothderiv(float('.'+e))
1164
g.save(dir + '/jump-smooth-deriv-%s.%s'%(e,ext), fontsize=20, typeset='latex')
1165
1166
1167
def smoothderiv(e):
1168
def deriv(f, delta):
1169
# return very approximate numerical derivative of callable f, using
1170
# a given choice of delta
1171
def g(x):
1172
return (f(x + delta) - f(x))/delta
1173
return g
1174
1175
v = line( [(0,1), (3-e,1)], thickness=2) + line([(3+e,2), (6,2)], thickness=2)
1176
v.ymin(0)
1177
S = spline( [(3-e,1), (3-e+e/20, 1), (3,1.5), (3+e-e/20, 2), (3+e,2)] )
1178
v += plot(S, (3-e, 3+e), thickness=2)
1179
D = (line( [(0,0), (3-e,0)], rgbcolor='red', thickness=2) +
1180
line([(3+e,0), (6,0)], rgbcolor='red', thickness=2))
1181
D += plot( deriv(S, e/30), (3-e, 3+e), rgbcolor='red', thickness=2)
1182
v += D
1183
return v
1184
1185
1186
def fig_dirac(dir,ext):
1187
g = line([(0,0),(0,100)], rgbcolor='black')
1188
g += arrow((0,-1),(0,50), width=3)
1189
g += line([(-1.2,0), (1.25,0)], thickness=3)
1190
g.save(dir+'/dirac_delta.%s'%ext,
1191
frame=False, xmin=-1, xmax=1, ymax=50, axes=False, gridlines=True, fontsize=20, typeset='latex')
1192
1193
def fig_two_delta(dir,ext):
1194
g = line([(0,0),(0,100)], rgbcolor='black')
1195
g += arrow((-1/2,-1),(-1/2,50),width=3)
1196
g += arrow((1/2,-1),(1/2,50),width=3)
1197
g += line([(-1.2,0), (1.25,0)], thickness=3)
1198
g += text("$-x$", (-1/2-1/20,-4), rgbcolor='black', fontsize=35, horizontal_alignment='center')
1199
g += text("$x$", (1/2,-4), rgbcolor='black', fontsize=35, horizontal_alignment='center')
1200
g.save(dir+'/two_delta.%s'%ext,
1201
frame=False, xmin=-1, xmax=1, ymax=50, axes=False, gridlines=True, fontsize=20, typeset='latex')
1202
1203
##############################################################
1204
# Cosine sums
1205
##############################################################
1206
1207
def fig_phi(dir,ext):
1208
fontsize = 18
1209
g = phi_approx_plot(2,30,1000, fontsize=fontsize)
1210
g.save(dir+'/phi_cos_sum_2_30_1000.%s'%ext, axes=False, ymin=-5,ymax=160, fontsize=fontsize, typeset='latex')
1211
1212
g = phi_interval_plot(26, 34, fontsize=fontsize)
1213
g.save(dir+'/phi_cos_sum_26_34_1000.%s'%ext, axes=False, fontsize=fontsize, typeset='latex')
1214
1215
g = phi_interval_plot(1010,1026,15000,drop=60, fontsize=fontsize)
1216
g.save(dir+'/phi_cos_sum_1010_1026_15000.%s'%ext, axes=False, ymin=-50, fontsize=fontsize, typeset='latex')
1217
1218
def phi_interval_plot(xmin, xmax, zeros=1000,fontsize=12,drop=20):
1219
g = phi_approx_plot(xmin,xmax,zeros=zeros,fontsize=fontsize,drop=drop)
1220
g += line([(xmin,0),(xmax,0)],rgbcolor='black')
1221
return g
1222
1223
1224
def phi_approx(m, positive_only=False, symbolic=False):
1225
if symbolic:
1226
assert not positive_only
1227
s = var('s')
1228
return -sum(cos(log(s)*t) for t in zeta_zeros()[:m])
1229
1230
from math import cos, log
1231
v = [float(z) for z in zeta_zeros()[:m]]
1232
def f(s):
1233
s = log(float(s))
1234
return -sum(cos(s*t) for t in v)
1235
if positive_only:
1236
z = float(0)
1237
def g(s):
1238
return max(f(s),z)
1239
return g
1240
else:
1241
return f
1242
1243
def phi_approx_plot(xmin, xmax, zeros, pnts=2000, dots=True, positive_only=False,
1244
fontsize=7, drop=10, **args):
1245
phi = phi_approx(zeros, positive_only)
1246
g = plot(phi, xmin, xmax, alpha=0.7,
1247
plot_points=pnts, adaptive_recursion=0, **args)
1248
g.xmin(xmin); g.xmax(xmax)
1249
if dots: g += primepower_dots(xmin,xmax, fontsize=fontsize,drop=drop)
1250
return g
1251
1252
def primepower_dots(xmin, xmax, fontsize=7, drop=10):
1253
"""
1254
Return plot with dots at the prime powers in the given range.
1255
"""
1256
g = Graphics()
1257
for n in range(max(xmin,2),ceil(xmax)+1):
1258
F = factor(n)
1259
if len(F) == 1:
1260
g += point((n,0), pointsize=50*log(F[0][0]), rgbcolor=(1,0,0))
1261
if fontsize>0:
1262
g += text(str(n),(n,-drop),fontsize=fontsize, rgbcolor='black')
1263
g.xmin(xmin)
1264
g.xmax(xmax)
1265
return g
1266
1267
##############################################################
1268
# psi waves
1269
##############################################################
1270
1271
def fig_psi_waves(dir, ext):
1272
theta, t = var('theta, t')
1273
f = (theta*sin(t*theta) + 1/2 * cos(t*theta)) / (theta^2 + 1/4)
1274
g = plot(f(theta=zeta_zeros()[0]),(t,0,pi))
1275
g.save(dir + '/psi_just_waves1.%s'%ext, typeset='latex')
1276
1277
g += plot(f(theta=zeta_zeros()[1]),(t,0,pi), rgbcolor='red')
1278
g.save(dir + '/psi_2_waves.%s'%ext, typeset='latex')
1279
1280
f = (e^(t/2)*theta*sin(t*theta) + 1/2 * e^(t/2) * cos(t*theta))/(theta^2 + 1/4)
1281
g = plot(f(theta=zeta_zeros()[0]),(t,0,pi))
1282
g.save(dir + '/psi_with_first_zero.%s'%ext, typeset='latex')
1283
1284
g += plot(f(theta=zeta_zeros()[1]),(t,0,pi), rgbcolor='red')
1285
g.save(dir + '/psi_with_exp_2.%s'%ext, typeset='latex')
1286
1287
1288
##############################################################
1289
# Riemann's R_k
1290
##############################################################
1291
1292
def fig_moebius(dir,ext):
1293
g = plot(moebius,0, 50)
1294
g.save(dir+'/moebius.%s'%ext,figsize=[10,2], axes_pad=.1, fontsize=14, typeset='latex')
1295
1296
1297
def riemann_R(terms):
1298
c = [0] + [float(moebius(n))/n for n in range(1, terms+1)]
1299
def f(x):
1300
x = float(x)
1301
s = float(0)
1302
for n in range(1,terms+1):
1303
y = x^float(1/n)
1304
if y < 2:
1305
break
1306
s += c[n] * Li(y)
1307
return s
1308
return f
1309
1310
def plot_pi_riemann_gauss(xmin, xmax, terms):
1311
R = riemann_R(terms)
1312
g = plot(R, xmin, xmax)
1313
g += plot(prime_pi, xmin, xmax, rgbcolor='red')
1314
g += plot(Li, xmin, xmax, rgbcolor='purple')
1315
#x = var('x'); g += plot(x/(log(x)-1), xmin, xmax, rgbcolor='green')
1316
return g
1317
1318
def fig_pi_riemann_gauss(dir,ext):
1319
for m in [100,1000]:
1320
g = plot_pi_riemann_gauss(2,m, 100)
1321
g.save(dir+'/pi_riemann_gauss_%s.%s'%(m,ext), fontsize=20, typeset='latex')
1322
g = plot_pi_riemann_gauss(10000,11000, 100)
1323
g.save(dir +'/pi_riemann_gauss_10000-11000.%s'%ext, axes=False, frame=True, fontsize=20, typeset='latex')
1324
1325
1326
class RiemannPiApproximation:
1327
r"""
1328
Riemann's explicit formula for `\pi(X)`.
1329
1330
EXAMPLES::
1331
1332
We compute Riemann's analytic approximatin to `\pi(25)` using `R_{10}(x)`:
1333
1334
sage: R = RiemannPiApproximation(10, 100); R
1335
Riemann explicit formula for pi(x) for x <= 100 using R_k for k <= 10
1336
sage: R.Rk(100, 10)
1337
25.3364299527
1338
sage: prime_pi(100)
1339
25
1340
1341
"""
1342
def __init__(self, kmax, xmax, prec=50):
1343
"""
1344
INPUT:
1345
1346
- ``kmax`` -- (integer) large k allowed
1347
1348
- ``xmax`` -- (float) largest input x value allowed
1349
1350
- ``prec`` -- (default: 50) determines precision of certain series approximations
1351
1352
"""
1353
from math import log
1354
self.xmax = xmax
1355
self.kmax = kmax
1356
self.prec = prec
1357
self.N = int(log(xmax)/log(2))
1358
self.rho_k = [0] + [CDF(0.5, zeta_zeros()[k-1]) for k in range(1,kmax+1)]
1359
self.rho = [[0]+[rho_k / n for n in range(1, self.N+1)] for rho_k in self.rho_k]
1360
self.mu = [float(x) for x in moebius.range(0,self.N+2)]
1361
self.msum = sum([moebius(n) for n in xrange(1,self.N+1)])
1362
self._init_coeffs()
1363
1364
def __repr__(self):
1365
return "Riemann explicit formula for pi(x) for x <= %s using R_k for k <= %s"%(self.xmax, self.kmax)
1366
1367
def _init_coeffs(self):
1368
self.coeffs = [1]
1369
n_factorial = 1.0
1370
for n in xrange(1, self.prec):
1371
n_factorial *= n
1372
zeta_value = float(abs(zeta(n+1)))
1373
self.coeffs.append(float(1.0/(n_factorial*n*zeta_value)))
1374
1375
def _check(self, x, k):
1376
if x > self.xmax:
1377
raise ValueError, "x (=%s) must be at most %s"%(x, self.xmax)
1378
if k > self.kmax:
1379
raise ValueError, "k (=%s) must be at most %s"%(k, self.kmax)
1380
1381
@cached_method
1382
def R(self, x):
1383
from math import log
1384
y = log(x)
1385
z = y
1386
a = float(1)
1387
for n in xrange(1,self.prec):
1388
a += self.coeffs[n]*z
1389
z *= y
1390
return a
1391
1392
@cached_method
1393
def Rk(self, x, k):
1394
return self.R(x) + self.Sk(x, k)
1395
1396
@cached_method
1397
def Sk(self, x, k):
1398
"""
1399
Compute approximate correction term, so Rk(x,k) = R(x) + Sk(x,k)
1400
"""
1401
self._check(x, k)
1402
1403
from math import atan, pi, log
1404
log_x = log(x) # base e
1405
# This is from equation 32 on page 978 of Riesel-Gohl.
1406
term1 = self.msum / (2*log_x) + \
1407
(1/pi) * atan(pi/log_x)
1408
1409
# This is from equation 19 on page 975
1410
term2 = sum(self.Tk(x, v) for v in xrange(1,k+1))
1411
return term1 + term2
1412
1413
@cached_method
1414
def Tk(self, x, k):
1415
"""
1416
Compute sum from 1 to N of
1417
mu(n)/n * ( -2*sqrt(x) * cos(im(rho_k/n)*log(x) \
1418
- arg(rho_k/n)) / ( pi_over_2 * log(x) )
1419
"""
1420
self._check(x, k)
1421
x = float(x)
1422
log_x = log(x)
1423
val = float(0)
1424
rho_k = self.rho_k[k]
1425
rho = self.rho[k]
1426
for n in xrange(1, self.N+1):
1427
rho_k_over_n = rho[n]
1428
mu_n = self.mu[n]
1429
if mu_n != 0:
1430
z = Ei( rho_k_over_n * log_x)
1431
val += (mu_n/float(n)) * (2*z).real()
1432
return -val
1433
1434
def plot_Rk(self, k, xmin=2, xmax=None, **kwds):
1435
r"""
1436
Plot `\pi(x)` and `R_k` between ``xmin`` and ``xmax``. If `k`
1437
is a list, also plot every `R_k`, for `k` in the list.
1438
1439
The **kwds are passed onto the line function, which is used
1440
to draw the plot of `R_k`.
1441
"""
1442
if not xmax:
1443
xmax = self.xmax
1444
else:
1445
if xmax > self.xmax:
1446
raise ValueError, "xmax must be at most %s"%self.xmax
1447
xmax = min(self.xmax, xmax)
1448
if kwds.has_key('plot_points'):
1449
plot_points = kwds['plot_points']
1450
del kwds['plot_points']
1451
else:
1452
plot_points = 100
1453
eps = float(xmax-xmin)/plot_points
1454
if not isinstance(k, list):
1455
k = [k]
1456
f = sum(line([(x,self.Rk(x,kk)) for x in [xmin,xmin+eps,..,xmax]], **kwds)
1457
for kk in k)
1458
g = prime_pi.plot(xmin, xmax, rgbcolor='red')
1459
return g+f
1460
1461
1462
def fig_Rk(dir, ext):
1463
R = RiemannPiApproximation(50, 500, 50)
1464
for k,xmin,xmax in [(1,2,100), (10,2,100), (25,2,100),
1465
(50,2,100), (50,2,500) ]:
1466
print "plotting k=%s"%k
1467
g = R.plot_Rk(k, xmin, xmax, plot_points=300, thickness=0.65)
1468
g.save(dir + '/Rk_%s_%s_%s.%s'%(k,xmin,xmax,ext), typeset='latex')
1469
#
1470
g = R.plot_Rk(50, 350, 400, plot_points=200)
1471
g += plot(Li,350,400,rgbcolor='green')
1472
g.save(dir + '/Rk_50_350_400.%s'%ext, aspect_ratio=1, typeset='latex')
1473
1474
1475
##############################################################
1476
# Gaussian Primes
1477
##############################################################
1478
1479
def gaussian_primes(B):
1480
K.<i> = QuadraticField(-1)
1481
v = K.primes_of_bounded_norm(2*B*B)
1482
w = []
1483
for I in v:
1484
a = I.gens_reduced()[0]
1485
if abs(a[0])<=B and abs(a[1]) <= B:
1486
w.extend([a,-a,i*a,-i*a])
1487
w = [z for z in w if z[0]>0 and z[1]>=0]
1488
return w
1489
1490
def fig_gaussian_primes(dir, ext):
1491
B = 10
1492
w = gaussian_primes(B)
1493
points(w, pointsize=90, zorder=10).save("%s/gaussian_primes-%s.%s"%(dir,B,ext), aspect_ratio=1, gridlines=True, ticks=[[-1..B],[-1..B]], xmin=-.5, ymin=-.5, typeset='latex')
1494
B = 100
1495
w = gaussian_primes(B)
1496
p1 = points(w, pointsize=10, zorder=10)
1497
p2 = points(list(cartesian_product_iterator([[0..B],[0..B]])), pointsize=1, color='grey', zorder=15)
1498
(p1 + p2).save("%s/gaussian_primes-%s.%s"%(dir,B,ext), aspect_ratio=1, frame=True, axes=False, typeset='latex')
1499
1500
1501
1502
##############################################################
1503
# Random Walks
1504
##############################################################
1505
def random_walk(n):
1506
import random
1507
return stats.TimeSeries([random.choice([-1r,1r]) for _ in range(n)]).sums()
1508
1509
def random_walks(path, ext, B, n=1000, seed=1):
1510
set_random_seed(seed)
1511
path = path + '-%s'%B
1512
v = [random_walk(n) for i in range(B)]
1513
g = sum([z.plot(thickness=.3) for z in v])
1514
g.save(path + '.' + ext, fontsize=18, typeset='latex')
1515
s = sum([z.abs().vector()/B for z in v])
1516
avg = stats.TimeSeries(list(s))
1517
h = avg.plot() + plot(sqrt(2/pi)*sqrt(x), (0,g.xmax()), color='red', thickness=2)
1518
h.save(path + '-mean.' + ext, fontsize=18, typeset='latex')
1519
1520
def fig_random_walks(dir, ext):
1521
random_walks(dir + '/random_walks', ext, 3)
1522
random_walks(dir + '/random_walks', ext, 10)
1523
random_walks(dir + '/random_walks', ext, 100)
1524
random_walks(dir + '/random_walks', ext, 1000)
1525
1526
1527
##############################################################
1528
# Theta_C
1529
##############################################################
1530
1531
def plot_theta_C(C, xmax):
1532
theta = var('theta')
1533
T = SR(0)
1534
for q in prime_powers(C+1):
1535
if q > 1:
1536
p, n = factor(q)[0]
1537
T += 2 * p^(-n/2) * log(p) * cos(n*log(p)*theta)
1538
g = plot(T, 0, xmax)
1539
if g.ymax() > 10:
1540
g.ymax(10)
1541
ym = g.ymax()
1542
Z = []
1543
for y in zeta_zeros():
1544
if y > xmax: break
1545
Z.append(y)
1546
zeros = sum([arrow((x,ym),(x,0),rgbcolor='red',width=1.5,arrowsize=2)
1547
for i, x in enumerate(Z)])
1548
g += zeros
1549
g += plot(T.derivative(), 0, xmax, color='grey', alpha=.5)
1550
g.ymax(min(ym,10))
1551
g.ymin(max(-15,g.ymin()))
1552
return g
1553
1554
def fig_theta_C(dir, ext):
1555
def f(C,xmax):
1556
plot_theta_C(C,xmax).save(dir+'/theta_C-%s.%s'%(C,ext), fontsize=20, typeset='latex')
1557
f(3, 40)
1558
f(5, 40)
1559
f(10, 40)
1560
f(20, 40)
1561
f(100, 40)
1562
f(200, 40)
1563
f(500, 40)
1564
1565
def fig_theta_C_intro(dir, ext):
1566
theta = var('theta')
1567
def f(C):
1568
T = SR(0)
1569
for q in prime_powers(C+1):
1570
if q > 1:
1571
p, n = factor(q)[0]
1572
T += 2 * p^(-n/2) * log(p) * cos(n*log(p)*theta)
1573
return T
1574
h = f(3)
1575
plot(h, 5, 40).save("%s/theta_3_intro-1.%s"%(dir, ext), fontsize=18, typeset='latex')
1576
1577
roots = [(h.derivative().find_root(a,b),0) for a,b in [(5,7), (7,9), (10,12), (13,15), (16,17.8), (18,21), (21,24), (24.5,26), (27,30),
1578
(30, 33), (33,36), (36, 38), (38,42)]]
1579
1580
save(plot(h, 5, 40) + plot(h.derivative(),5,40,color='grey') + points(roots, color='red', pointsize=50, zorder=10),
1581
"%s/theta_3_intro-2.%s"%(dir, ext), fontsize=18)
1582
1583
1584
##############################################################
1585
# Cesaro Sum
1586
##############################################################
1587
1588
def fig_cesaro(dir, ext):
1589
theta = var('theta')
1590
def f(C):
1591
T = SR(0)
1592
for q in prime_powers(C+1):
1593
if q > 1:
1594
p, n = factor(q)[0]
1595
T += 2 * p^(-n/2) * log(p) * cos(n*log(p)*theta)
1596
return T
1597
k = f(5)^2
1598
T = var('T')
1599
assume(T>0)
1600
g(T) = 1/T * k.integrate(theta, 0, T)
1601
h = plot(g, 0, 40, color='red') + plot(k,0,40)
1602
h.save('%s/cesaro.%s'%(dir, ext), fontsize=16, typeset='latex')
1603
1604
1605
##############################################################
1606
# Spacing of Zeros
1607
##############################################################
1608
1609
def fig_zero_spacing(dir, ext):
1610
v = zeta_zeros()
1611
1612
def rmod(a,b):
1613
return a - a//b *b
1614
1615
def f(per):
1616
per = float(per)
1617
w = [rmod(a,per) for a in v]
1618
return stats.TimeSeries(w).plot_histogram()
1619
1620
f(2*pi).save('%s/zero-spacing-mod2pi.%s'%(dir, ext), fontsize=24, typeset='latex')
1621
1622
f(1).save('%s/zero-spacing-mod1.%s'%(dir, ext), fontsize=24, typeset='latex')
1623
1624
1625
##############################################################
1626
# Staircase of the Riemann spectrum
1627
##############################################################
1628
1629
def fig_staircase_riemann_spectrum(dir, ext):
1630
v = [float(a) for a in zeta_zeros()]
1631
t = stats.TimeSeries(v)
1632
def zeta_pi(X): # not efficient -- but lets us use general plotting machinery, and is fast enough for a book illustration!
1633
return len(t.clip_remove(max=X))
1634
def g(n, curve=False, **kwds):
1635
p = plot(zeta_pi, 1,n,
1636
plot_points=1000,rgbcolor='red',
1637
fillcolor=(.9,.9,.9),fill=True, **kwds)
1638
if curve:
1639
T = var('T')
1640
p += plot(T/(2*pi) * log(T/(2*pi*e)), 1, n, thickness=.5)
1641
return p
1642
g(30, curve=False, thickness=3).save('%s/staircase-riemann-spectrum-30.%s'%(dir, ext), fontsize=18, typeset='latex')
1643
g(50, curve=False, thickness=3).save('%s/staircase-riemann-spectrum-50.%s'%(dir, ext), fontsize=18, typeset='latex')
1644
g(100, curve=True, thickness=3).save('%s/staircase-riemann-spectrum-100.%s'%(dir, ext), fontsize=18, typeset='latex')
1645
g(1000, curve=True, thickness=3).save('%s/staircase-riemann-spectrum-1000.%s'%(dir, ext), fontsize=18, typeset='latex')
1646
1647
1648
##############################################################
1649
# Distribution of Riemann Spectrum Gaps
1650
##############################################################
1651
1652
def fig_riemann_spectrum_gaps(dir, ext):
1653
t = stats.TimeSeries(zeta_zeros())
1654
g = t.diffs().plot_histogram(bins=500, normalize=False)
1655
g.save("%s/riemann_spectrum_gaps.%s"%(dir, ext), xmax=2.3, frame=True, gridlines=True, axes=False, fontsize=18, typeset='latex')
1656
1657
1658
1659
##############################################################
1660
# Cesaro Sum of Li - Pi
1661
##############################################################
1662
1663
def prime_pi_time_series(B):
1664
"""
1665
Return the time series stats.TimeSeries([prime_pi(n) for n in [0..B-1]])
1666
but computed (vastly) more efficiently.
1667
"""
1668
x = 0
1669
w = []
1670
pp = 0
1671
for p in prime_range(B):
1672
w.extend([x]*(p-pp))
1673
pp = p
1674
x += 1
1675
w.extend([x]*(B-pp))
1676
return stats.TimeSeries(w)
1677
1678
def running_average(v):
1679
"""
1680
Return the running average of the time series... i.e., the Cesaro sum.
1681
i.e., stats.TimeSeries([0]+[v[:i].mean() for i in range(1,len(v))]),
1682
but computed (vastly) more efficiently.
1683
"""
1684
s = v.sums()
1685
# now divide s[i] by i+1 for i>=1.
1686
for i in range(1,len(s)):
1687
s[i] /= i+1
1688
return s
1689
1690
def fig_li_minus_pi(dir, ext):
1691
B = 250000
1692
v = prime_pi_time_series(B)
1693
from mpmath import li
1694
# This takes about 30 seconds...
1695
li_minus_pi = stats.TimeSeries([0,0] + [li(i,offset=True)-v[i] for i in range(2,B)])
1696
v2 = running_average(li_minus_pi)
1697
g = li_minus_pi.plot(plot_points=B) + v2.plot(color='red', plot_points=B) + plot(sqrt(2/pi)*sqrt(x/log(x)), (2, B), color='purple')
1698
g.save('%s/li-minus-pi-250000.%s'%(dir, ext), fontsize=18, ticks=[[10000, 100000, 250000], None], tick_formatter=['latex', None], typeset='latex')
1699
1700
1701
1702