Elementary Number Theory in Sage

In fact, "mod(a,n)" will create an integer modulo $n$, represented as the least nonnegative residue of that integer $a$.

Note two things about the previous cell.

• It just taught the student that the semicolon is a good way to separate Sage statements that are still all shown.
• It verifies Fermat's "Little Theorem" for the prime $p=11$ and input $a=2$.  

This highlights two points about how to use Sage.

• A very basic way to use Sage is just as a big calculator for these kinds of things.
• If you want to, it's not hard to introduce (or reinforce) very basic programming concepts in a natural way.

Here's another example of this.  Recalling the basic programming construct called a loop, we can verify this for all $a$ in the integers modulo $p$.  

• Here, instead of looping over an explicit list like '[0,1,2,3,...8,9,10]', we loop over a Sage object.  You can do either.

Notice that 'Integers(11)' gave us an algebraic object which is the ring of integers modulo the prime ideal generated by the element 11.  This is handy later.

This works for much bigger numbers too, of course.  We won't want to print out all the integers here, but we can check things, at least!

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

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

A slight interlude:

Whenever you encounter a new object, you should try tab-completion to see what you can do to it.   Let's try it here, and pick one (hopefully one that won't be too long!).

Here's one that sounds interesting.

And we use it to find the fifth roots of unity in this field.  

Are these really fifth roots of unity?  We use the list comprehension (set-builder) notation this time - another opportunity.

Luckily, it checked out.  

• (If you didn't get any, then your random prime ended up being one for which there are no fifth roots of unity - try doing the sequence over again!)

I can also try to make this interactive, to remove the 'programmingese'.

Incidentally, it seems like there is either a bug here or a documentation error.  Shouldn't ONE be a fifth root of unity?

This also exemplifies that it's important to know the level you are aiming at.  

• Even list comprehensions could be challenging.
• Will they be confused by long printed-out lists coming from list comprehensions, or from printing out a loop?  This will vary a lot by the specific group of students and institution and prerequisites.

Similarly, what background do they have mathematically?

• In my class, I expect them to have seen induction before and be somewhat able to use it.
• But I don't expect algebra.  I teach a very small amount about groups just so we can use that language, which is much more concise.  We prove Lagrange's theorem about orders of elements in a very concrete way, and that's about it.
• Some classes might have more background.
• If you are just doing things in a proof transition course, you would change things correspondingly. 

So the previous example may seem like it wasn't elementary, but let's think about the connection.

• A fifth root of unity is just a solution to $x^5=1$.  
• So these are solutions to $x^5\equiv 1\text{ mod }(p)$.

And that is something we constantly should emphasize - that things are different modulo $p$!

There is a lot more we can access at the elementary level, such as

• primitive roots, 
• ways to write a number as a sum of squares, 
• Legendre symbols,
• modular solving of basic equation,
• etc.  

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.

By the way, there are lots of pedagogical opportunities to improve Sage in number theory.  For instance, someone just added the Jacobi symbol to fill in the gap between Legendre and Kronecker symbols.

Another place where someone is working is in the 'naive' modular solver.

All of these capabilities makes it fairly easy to construct elementary cryptographic examples as well.  

Here is a standard example of a Diffie-Hellman key exchange, for instance.  

• Remember - if we don't do the second line, exponentiation would be impractical.

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).

• (We don't need $q$ to be a primitive root, but it's convenient to make sure it has the right order.)

The last thing verifies that the private keys they get are the same.

Most number-theoretic functions are pretty easy to find with tab-completion.  Let's try this now - what's something we should have available?

For those who are interested, more advanced number-theoretic objects are easy to come by.

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.

= 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 }}}

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

