Sage Quickstart for Number Theory

This Sage quickstart tutorial was developed for the MAA PREP Workshop “Sage: Using Open-Source Mathematics Software with Undergraduates” (funding provided by NSF DUE 0817071). It is licensed under the Creative Commons Attribution-ShareAlike 3.0 license (CC BY-SA).

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. And it’s also just enjoyable to explore elementary number theory by computer.

Modular Arithmetic

Conveniently, the ring of integers modulo \(n\) is always available in Sage, so we can do modular arithmetic very easily. For instance, we can create a number in \(\ZZ/11\ZZ\). The type command tells us that \(a\) is not a regular integer.

sage: a = mod(2,11); a; type(a); a^10; a^1000000
2
<class 'sage.rings.finite_rings.integer_mod.IntegerMod_int'>
1
1

>>> from sage.all import *
>>> a = mod(Integer(2),Integer(11)); a; type(a); a**Integer(10); a**Integer(1000000)
2
<class 'sage.rings.finite_rings.integer_mod.IntegerMod_int'>
1
1

a = mod(2,11); a; type(a); a^10; a^1000000

Note that we verified Fermat’s “Little Theorem” in the previous cell, for the prime \(p=11\) and input \(a=2\).

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.

sage: for a in Integers(11):
....:     print("{} {}".format(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

>>> from sage.all import *
>>> for a in Integers(Integer(11)):
...     print("{} {}".format(a, a**Integer(10)))
0 0
1 1
2 1
3 1
4 1
5 1
6 1
7 1
8 1
9 1
10 1

for a in Integers(11):
    print("{} {}".format(a, a^10))

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 works for much bigger numbers too, of course.

sage: p=random_prime(10^200,proof=True)
sage: Zp=Integers(p) # Here we give ourselves shorthand for the modular integers
sage: a=Zp(2) # Here we ask for 2 as an element of that ring
sage: p; a; a^(p-1); a^(10^400)
83127880126958116183078749835647757239842289608629126718968330784897646881689610705533628095938824110164366160161355539845499311180100402016248362566463907409939681883876411550651284088712896660589151
2
1
21428062940380156097538386722913390559445420264965409656232714664156136144920173818936415427844953758774658350253363113516712541610660591925149144205368271119123211091215746697984955519927521190733305

>>> from sage.all import *
>>> p=random_prime(Integer(10)**Integer(200),proof=True)
>>> Zp=Integers(p) # Here we give ourselves shorthand for the modular integers
>>> a=Zp(Integer(2)) # Here we ask for 2 as an element of that ring
>>> p; a; a**(p-Integer(1)); a**(Integer(10)**Integer(400))
83127880126958116183078749835647757239842289608629126718968330784897646881689610705533628095938824110164366160161355539845499311180100402016248362566463907409939681883876411550651284088712896660589151
2
1
21428062940380156097538386722913390559445420264965409656232714664156136144920173818936415427844953758774658350253363113516712541610660591925149144205368271119123211091215746697984955519927521190733305

p=random_prime(10^200,proof=True)
Zp=Integers(p) # Here we give ourselves shorthand for the modular integers
a=Zp(2) # Here we ask for 2 as an element of that ring
p; a; a^(p-1); a^(10^400)

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

sage: Zp.

>>> from sage.all import *
>>> Zp.

Zp.

Here’s one that sounds interesting.

sage: Zp.zeta?

>>> from sage.all import *
>>> Zp.zeta?

Zp.zeta?

And we use it to find the fifth roots of unity in this field. We use the list comprehension (set-builder) notation from the programming tutorial this time.

sage: root_list = Zp.zeta(5,all=True); root_list
[80199770388563324100334548626345240081294273289109866436996652525328268652922508892946068538641538316054373187019168781211876036849531337824832226216684677717580165592175377569174402189281574130719978, 69839783895572286297568834485025073557885364348071061715465477061873400359794989367423407683971299361817245213947182344090739843367076197016322541936552333837227080274674865687645877633828974738751695, 57407444219199061498801298672323590163238592201574572482836619025676869537007609315386800852204337587805249250896651467970585450518157701115893749407382500580682168292929753154872678880962261809848942, 41936641877539676652531567723249367917108638987131879521606243741814402095343724540844607212999297064816230828621064026263296602805535970091696570138772210094329631491849238240260893562065879302446837]

>>> from sage.all import *
>>> root_list = Zp.zeta(Integer(5),all=True); root_list
[80199770388563324100334548626345240081294273289109866436996652525328268652922508892946068538641538316054373187019168781211876036849531337824832226216684677717580165592175377569174402189281574130719978, 69839783895572286297568834485025073557885364348071061715465477061873400359794989367423407683971299361817245213947182344090739843367076197016322541936552333837227080274674865687645877633828974738751695, 57407444219199061498801298672323590163238592201574572482836619025676869537007609315386800852204337587805249250896651467970585450518157701115893749407382500580682168292929753154872678880962261809848942, 41936641877539676652531567723249367917108638987131879521606243741814402095343724540844607212999297064816230828621064026263296602805535970091696570138772210094329631491849238240260893562065879302446837]

root_list = Zp.zeta(5,all=True); root_list

Are these really fifth roots of unity?

sage: [root^5 for root in root_list]
[1, 1, 1, 1]

>>> from sage.all import *
>>> [root**Integer(5) for root in root_list]
[1, 1, 1, 1]

[root^5 for root in root_list]

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

More basic functionality

Similarly to the situation in linear algebra, there is much 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 equations,

  • 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 attempting to prove them.

sage: p = 13
sage: primitive_root(p); two_squares(p); is_prime(p)
2
(2, 3)
True

>>> from sage.all import *
>>> p = Integer(13)
>>> primitive_root(p); two_squares(p); is_prime(p)
2
(2, 3)
True

p = 13
primitive_root(p); two_squares(p); is_prime(p)

This makes it easy to construct elementary cryptographic examples as well. Here is a standard example of a Diffie-Hellman key exchange, for instance. If we didn’t do the second line, exponentiation would be impractical.

sage: p=random_prime(10^20,10^30) # a random prime between these numbers
sage: q=mod(primitive_root(p),p) # makes the primitive root a number modulo p, not an integer
sage: n=randint(1,p) # Alice's random number
sage: m=randint(1,p) # Bob's random number
sage: x=q^n; y=q^m
sage: x; y; x^m; y^n
66786436189350477660
77232558812003408270
45432410008036883324
45432410008036883324

>>> from sage.all import *
>>> p=random_prime(Integer(10)**Integer(20),Integer(10)**Integer(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(Integer(1),p) # Alice's random number
>>> m=randint(Integer(1),p) # Bob's random number
>>> x=q**n; y=q**m
>>> x; y; x**m; y**n
66786436189350477660
77232558812003408270
45432410008036883324
45432410008036883324

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
x; y; x^m; y^n

The final line of the cell first requests Alice and Bob’s (possibly) public information, and then verifies that the private keys they get are the same.

It is hard to resist including just one interact. How many theorems do you see here?

sage: @interact
sage: 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)

>>> from sage.all import *
>>> @interact
>>> def power_table_plot(p=(Integer(7),prime_range(Integer(50)))):
...     P=matrix_plot(matrix(p-Integer(1),[mod(a,p)**b for a in range(Integer(1),p) for b in srange(p)]),cmap='jet')
...     show(P)

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

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

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

sage: prime_pi(100); plot(prime_pi,1,100)
25
Graphics object consisting of 1 graphics primitive

>>> from sage.all import *
>>> prime_pi(Integer(100)); plot(prime_pi,Integer(1),Integer(100))
25
Graphics object consisting of 1 graphics primitive

prime_pi(100); plot(prime_pi,1,100)

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

sage: var('x')
x
sage: plot(prime_pi,2,10^6,thickness=2)+plot(Li,2,10^6,color='red')+plot(x/ln(x),2,10^6,color='green')
Graphics object consisting of 3 graphics primitives

>>> from sage.all import *
>>> var('x')
x
>>> plot(prime_pi,Integer(2),Integer(10)**Integer(6),thickness=Integer(2))+plot(Li,Integer(2),Integer(10)**Integer(6),color='red')+plot(x/ln(x),Integer(2),Integer(10)**Integer(6),color='green')
Graphics object consisting of 3 graphics primitives

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

Advanced Number Theory

For those who are interested, more advanced number-theoretic objects are easy to come by; we end with a brief sampler of these.

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.

sage: K.<a> = NumberField(x^2+14); K
Number Field in a with defining polynomial x^2 + 14

>>> from sage.all import *
>>> K = NumberField(x**Integer(2)+Integer(14), names=('a',)); (a,) = K._first_ngens(1); K
Number Field in a with defining polynomial x^2 + 14

K.<a> = NumberField(x^2+14); K

sage: K.discriminant(); K.class_group().order(); K.class_group().is_cyclic()
-56
4
True

>>> from sage.all import *
>>> K.discriminant(); K.class_group().order(); K.class_group().is_cyclic()
-56
4
True

K.discriminant(); K.class_group().order(); K.class_group().is_cyclic()

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

sage: complex_plot(zeta, (-30,30), (-30,30))
Graphics object consisting of 1 graphics primitive

>>> from sage.all import *
>>> complex_plot(zeta, (-Integer(30),Integer(30)), (-Integer(30),Integer(30)))
Graphics object consisting of 1 graphics primitive

complex_plot(zeta, (-30,30), (-30,30))

Cryptography

Sage supports various more advanced cryptographic procedures as well as some basic pedagogical ones natively. This example is adapted from the documentation.

sage: from sage.crypto.block_cipher.sdes import SimplifiedDES
sage: sdes = SimplifiedDES(); sdes
Simplified DES block cipher with 10-bit keys

>>> from sage.all import *
>>> from sage.crypto.block_cipher.sdes import SimplifiedDES
>>> sdes = SimplifiedDES(); sdes
Simplified DES block cipher with 10-bit keys

from sage.crypto.block_cipher.sdes import SimplifiedDES
sdes = SimplifiedDES(); sdes

sage: bin = BinaryStrings()
sage: P = [0,1,0,0,1,1,0,1] # our message
sage: K = sdes.random_key() # generate a random key
sage: C = sdes.encrypt(P, K) # encrypt our message
sage: plaintxt = sdes.decrypt(C, K) # decrypt it
sage: plaintxt # print it
[0, 1, 0, 0, 1, 1, 0, 1]

>>> from sage.all import *
>>> bin = BinaryStrings()
>>> P = [Integer(0),Integer(1),Integer(0),Integer(0),Integer(1),Integer(1),Integer(0),Integer(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]

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

See also the cryptography example in the discrete math quickstart.