Number Theory - Sage Education Days 3
system:sage


<h1 style="text-align: center;">Elementary Number Theory in Sage</h1>
<p style="text-align: center;"><span style="font-size: large;">Karl-Dieter Crisman, Gordon College</span></p>
<p style="text-align: center;"><span style="font-size: large;">Sage Education Days 3</span></p>
<p style="text-align: left;">&nbsp;</p>
<p>This&nbsp;<a href="http://www.sagemath.org" target="_blank">Sage</a>&nbsp;worksheet is partly based on materials developed for the MAA PREP Workshop "Sage: Using Open-Source Mathematics Software with Undergraduates" (funding provided by NSF DUE 0817071) and was developed for the Sage Education Days 3 as part of the UTMOST Grant (several NSF DUE grants).</p>
<h3>Why Sage and elementary number theory?</h3>
<p>Elementary number theory shows up in various places in the curriculum.</p>
<ul>
<li>"Proof Transitions"-type course</li>
<li>Discrete Mathematics</li>
<li>Abstract algebra</li>
<li>Number Theory (of course)</li>
<li>Cryptography</li>
<li>All of the above ideas in a non-major course</li>
</ul>
<p>So why use Sage (or any other computer system) in these contexts?</p>
<ul>
<li>It's not just about them seeing big numbers in action, though that was my original motivation.</li>
<li>It's about being able to take exploration 'by hand' to a different level.</li>
<li>It's about <em>visualizing</em> number theory in as many ways as possible.</li>
<li>It's about learning how to use a new tool with concepts we can still compute on pencil and paper.</li>
</ul>
<p>Since Sage began life as a project in algebraic and analytic number theory (and this continues to be a big emphasis), it is no surprise that functionality in this area is extremely comprehensive. &nbsp;</p>
<ul>
<li>But we don't have to be intimidated by that! &nbsp;</li>
<li>Instead, we should be glad that a free resource is truly state-of-the art.</li>
<li>Plus, it's also just <em>enjoyable</em> to explore elementary number theory by computer. &nbsp;</li>
</ul>
<p>The rest of this talk will show various basic topics that I cover in my own number theory course, highlighting how I use Sage to discuss them, especially to visualize them. &nbsp;Along the way, there will be plenty of '@interacts' or 'Sagelets' to (hopefully) inspire.</p>
<p>&nbsp;</p>
<p>One of the most important things in Sage is that the ring of integers modulo $n$ is available without further effort, so one can do examples in modular arithmetic very easily. &nbsp;</p>
<p>For instance, we can create a number in $\mathbb{Z}_n$. &nbsp;The 'type' command tells us that $a$ is not a regular integer.</p>

{{{id=109|
a = mod(2,11); a; type(a); a^10; a^1000000
///
2
<type 'sage.rings.finite_rings.integer_mod.IntegerMod_int'>
1
1
}}}

<p>In fact, "mod(a,n)" will create an integer modulo $n$, represented as the least nonnegative residue of that integer $a$.</p>
<p>Note two things about the previous cell.</p>
<ul>
<li>It just taught the student that the semicolon is a good way to separate Sage statements that are still all shown.</li>
<li>It verifies Fermat's "Little Theorem" for the prime $p=11$ and input $a=2$. &nbsp;</li>
</ul>
<p>This highlights two points about how to use Sage.</p>
<ul>
<li>A very basic way to use Sage is just as a big calculator for these kinds of things.</li>
<li><em>If you want to</em>, it's not hard to introduce (or reinforce) very basic programming concepts in a natural way.</li>
</ul>
<p>Here's another example of this. &nbsp;Recalling the basic programming construct called a <em>loop</em>, we can verify this for all $a$ in the integers modulo $p$. &nbsp;</p>
<ul>
<li>Here, instead of looping over an explicit list like '[0,1,2,3,...8,9,10]', we loop over a Sage object. &nbsp;You can do either.</li>
</ul>

{{{id=110|
for a in Integers(11):
    print a, a^10
///
0 0
1 1
2 1
3 1
4 1
5 1
6 1
7 1
8 1
9 1
10 1
}}}

<p>Notice that 'Integers(11)' gave us an algebraic object which is the ring of integers modulo the prime ideal generated by the element 11. &nbsp;This is handy later.</p>
<p>This works for much bigger numbers too, of course. &nbsp;We won't want to print out all the integers here, but we can check things, at least!</p>

