CoCalc Public Filesbook / rh / code / code.sage
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
8def 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##############################################################
28def 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
49class 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
212def 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##############################################################
231##############################################################
232def fig_questions(dir, ext):
233    g = questions(100,17,20)
234    g.save(dir + '/questions.%s'%ext, axes=False, typeset='latex')
235
236def 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
248def 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
256def 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
298def 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
339def 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
344def 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
370def 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
374def 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
379def 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
407def 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
436def 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
449def 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
466def 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
474def 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
479for B in [1..8]:
480    f(B)
481"""
482
483
484def 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
506def 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##############################################################
528def 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
556def 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
566def 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
598def 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##############################################################
636def 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
648def 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
656def 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
667def 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
675def 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
684def 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
693def 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
701def 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)
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
747    s1 = plot_sieve(1, m, poly={'rgbcolor':(.85,.9,.7)},
749    s2 = plot_sieve(2, m, poly={'rgbcolor':(.75,.8,.6)},
751    s3 = plot_sieve(0, m, poly={'rgbcolor':(1,.7,.5)},
753    return s1+s2+s3
754
755def 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},
774        g += s
775    return g
776
778    P = [1] + prime_range(int(sqrt(x))+1) + [0]
780    return G
781
782
783##############################################################
784# Area under plot of 1/log(x)
785##############################################################
786
787#auto
788def 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
817def 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##############################################################
827
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
840def 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
854def 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
864def plot_psi(xmax, **kwds):
865    v = psi_data(xmax)
866    return plot_step_function(v, **kwds)
867
868def 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
876def 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
881def 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
887def 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
907def 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
939def 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
988def 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
995def 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
1001def 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##############################################################
1009def 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##############################################################
1019def 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##############################################################
1028def 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
1036def 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
1043def 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
1057def 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)
1066
1067
1068##############################################################
1069# Distribution section
1070##############################################################
1071
1072def 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
1082def 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
1086def 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
1092def 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
1098def 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
1107def 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
1125def 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
1130def 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
1147def 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
1167def 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
1186def 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
1193def 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
1207def 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
1218def 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
1224def 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
1243def 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,
1248    g.xmin(xmin); g.xmax(xmax)
1249    if dots: g += primepower_dots(xmin,xmax, fontsize=fontsize,drop=drop)
1250    return g
1251
1252def 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
1271def 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
1292def fig_moebius(dir,ext):
1293    g = plot(moebius,0, 50)
1295
1296
1297def 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
1310def 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
1318def 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
1326class 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
1462def 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
1479def gaussian_primes(B):
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
1490def 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##############################################################
1505def random_walk(n):
1506    import random
1507    return stats.TimeSeries([random.choice([-1r,1r]) for _ in range(n)]).sums()
1508
1509def 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
1520def 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
1531def 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
1554def 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
1565def 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
1588def 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
1609def 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
1629def 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
1652def 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
1663def 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
1678def 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
1690def 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