{{{id=132| complex_plot(zeta, (-30,30), (-30,30)) /// }}}

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

{{{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| /// }}}

 

 

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

{{{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("

") for m in [0..top]: html("
• $%s^{%s}\equiv %s\\text{ mod }(%s)"%(a,m,b^m,p)) html("

") /// }}}

Did you get a chance to come up with something?  What do you think?  Did you come up with some potential theorems? Notice that we use a little string formatting - totally optional, of course.

Now, how many theorems do you see here?

{{{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) /// }}}

This is a graphic giving the various powers of integers modulo $p$ as colors, not numbers.  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$).  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$.  For instance, $(3,4)$ corresponds to $2^4$ mod ($7$) in the initial example below.

Theorems:

• Fermat's Little
• Square Root of Fermat's Little
• Existence of Primitive Roots
• All powers of $-1$ are $\pm 1$
• More?

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.  Let's try this.

 

But with just a little work, we can make this interactive.  Now the programming is hidden, and it's far more interesting to the students.  This is my pedagogy.

{{{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) /// }}}

My contention is that combining interaction with visualization is a great way to teach this discipline.

 

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.  Some others you could try are 'Oranges' or  'hsv' or $\ldots$ 

 

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

{{{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) /// }}}

Another very useful object is the prime counting function $\pi(x)$.  This comes with its own custom plotting.

{{{id=119| prime_pi(100); plot(prime_pi,1,100) /// 25 }}}

A very nice aspect of Sage is combining several aspects of mathematics together.  It can be very eye-opening to students to see analytic aspects of number theory early on.  

• (Note that we have to reassign $x$ to a variable, since above it was a cryptographic key!)

{{{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') /// }}}

We'll return to this later.

{{{id=138| /// }}}

 

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.  The following is part of the proof of quadratic reciprocity.

 

Recall the statement of quadratic reciprocity.  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}}$$  That is to say, the Legendre symbols are the same unless $p$ and $q$ are both of the form $4k+3$.

Proof:

We will use a criterion of Eisenstein's that I introduced in class (see this article).  

Recall that $p$ and $q$ are odd primes. Now let  $$R=\sum_{\text{even }e,\; 0<e<p}\left\lfloor\frac{qe}{p}\right\rfloor\; $$ be the exponent in question, so that $$\left(\frac{q}{p}\right)=(-1)^R\; .$$  

 The key will be geometrically interpreting $\left\lfloor\frac{qe}{p}\right\rfloor$ as being the biggest integer less than $\frac{qe}{p}$.  The following features are present in the next graphic, which should clarify how we'll think of it geometrically.

• The line through the origin with slope $q/p$ (dotted blue).
• All the grid points in the box of width $p$ and height $q$ (box red, points black).
• Points with even $x$-coordinate corresponding to the highest one can get but stay under the line of slope $q/p$ (points blue).
• The box of width $\frac{p-1}{2}$ and height $\frac{q-1}{2}$ (green), which we'll need in a moment.

It should be clear that each blue stack has the same height as $\left\lfloor\frac{qe}{p}\right\rfloor$ for some even $e$.

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) /// }}}

What we will now do is try to convince ourselves that the number of blue points has the same parity as the total number of (positive) points in and on the green box under the dotted line.   If we can do that, the following steps finish the proof of quadratic reciprocity.

• 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.
• 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 to the right of the dotted line.  (We've switched the role of $x$ and $y$, so 'above the $y$-axis' means to its right.)
• 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.  That is, the exponent is the number of points in and on the green box.
• 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<e<p}\left\lfloor\frac{qe}{p}\right\rfloor+\sum_{\text{even }f,\; 0<f<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}}\; .$$ 

This would definitely complete the theorem!  And I prove this in the rest of class.

Why is this better? 

• Many proofs of QR are visual.  But there is one picture in the book, which is very hard to interpret for students.
• By having color and interaction, students can really see what changes, and how it works.
• This leads to students not simply mechanically verifying the proofs, but actually understanding the proof as a whole.
• 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.

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...  here's another excerpt, with some final graphics.  Of course, this is toward the end of class.  

Somewhat remarkably, the first people we know of compiling substantial data about this are Gauss and Legendre.  See the handout.

• Legendre tries to estimate that $\pi(x)\approx \frac{x}{\ln(x)-A}$, where he fudges the constant $A\approx 1.08366$.
  • In fact, he says something more precise - that $\pi(x)$ is asymptotic to this function.  That is, he suggested that $$\lim_{x\to\infty}\frac{\pi(x)}{x/(\ln(x)-A)}=1$$
  • 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}$. 
• Naturally, Gauss came up with a solution that was both more elegant and correct.  And he didn't tell anyone for over fifty years!
  • 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)}$.
  • So, roughly $1/\ln(x)$ integers near $x$ would be prime.
• In fact, Gauss makes this more precise.  
  • 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 integrate the function to get the cumulative amount!
  • That is, we should expect that $\pi(x)\approx \int_0^x\frac{dt}{\ln(t)}$.
  • Or, that $$\lim_{x\to\infty}\pi(x)/\int_0^x\frac{dt}{\ln(t)}=1\; .$$

This should sound 100% crazy!  But Gauss and Legendre were no fools, and the accuracy is astounding.  Let's call $Li(x)=\int_0^\infty \frac{dt}{\ln(t)}$ (yes, this is a convergent integral).

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

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.

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

That's just for fun, really.  Thank you!