{{{id=145|
p=random_prime(10^200,proof=True)
Zp=Integers(p) # Here we give ourselves shorthand for the modular integers
///
}}}

<p>If I'm going to do things with my integers modulo $p$, I should name them. &nbsp;I can name individual elements, too, as we saw above.</p>

{{{id=117|
a=Zp(2) # Here we ask for 2 as an element of that ring
a
p; a^(p-1)
///
12175749457485261323484342548813090831853594104614993416363962307255475872160716311272728940078211326045601051983589414266918386755653479471852846103058479931430533099115465341959637292476842332898283
1
}}}

<p>Here the student sees that it's the last line that counts when it comes to printing.</p>

{{{id=9|
a^(10^400)
///
667341860419940833796129083121003719719303707438654397250300236758579316812430945254713755146961572305624053487037532854247248289759595685250040511759445620464084950052304753065232782552186716875365
}}}

<p>A slight interlude:</p>
<p>Whenever you encounter a new object, you should try tab-completion to see what you can do to it. &nbsp; Let's try it here, and pick one (hopefully one that won't be too long!).</p>

{{{id=121|
Zp.
///
}}}

<p>Here's one that sounds interesting.</p>

{{{id=124|
Zp.zeta?
///
<html><!--notruncate-->

<div class="docstring">
    
  <p><strong>File:</strong> /home/sageserver/sage/devel/sage/sage/rings/ring.pyx</p>
<p><strong>Type:</strong> &lt;type &#8216;builtin_function_or_method&#8217;&gt;</p>
<p><strong>Definition:</strong> Zp.zeta(n=2, all=False)</p>
<p><strong>Docstring:</strong></p>
<blockquote>
<p>Return an n-th root of unity in self if there is one,
or raise an ArithmeticError otherwise.</p>
<p>INPUT:</p>
<ul class="simple">
<li><tt class="docutils literal"><span class="pre">n</span></tt> &#8211; positive integer</li>
<li><tt class="docutils literal"><span class="pre">all</span></tt> &#8211; bool, default: False.  If True, return a list of all n-th
roots of 1.</li>
</ul>
<p>OUTPUT:</p>
<p>element of self of finite order</p>
<p>EXAMPLES:</p>
<div class="highlight-python"><div class="highlight"><pre class="literal-block"><span class="gp">sage: </span><span class="n">QQ</span><span class="o">.</span><span class="n">zeta</span><span class="p">()</span>
<span class="go">-1</span>
<span class="gp">sage: </span><span class="n">QQ</span><span class="o">.</span><span class="n">zeta</span><span class="p">(</span><span class="mi">1</span><span class="p">)</span>
<span class="go">1</span>
<span class="gp">sage: </span><span class="n">CyclotomicField</span><span class="p">(</span><span class="mi">6</span><span class="p">)</span><span class="o">.</span><span class="n">zeta</span><span class="p">()</span>
<span class="go">zeta6</span>
<span class="gp">sage: </span><span class="n">CyclotomicField</span><span class="p">(</span><span class="mi">3</span><span class="p">)</span><span class="o">.</span><span class="n">zeta</span><span class="p">()</span>
<span class="go">zeta3</span>
<span class="gp">sage: </span><span class="n">CyclotomicField</span><span class="p">(</span><span class="mi">3</span><span class="p">)</span><span class="o">.</span><span class="n">zeta</span><span class="p">()</span><span class="o">.</span><span class="n">multiplicative_order</span><span class="p">()</span>
<span class="go">3</span>
<span class="gp">sage: </span><span class="n">a</span> <span class="o">=</span> <span class="n">GF</span><span class="p">(</span><span class="mi">7</span><span class="p">)</span><span class="o">.</span><span class="n">zeta</span><span class="p">();</span> <span class="n">a</span>
<span class="go">3</span>
<span class="gp">sage: </span><span class="n">a</span><span class="o">.</span><span class="n">multiplicative_order</span><span class="p">()</span>
<span class="go">6</span>
<span class="gp">sage: </span><span class="n">a</span> <span class="o">=</span> <span class="n">GF</span><span class="p">(</span><span class="mi">49</span><span class="p">,</span><span class="s">&#39;z&#39;</span><span class="p">)</span><span class="o">.</span><span class="n">zeta</span><span class="p">();</span> <span class="n">a</span>
<span class="go">z</span>
<span class="gp">sage: </span><span class="n">a</span><span class="o">.</span><span class="n">multiplicative_order</span><span class="p">()</span>
<span class="go">48</span>
<span class="gp">sage: </span><span class="n">a</span> <span class="o">=</span> <span class="n">GF</span><span class="p">(</span><span class="mi">49</span><span class="p">,</span><span class="s">&#39;z&#39;</span><span class="p">)</span><span class="o">.</span><span class="n">zeta</span><span class="p">(</span><span class="mi">2</span><span class="p">);</span> <span class="n">a</span>
<span class="go">6</span>
<span class="gp">sage: </span><span class="n">a</span><span class="o">.</span><span class="n">multiplicative_order</span><span class="p">()</span>
<span class="go">2</span>
<span class="gp">sage: </span><span class="n">QQ</span><span class="o">.</span><span class="n">zeta</span><span class="p">(</span><span class="mi">3</span><span class="p">)</span>
<span class="gt">Traceback (most recent call last):</span>
<span class="c">...</span>
<span class="nc">ValueError</span>: <span class="n-Identifier">no n-th root of unity in rational field</span>
<span class="gp">sage: </span><span class="n">Zp</span><span class="p">(</span><span class="mi">7</span><span class="p">,</span> <span class="n">prec</span><span class="o">=</span><span class="mi">8</span><span class="p">)</span><span class="o">.</span><span class="n">zeta</span><span class="p">()</span>
<span class="go">3 + 4*7 + 6*7^2 + 3*7^3 + 2*7^5 + 6*7^6 + 2*7^7 + O(7^8)</span>
</pre></div>
</div>
</blockquote>


</div>
</html>
}}}

<p>And we use it to find the fifth roots of unity in this field. &nbsp;</p>

{{{id=136|
root_list = Zp.zeta(5,all=True); root_list
///
[80199770388563324100334548626345240081294273289109866436996652525328268652922508892946068538641538316054373187019168781211876036849531337824832226216684677717580165592175377569174402189281574130719978, 69839783895572286297568834485025073557885364348071061715465477061873400359794989367423407683971299361817245213947182344090739843367076197016322541936552333837227080274674865687645877633828974738751695, 57407444219199061498801298672323590163238592201574572482836619025676869537007609315386800852204337587805249250896651467970585450518157701115893749407382500580682168292929753154872678880962261809848942, 41936641877539676652531567723249367917108638987131879521606243741814402095343724540844607212999297064816230828621064026263296602805535970091696570138772210094329631491849238240260893562065879302446837]
}}}

<p>Are these really fifth roots of unity? &nbsp;We use the&nbsp;<em>list comprehension</em>&nbsp;(set-builder) notation this time - another opportunity.</p>

{{{id=125|
[root^5 for root in root_list]
///
[1, 1, 1, 1]
}}}

<p>Luckily, it checked out. &nbsp;</p>
<ul>
<li><span style="font-size: small;">(If you didn't get any, then your random prime ended up being one for which there </span><em><span style="font-size: small;">are</span></em><span style="font-size: small;">&nbsp;no fifth roots of unity - try doing the sequence over again!)</span></li>
</ul>
<p><span style="font-size: medium;">I can also try to make this interactive, to remove the 'programmingese'.</span></p>

{{{id=16|
@interact
def _(update=''):
    p=random_prime(10^20,proof=True)
    Zp=Integers(p) 
    root_list = Zp.zeta(5,all=True)
    html("The prime is $%s$"%p)
    html("The roots are:")
    print root_list
    html("The fifth powers are:")
    print [root^5 for root in root_list]
///
}}}

<p>Incidentally, it seems like there is either a bug here or a documentation error. &nbsp;Shouldn't ONE be a fifth root of unity?﻿</p>
<p>This also exemplifies that it's important to know the level you are aiming at. &nbsp;</p>
<ul>
<li>Even list comprehensions could be challenging.</li>
<li>Will they be confused by long printed-out lists coming from list comprehensions, or from printing out a loop? &nbsp;This will vary a lot by the specific group of students and institution and prerequisites.</li>
</ul>
<p>Similarly, what background do they have mathematically?</p>
<ul>
<li>In my class, I expect them to have seen induction before and be somewhat able to use it.</li>
<li>But I don't expect algebra. &nbsp;I teach a very small amount about groups just so we can use that language, which is much more concise. &nbsp;We prove Lagrange's theorem about orders of elements in a very concrete way, and that's about it.</li>
<li>Some classes might have more background.</li>
<li>If you are just doing things in a proof transition course, you would change things correspondingly.&nbsp;</li>
</ul>
<p>&nbsp;</p>
<p>So the previous example may seem like it wasn't elementary, but let's think about the connection.</p>
<ul>
<li>A fifth root of unity is just a solution to $x^5=1$. &nbsp;</li>
<li>So these are solutions to $x^5\equiv 1\text{ mod }(p)$.</li>
</ul>
<p>And that <em>is</em>&nbsp;something we constantly should emphasize - that things are different modulo $p$!</p>
<p>There is a lot more we can access at the elementary level, such as</p>
<ul>
<li>primitive roots,&nbsp;</li>
<li>ways to write a number as a sum of squares,&nbsp;</li>
<li>Legendre symbols,</li>
<li>modular solving of basic equation,</li>
<li>etc. &nbsp;</li>
</ul>
<p>A good way to use Sage in this context is to allow students to experiment with pencil and paper first, then use Sage to see whether patterns they discover hold true before one attempts to prove them.</p>

{{{id=62|
p = 13
primitive_root(p); two_squares(p); is_prime(p); legendre_symbol(p,next_prime(p))
///
2
(2, 3)
True
1
}}}

<p>By the way, there are lots of pedagogical opportunities to improve Sage in number theory. &nbsp;For instance, someone just<a href="http://trac.sagemath.org/sage_trac/ticket/11138" target="_blank"> added the Jacobi symbol</a> to fill in the gap between Legendre and Kronecker symbols.</p>
<p>Another place where someone is working is in the 'naive' modular solver.</p>

{{{id=19|
var('x,y')
solve_mod([x^2 + 2 == x, x^2 + y == y^2], 14)
///
[(4, 2), (4, 6), (4, 9), (4, 13)]
}}}

<p>All of these capabilities makes it fairly easy to&nbsp;construct elementary cryptographic examples as well. &nbsp;</p>
<p>Here is a standard example of a Diffie-Hellman key exchange, for instance. &nbsp;</p>
<ul>
<li>Remember - if we don't do the second line, exponentiation would be impractical.</li>
</ul>

{{{id=128|
p=random_prime(10^20,10^30) # a random prime between these numbers
q=mod(primitive_root(p),p) # makes the primitive root a number modulo p, not an integer
n=randint(1,p) # Alice's random number
m=randint(1,p) # Bob's random number
x=q^n; y=q^m
p;q;x; y; x^m; y^n
///
17686401212520143419
2
11357069553300816796
6155949584214785623
7166503060918115829
7166503060918115829
}}}

<p>So Alice picks $n$ and sends $q^n$, Bob picks $m$ and sends $q^m$, and then they jointly have the key $$q^{mn}=q^{nm}$$ which an "Eve"sdropper can't discover (until they send lots of messages, at least).</p>
<ul>
<li><span style="font-size: small;">(We don't need $q$ to be a primitive root, but it's convenient to make sure it has the right order.)</span></li>
</ul>
<p>The last thing verifies that the private keys they get are the same.</p>
<p>&nbsp;</p>
<p>Most number-theoretic functions are pretty easy to find with tab-completion. &nbsp;Let's try this now - what's something we should have available?</p>

{{{id=22|

///
}}}

<p>For those who are interested, more advanced number-theoretic objects are easy to come by.</p>
<p>In the first example, K is the field extension $\QQ(\sqrt{-14})$, where the symbol 'a' plays the role of $\sqrt{-14}$; we discover several basic facts about $K$ in the next several cells.</p>

{{{id=111|
K.<a> = NumberField(x^2+14); K
///
Number Field in a with defining polynomial x^2 + 14
}}}

{{{id=113|
K.discriminant(); K.class_group().order(); K.class_group().is_cyclic()
///
-56
4
True
}}}

<p>Various zeta functions are also available; here is a complex plot of the Riemann zeta.</p>

{{{id=132|
complex_plot(zeta, (-30,30), (-30,30))
///
<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

<p>Sage supports various more advanced cryptographic procedures as well as some basic pedagogical ones natively. &nbsp;This example is adapted from the documentation. &nbsp;The pycrypto module has serious versions.</p>

{{{id=116|
from sage.crypto.block_cipher.sdes import SimplifiedDES
sdes = SimplifiedDES(); sdes
///
Simplified DES block cipher with 10-bit keys
}}}

{{{id=133|
bin = BinaryStrings()
P = [0,1,0,0,1,1,0,1] # our message
K = sdes.random_key() # generate a random key
C = sdes.encrypt(P, K) # encrypt our message
plaintxt = sdes.decrypt(C, K) # decrypt it
plaintxt # print it
///
[0, 1, 0, 0, 1, 1, 0, 1]
}}}

{{{id=26|

///
}}}

<p>&nbsp;</p>
<p>&nbsp;</p>
<p>Now it's time to start <em>visualizing</em> mathematics. First, let's see the traditional way. Look at powers of $a\text{ mod }(p)$; students are to look for patterns.</p>

{{{id=0|
@interact(layout=[['p','a']])
def _(p=(7,prime_range(50)),a=(3,[0..50])):
    b=mod(a,p)
    top=ceil(2*p/10)*10
    html("If we look at some of the powers of $%s$"%(a,))
    html("modulo the prime $%s$, we get:"%(p,))
    html("<ul>")
    for m in [0..top]:
        html("<li>$%s^{%s}\equiv %s\\text{ mod }(%s)"%(a,m,b^m,p))
    html("</ul>")
///
}}}

<p>Did you get a chance to come up with something? &nbsp;What do you think? &nbsp;Did you come up with some potential theorems? Notice that we use a little string formatting - totally optional, of course.</p>
<p>Now, how many theorems do you see here?</p>

{{{id=135|
p=11
P=matrix_plot(matrix(p-1,[mod(a,p)^b for a in range(1,p) for b in srange(p)]),cmap='jet')
show(P)
///
<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

<p>This is a graphic giving the various powers of integers modulo $p$ as colors, not numbers. &nbsp;The columns are the powers, so the first column is the zeroth power (always 1) and the second column gives the colors for the numbers modulo the given prime (first power).  The $a-1$ row and $b$ column gives the color corresponding to $(a-1)^b$ mod ($p$). &nbsp;That means the first ($0$th) column is the color for $1$ and the second ($1$th) column gives the colors of each element of $\mathbb{Z}_p$. &nbsp;For instance, $(3,4)$ corresponds to $2^4$ mod ($7$) in the initial example below.</p>
<p>Theorems:</p>
<ul>
<li>Fermat's Little</li>
<li>Square Root of Fermat's Little</li>
<li>Existence of Primitive Roots</li>
<li>All powers of $-1$ are $\pm 1$</li>
<li>More?</li>
</ul>
<p>It's certainly possible to have students or you now try changing $p$ to see what else happens. After all, they won't necessarily see any patterns immediately. &nbsp;Let's try this.</p>
<p>&nbsp;</p>
<p>But with just a little work, we can make this interactive. &nbsp;Now the programming is hidden, and it's far more interesting to the students. &nbsp;This is my pedagogy.</p>

{{{id=146|
@interact
def power_table_plot(p=(7,prime_range(50))):
    P=matrix_plot(matrix(p-1,[mod(a,p)^b for a in range(1,p) for b in srange(p)]),cmap='jet')
    show(P)
///
}}}

<p>My contention is that combining interaction with visualization is a great way to teach this discipline.</p>
<p>&nbsp;</p><p>Note that if you don't like the colors, you can change the word in the quotes after the word 'cmap'; if you get rid of that, it will be a grayscale plot. &nbsp;Some others you could try are 'Oranges' or &nbsp;'hsv' or $\ldots$&nbsp;</p>
<p>&nbsp;</p>


<p>Another nice exercise is to start trying to find patterns in the primes $p$ such that &nbsp;$$x^2\equiv -1\text{ mod }(p)$$ has a solution. &nbsp;Students would start by hand, but then we create a table and see if we can find something. &nbsp; Here, the visualization is more traditional - via a table - but still interactive, with immediate feedback.</p>

{{{id=31|
@interact
def _(n=20):
    yeslist=[]
    nolist=[]
    for p in prime_range(3,n):
        res = 0
        for res in [0..p]:
            if mod(res,p)^2+1 == 0:
                yeslist.append(p)
                break
        else:
            nolist.append(p)
    t = [['exist', 'do not exist']] + [[a,b] for (a,b) in map(None,yeslist,nolist)]
    for item in t:
        for i in range(len(item)):
            if item[i] is None:
                item[i]=''
    html("Solutions to $x^2\equiv -1$ mod $(p)$ for $2< p <%s$:"%n)
    html.table(t, header = True)
///
}}}

<p>Another very useful object is the prime counting function $\pi(x)$. &nbsp;This comes with its own custom plotting.</p>

{{{id=119|
prime_pi(100); plot(prime_pi,1,100)
///
25
<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

<p>A very nice aspect of Sage is combining several aspects of mathematics together. &nbsp;It can be very eye-opening to students to see analytic aspects of number theory early on. &nbsp;</p>
<ul>
<li><span style="font-size: small;">(Note that we have to reassign $x$ to a variable, since above it was a cryptographic key!)</span></li>
</ul>

{{{id=127|
var('x')
plot(prime_pi,2,10^6,thickness=2)+plot(Li,2,10^6,color='red')+plot(x/ln(x),2,10^6,color='green')
///
<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

<p>We'll return to this later.</p>

{{{id=138|

///
}}}

<p>&nbsp;</p>
<p>I think that one of the most powerful ways to use Sage is to illuminate visual proofs that are often hard to understand with a single static picture. &nbsp;The following is part of the proof of quadratic reciprocity.</p>
<p>&nbsp;</p>
<p>Recall the statement of quadratic reciprocity. &nbsp;For odd primes $p$ and $q$, we have that $$\left(\frac{p}{q}\right)\left(\frac{q}{p}\right)=(-1)^{\frac{p-1}{2}\frac{q-1}{2}}$$ &nbsp;That is to say, the Legendre symbols are the same unless $p$ and $q$ are both of the form $4k+3$.</p>
<p><strong>Proof</strong>:</p>
<p>We will use a criterion of Eisenstein's that I introduced in class (see <a href="http://sofia.nmsu.edu/~history/eisenstein/eisenstein.html" target="_blank">this article</a>). &nbsp;</p>
<p>Recall that $p$ and $q$ are odd primes. Now let &nbsp;$$R=\sum_{\text{even }e,\; 0&lt;e&lt;p}\left\lfloor\frac{qe}{p}\right\rfloor\; $$ be the exponent in question, so that $$\left(\frac{q}{p}\right)=(-1)^R\; .$$ &nbsp;</p>
<p>&nbsp;The key will be <em>geometrically</em> interpreting $\left\lfloor\frac{qe}{p}\right\rfloor$ as being the biggest integer less than $\frac{qe}{p}$. &nbsp;The following features are present in the next graphic, which should clarify how we'll think of it geometrically.</p>
<ul>
<li>The line through the origin with slope $q/p$ (dotted blue).</li>
<li>All the grid points in the box of width $p$ and height $q$ (box red, points black).</li>
<li>Points with <em>even</em> $x$-coordinate corresponding to the highest one can get but stay under the line of slope $q/p$ (points blue).</li>
<li>The box of width $\frac{p-1}{2}$ and height $\frac{q-1}{2}$ (green), which we'll need in a moment.</li>
</ul>
<p>It should be clear that each blue stack has the same height as $\left\lfloor\frac{qe}{p}\right\rfloor$ for some even $e$.</p>

{{{id=1|
var('x,y')
@interact
def _(p=(11,prime_range(3,100)),q=(7,prime_range(3,100))):
    E = [2,4..p-1]
    plot4 = plot((q/p)*x,(x,0,p),linestyle='--')
    plot3 = line([[0,0],[p,0],[p,q],[0,q],[0,0]],rgbcolor=(1,0,0))
    plot2 = line([[0,0],[(p-1)/2,0],[(p-1)/2,(q-1)/2],[0,(q-1)/2],[0,0]],color='green')
    grid_pts_1 = [[i,j] for i in [1..p] for j in [1..q]]
    grid_pts_2 = [[i,j] for i in [1..(p-1)/2] for j in [1..(q-1)/2]]
    plot_grid_pts = points(grid_pts_1,rgbcolor=(0,0,0),pointsize=10)
    lattice_pts1 = [coords for coords in grid_pts_1 if (coords[0]*q-coords[1]*p>0 and coords[0]<p and coords[0] in E)]
    lattice_pts2 = [coords for coords in grid_pts_2 if (coords[0]*q-coords[1]*p>0 and coords[0]>p/2)]
    num1, num2 = len(lattice_pts1), len(lattice_pts2)
    if len(lattice_pts1)!=0:
        plot_lattice_pts1 = points(lattice_pts1, rgbcolor = (0,0,1),pointsize=20)
    else:
        plot_lattice_pts1 = Graphics()
    if len(lattice_pts2)!=0:
        plot_lattice_pts2 = points(lattice_pts2, rgbcolor = (0,.5,0),pointsize=20)
    else:
        plot_lattice_pts2 = Graphics()
    show(plot2+plot3+plot4+plot_grid_pts+plot_lattice_pts1,xmax=p,ymax=q,ymin=0)
    forms = '$'+'+'.join(['\left\lfloor\\frac{%s\cdot %s}{%s}\\right\\rfloor'%(q,e,p) for e in E])+'$'
    html("The blue dots represent "+forms)
    forms2 = '$'+'+'.join(['\left\lfloor\\frac{%s}{%s}\\right\\rfloor'%(q*e,p) for e in E])
    forms3 = '+'.join(['%s'%(floor(q*e/p)) for e in E])+'=%s\equiv%s\\text{ mod }(2)$'%(sum([floor(q*e/p) for e in E]),sum([floor(q*e/p) for e in E])%2)
    html("This simplifies to "+forms2+'='+forms3)
///
}}}

<p>What we will now do is try to convince ourselves that the number of blue points has the same parity as the <em>total</em> number of (positive) points in and on the green box under the dotted line. &nbsp; If we can do that, the following steps finish the proof of quadratic reciprocity.</p>
<ul>
<li>First, to get $\left(\frac{q}{p}\right)$, we can ignore $R$ and just focus on the number of positive lattice points in that more-or-less triangle.</li>
<li>Then the same argument applies to $\left(\frac{p}{q}\right)$; we could ignore its crazy exponent (call it $R'$, if you must) and instead focus on the number of positive lattice points in and on the green box <em>to the right of</em>&nbsp;the dotted line. &nbsp;(We've switched the role of $x$ and $y$, so 'above the $y$-axis' means to its right.)</li>
<li>So the exponent of $-1$ we expect from $\left(\frac{q}{p}\right)\left(\frac{p}{q}\right)$ is just the sum of those two amounts. &nbsp;That is, the exponent is the number of points in and on the green box.</li>
<li>You'll note that this box has dimensions $\frac{p-1}{2}$ and $\frac{q-1}{2}$, so that would mean $$\frac{p-1}{2}\cdot\frac{q-1}{2}\equiv \sum_{\text{even }e,\; 0&lt;e&lt;p}\left\lfloor\frac{qe}{p}\right\rfloor+\sum_{\text{even }f,\; 0&lt;f&lt;q}\left\lfloor\frac{pf}{q}\right\rfloor\text{ mod }(2)\; ,$$ so that $$\left(\frac{q}{p}\right)\left(\frac{p}{q}\right)=(-1)^{R+R'}=(-1)^{\frac{p-1}{2}\cdot\frac{q-1}{2}}\; .$$&nbsp;</li>
</ul>
<p>This would definitely complete the theorem! &nbsp;And I prove this in the rest of class.</p>
<p>&nbsp;</p>
<p>Why is this better?&nbsp;</p>
<ul>
<li>Many proofs of QR are visual. &nbsp;But there is <strong>one</strong>&nbsp;picture in the book, which is very hard to interpret for students.</li>
<li>By having color and interaction, students can really see what changes, and how it works.</li>
<li>This leads to students not simply mechanically <em>verifying</em>&nbsp;the proofs, but actually understanding the proof as a whole.</li>
<li>And since (for me, at least) I do not expect students to repeat proofs of these difficult theorems, now we have at least achieved them understanding a big idea in mathematics.</li>
</ul>
<p>&nbsp;</p>
<p>&nbsp;</p>
<p>Analytic number theory is another great place to look - which makes sense, because a lot of what it does is make computations related to L-functions... &nbsp;here's another excerpt, with some final graphics. &nbsp;Of course, this is toward the end of class. &nbsp;</p>
<p>&nbsp;</p>
<p>Somewhat remarkably, the first people we know of compiling substantial data about this are Gauss and Legendre. &nbsp;See the handout.</p>
<ul>
<li>Legendre tries to estimate that $\pi(x)\approx \frac{x}{\ln(x)-A}$, where he fudges the constant $A\approx 1.08366$.
<ul>
<li>In fact, he says something more precise - that $\pi(x)$ is <em>asymptotic</em> to this function. &nbsp;That is, he suggested that $$\lim_{x\to\infty}\frac{\pi(x)}{x/(\ln(x)-A)}=1$$</li>
<li>Another way to think of this is that the average chance of a number of size $x$ being prime is about $\frac{1}{\ln(x)-A}$.&nbsp;</li>
</ul>
</li>
<li>Naturally, Gauss came up with a solution that was both more elegant and correct. &nbsp;<em>And he didn't tell anyone for over fifty years!</em>
<ul>
<li>Gauss' conjecture was that $$\lim_{x\to\infty}\frac{\pi(x)}{x/\ln(x)}=1$$ - that $\pi(x)$ is asymptotic to $\frac{x}{\ln(x)}$.</li>
<li>So, roughly $1/\ln(x)$ integers near $x$ would be prime.</li>
</ul>
</li>
<li>In fact, Gauss makes this more precise. &nbsp;
<ul>
<li>If the probability density function is that $1/\ln(x)$ integers near $x$ are prime, then we should do as we always do with such functions, and <em>integrate</em>&nbsp;the function to get the cumulative amount!</li>
<li>That is, we should expect that $\pi(x)\approx \int_0^x\frac{dt}{\ln(t)}$.</li>
<li>Or, that $$\lim_{x\to\infty}\pi(x)/\int_0^x\frac{dt}{\ln(t)}=1\; .$$</li>
</ul>
</li>
</ul>
<p>This should sound 100% crazy! &nbsp;But Gauss and Legendre were no fools, and the accuracy is astounding. &nbsp;Let's call $Li(x)=\int_0^\infty \frac{dt}{\ln(t)}$ (yes, this is a convergent integral).</p>

{{{id=27|
@interact
def _(n=100):
    show(plot(prime_pi,3,n,color='black',legend_label='$\pi(x)$')+plot(x/ln(x),3,n,color='red',legend_label='$x/\ln(x)$')+plot(Li,3,n,color='green',legend_label='$Li(x)$'))
///
}}}

<p>Notice how much closer $Li(x)$ is than the $x/\ln(x)$ estimate. &nbsp;It's usually closer by several orders of magnitude.</p>

{{{id=28|
@interact
def _(n=[100,1000,1000000,1000000000]):
    P = prime_pi(n)
    html("$\pi(%s)=%s$"%(n,prime_pi(n)))
    html("The error with $%s/\ln(%s)$ is $\\approx %s$"%(n,n,P-(n/ln(n)).n()))
    html("The error with $Li(%s)$ is $\\approx %s$"%(n,P-Li(n)))
///
}}}

<p>Here's an even better estimate I talk about the last day or two, while students are working on take-home projects and don't want homework in addition.</p>

{{{id=139|
@interact
def _(n=(1000,(1000,10^6)),k=(3,[1..10])):
   P = plot(prime_pi,n-1000,n,color='black',legend_label='$\pi(x)$')
   P += plot(Li,n-1000,n,color='green',legend_label='$Li(x)$')
   F = lambda x: sum([Li(x^(1/j))*moebius(j)/j for j in [1..k]])
   P += plot(lambda x: Li(x)-.5*Li(sqrt(x)),n-1000,n,color='red',legend_label='$Li(x)-\\frac{1}{2}Li(\sqrt{x})$')
   P += plot(F,n-1000,n,color='blue',legend_label='$\sum_{j=1}^{%s}\\frac{\mu(j)}{j}Li(x^{1/j})$'%k)
   show(P)
///
}}}

{{{id=140|

///
}}}

<p>And here we see Riemann's exact formula for $\pi(x)$, through the zeros on the critical line up to $t=150$.</p>

{{{id=141|
import mpmath
L = lcalc.zeros_in_interval(10,150,0.1)
n=100
P = plot(prime_pi,n-50,n,color='black',legend_label='$\pi(x)$')
P += plot(Li,n-50,n,color='green',legend_label='$Li(x)$')
G = lambda x: sum([mpmath.li(x^(1/j))*moebius(j)/j for j in [1..3]])
P += plot(G,n-50,n,color='red',legend_label='$\sum_{j=1}^{%s}\\frac{\mu(j)}{j}Li(x^{1/j})$'%3)
F = lambda x: sum([(mpmath.li(x^(1/j))-log(2)+numerical_integral(1/(y*(y^2-1)*log(y)),x^(1/j),oo)[0])*moebius(j)/j for j in [1..3]])-sum([(mpmath.ei(log(x)*((0.5+l[0]*i)/j))+mpmath.ei(log(x)*((0.5-l[0]*i)/j))).real for l in L for j in [1..3]])
P += plot(F,n-50,n,color='blue',legend_label='Really good estimate',plot_points=50)
show(P)
///
<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

<p>That's just for fun, really. &nbsp;Thank you!</p>

{{{id=148|

///
}}}