4317 days ago by MathFest

Sage Mathfest I

Introduction to Sage

Neil Calkin and Dan Warner



Mission Statement

Create a viable free open source alternative to Magma, Maple, Mathematica, and Matlab.

History of Sage

  • William Stein started Sage at Harvard in January 2005.

  • Sage-1.0 released February 2006 at Sage Days 1 (UC San Diego).

  • Nearly 30 Sage Days Workshops (!)  at UCLA, UW, Cambridge, Bristol, Austin, France, San Diego, Seattle, MSRI, ..., Barcelona, ...

  • Sage won first prize in the Trophees du Libre (November 2007)

  • Funding from Microsoft, Univ of Washington, UC San Diego, NSF, DoD, Google, Sun, private donations, etc.


  • Hundreds of other people subsequently got involved in Sage's development, and the scope of the project has widened to cover all mathematical computation. There is interest from all areas of mathematics, physics, engineering, etc.  (Developer mailing list has 1184 subscribers and about 50 messages/day... and sage-support has 1718 subscribers.)


What is Sage?

  • Sage = Python + Math
  • A unified self-contained distribution of open source free mathematical software.
  • Nearly a half million lines of new Python (and Cython) code that implements new capabilities and algorithms.
  • A "cloud" application like GMail or Google Docs:  (currently about 30,000 users);  Of course Sage also runs on your desktop and supercomputer.


Sage is about
building the car instead of reinventing the wheel

  1. Sage uses Python,  a mainstream programming language, instead of inventing a custom mathematics language
  2. Use straightforward method to link programs together -- C library and pseudotty's, instead of XML servers/OpenMath.  We implement all conversion routines, instead of expecting upstream to do it: we make them communicate with Sage, whether they want to or not.
  3. Give copious credit to contributors and be very developer friendly (easily build from source).
  4. Reuse, improve, and contribute to existing libraries and projects (e.g., Scipy, Numpy, R, ATLAS, CVXopt, GSL), instead of starting over and competing with them.
  5. Make the GUI using a web browser: the world of java and javascript is available and Sage integrates with the web.



Key take-away points: 

  • FOSS: 100% free open source software: good for clusters, sharing, research, collaboration
  • Python: very mainstream, scientific-computing friendly programming language
  • Cython: optimized Python to C/C++ compiler we develop with support for C/C++ datatypes
  • Interfaces: to control Matlab, Octave, Mathematica, etc., from Sage
  • Self contained: can have many copies of Sage at once, so no "dependency hell"
  • Peer reviewed: get your code refereed
  • Web based: notebook interface

All aboard!

We want to pause and make sure that everyone can access Sage.

  1. Access the Sage Server at
  2. Install Sage on your machine, if not already downloaded get a flash drive with the correct version for your platform.



And now a short demo...

# This is a simple introduction to Sage. # # The Sage notebook uses blocks of text, which are essentially Python statements that get executed. # Lines that start with the pound sign are comments and do not get echoed when the block is evaluated. # # A block is evaluated by clicking on the evaluate link or by typing shift-enter. # # Pressing tab will give auto-complete options, or function descriptions after a parens ( # A '?' after a command gives detailed information, e.g. 'plot?' tells you about plot # Sage can do simple calculations 3+2*6^2 
# and more elaborate calculations factorial(100) 
# or 1000 digits of the square root of pi sqrt(pi).n(digits=1000) 
# Sage can factor integers factor(798) 
# and expand the previous result expand(_) 
# The same operations can be applied to algebraic expressions factor(x^5-1) 
p = expand(_) p 
# Algebraic expressions can be plotted P1 = plot( p , -0.5, 0.5, rgbcolor=(1,0,0) ) P2 = plot( p^3, -0.5, 0.5, rgbcolor=(0,1,0) ) show(P1+P2) 
# Sage handles rational functions with ease rf = x / (x^2 - 1) rf 
P3 = plot( rf , -0.5, 0.5, rgbcolor=(1,0,0) ) show(P3) 
# and partial fraction expansions pf = rf.partial_fraction(x) pf 
# Sage can display results in the more familiar mathematical format show(pf) 
# Sage uses LaTeX for the mathematical formatting latex(pf) 
# Sage can solve algebraic equations eq1 = x^2 - 2 == x show(eq1) solve([eq1],x) 
# Cubic polynomials are just as easy eq2 = x^3 + 2*x^2 - 5*x - 6 == 0 show(eq2) solve([eq2],x) 
# Sage can even find the ugly answers eq3 = x^3 - 2 == x show(eq3) solve([eq3],x) 

At this point we need to pause and look a little more carefully at how Sage works with Python.

First, Sage makes extensive use of classes

# What is pi? type(pi) 

# How do we tell Sage to compute pi to 4493 digits if we don't remember the command? # The tab key after "." of "(" generates a help display which can be used to # autocomplete the statement. pi. 
# Sage can solve equations for one variable in terms of others. # var('x b c') show(solve([x^2 + b*x + c == 0],x)) 

This example illustrated the fact that Python variables by default are "programming variables" and must be initialized. 

To create a symbolic variable we must declare the variable with the var command. 

The variable x is a special case for convenience.

# Sage can even find the ugly answers eq3 = x^3 - 2 == x show(eq3) sol1 = solve([eq3],x) for k in range(3): show(sol1[k]) 

This example illustrated the for loop in Python. 

First, the for statement walking through a list of objects and it ends with the colon. 

Second, the block inside the loop is indented.  There are no begin ... end constructs in Python. 

Good programming style strongly suggests that you indent subordinate blocks of code. 

Python insists that you indent and rewards you with shorter programs.

# Sage can solve linear systems of algebraic equations # Note the introduction of the symbolic variable y var('y') eq4 = x-1 == y+1 eq5 = x+1 == 2*(y-1) eq4;eq5 

This example underscored the fact that the single "=" is assignment and the double "==" is equality.

It also showed that multiple statements can be on the same line when separated by ";"

sol2 = solve([eq4,eq5],x,y) sol2 

In this example, the solution is returned in a list.  Lists are the principal data structure workhorse.

# Sage can also return the solutions in a dictionary # This provides access by name sol2 = solve([eq4,eq5],x,y,solution_dict=True) sd = sol2[0] print sd (sd[x].n(digits=2),sd[y].n(digits=3)) 

In this example we used the dictionary, which is the most sophisticated built-in data structure in Python.


Python Data Structures

The key data structures are: Lists, Tuples, and Dictionaries

Creating Lists I: [Square brackets]

L = [3, Permutation([5,1,4,2,3]), 17, 17, 3, 51] L 

Exercise: Create the list [63, 12, -10, 'a', 12], assign it to the variable L, and print the list.


Exercise: Create the empty list. (You will often need to do this.)


Creating Lists II: range

The range function provides an easy way to construct a list of integers. Here is the documentation of the range function.


range([start,] stop[, step]) -> list of integers

Return a list containing an arithmetic progression of integers.
range(i, j) returns [i, i+1, i+2, ..., j-1]; start (!) defaults to 0.
When step is given, it specifies the increment (or decrement).
For example, range(4) returns [0, 1, 2, 3].  The end point is omitted!
These are exactly the valid indices for a list of 4 elements.

Exercise: Use range to construct the list $[1,2,\ldots,50]$.


Exercise: Use range to construct the list of even numbers between 1 and 100 (including 100).


Exercise: The step argument for the range command can be negative. Use range to construct the list $[10, 7, 4, 1, -2]$.


Creating Lists III: list comprehensions

List comprehensions provide a concise way to create lists from other lists (or other data types). This technique is very similar to the mathematical notation for describing a set.

Example. We already know how to create the list $[1, 2, \ldots, 16]$:


Using a list comprehension, we can create the list of even numbers as follows:

E = [x for x in range(1,101) if x%2==0] E 

Using a list comprehension, we can also create the list $[1^2, 2^2, 3^2, ..., 16^2]$ as follows:

[i^2 for i in range(1,17)]
[i^2 for i in range(1,17)] 
sum([i^2 for i in range(1,17)]) 


The sum of the squares of the first ten natural numbers is,

1^(2) + 2^(2) + ... + 10^(2) = 385

The square of the sum of the first ten natural numbers is,

(1 + 2 + ... + 10)^(2) = 55^(2) = 3025

Hence the difference between the sum of the squares of the first ten natural numbers and the square of the sum is 3025 − 385 = 2640.

Find the difference between the sum of the squares of the first one hundred natural numbers and the square of the sum.

s1 = sum( 
s2 = ( 
s2 - s1 

Filtering lists with a list comprehension

A list can be filtered using a list comprehension.

Example: To create a list of the squares of the prime numbers between 1 and 100, we use a list comprehension as follows.

[p^2 for p in [1,2,..,100] if is_prime(p)] 
[p^2 for p in range(1,101) if is_prime(p)] 

Exercise: Use a list comprehension to list all the natural numbers below 20 that are multiples of 3 or 5.


  1. To get the remainder of 7 divided by 3 use 7%3.
  2. To test for equality use two equal signs (==); for example, 3 == 7.

Nested list comprehensions

List comprehensions can be nested!


[(x,y) for x in range(5) for y in range(3)] 
[[i^j for j in range(1,4)] for i in range(6)] 
matrix([[i^j for j in range(1,4)] for i in range(6)]) 


Pythagorean triple is a triple $(x,y,z)$ of positive integers satisfying $x^2+y^2=z^2$. 

The Pythagorean triples whose components are at most $10$ are :  $[(3, 4, 5), (4, 3, 5), (6, 8, 10), (8, 6, 10)]$.

Using a filtered list comprehension, construct the list of Pythagorean triples whose components are at most $50$.


Accessing individual elements of lists

To access an element of the list, use the syntax L[i], where i is the index of the item.


  1. Construct the list L = [1,2,3,4,3,5,6]. What is L[3]?
  2. What is L[1]?
  3. What is the index of the first element of L?
  4. What is L[-1]? What is L[-2]?
  5. What is L.index(2)? What is L.index(3)?


L = [1,2,..,6];L 

Modifying lists: changing an element in a list

To change the item in position i of a list L:

L = ["a", 4, 1, 8] L 
L[2] = 0 L 
L[1:3]=[1,2,-1] L 

Modifying lists: append and extend

To append an object to a list:

L = ["a", 4, 1, 8] L 
L.append(17) L 

To extend a list by another list:

L1 = [1,2,3] L2 = [7,8,9,0] 
print L1 print L2 
L1.extend(L2) L1 

Modifying lists: reverse, sort, ...

L = [4,2,5,1,3] L 
L.reverse() L 
L.sort() L 
L = [3,1,6,4] 

Concatenating Lists

To concatenate two lists, add them (+). This is not a commutative operation....

L1 = [1,2,3] L2 = [7,8,9,0] print L1 print L2 
L1 + L2 

Slicing Lists

You can slice a list using the syntax L[start : stop : step]. This will return a sublist of L.

Exercise: Below are some examples of slicing lists. Try to guess what the output will be before evaluating the cell.

L = range(20) L 
L[3:5] = [-2,'a'];L 

Exercise: The following function combines a loop with the some of the list operations above. What does the function do?

def f(number_of_iterations): L = [1] for n in range(2, number_of_iterations): L = [sum(L[:i]) for i in range(n-1, -1, -1)] return numerical_approx(2*L[0]*len(L)/sum(L), digits=50) 


A tuple is an immutable list. That is, it cannot be changed once it is created.

The syntax for creating a tuple is to the use parentheses. 

t = (3, 5, [3,1], (17,[2,3],17), 4) t 

We can create a tuple from a list, or vice-versa.


Tuples behave like lists in many respects:

Operation Syntax for lists Syntax for tuples
Accessing a letter
list[3] tuple[3]
Concatenation list1 + list2 tuple1 + tuple2
Slicing list[3:17:2] tuple[3:17:2]
A reversed copy
list[::-1] tuple[::-1]
Length len(list) len(tuple)

Trying to modify a tuple will fail.

t = (5, 'a', 6/5) t 
t[1] = 'b' 


"Tuple-comprehension" does not exist.  The syntax produces something called a generator.  A generator allows you to process a sequence of items one at a time.  Each item is created when it is needed, and then forgotten.  This can be very efficient if we only need to use each item once.

(i^2 for i in range(5)) 
g = (i^2 for i in range(5)) 
[x for x in g] 

g is now empty.

[x for x in g] 

A nice 'pythonic' trick is to use generators as the argument to functions.  We do not need double parentheses for this.

sum( i^2 for i in srange(100001) ) 


A dictionary is another built-in data type. Unlike lists, which are indexed by a range of numbers, dictionaries are "indexed" by keys, which can be any immutable object. Strings and numbers can always be keys (because they are immutable). Dictionaries are sometimes called "associative arrays" in other programming languages.

There are several ways to define dictionaries. One method is to use braces, {}, with comma-separated entries given in the form key:value.

d = {3:17, "key":[4,1,5,2,3], (3,1,2):"goo", 3/2 : 17} d 

Dictionaries behave as lists and tuples for several important operations.

Operation Syntax for lists Syntax for dictionaries
Accessing elements
L[3] D["key"]
Length len(L) len(D)
Modifying L[3] = 17 D["key"] = 17
Deleting items del L[3] del D["key"]

A dictionary can have the same value multiple times, but each key can only appear once must be immutable.

d = {3 : 14, 4: 14} d 
d = {3: 13, 3:14} d 
d = {[1,2,3] : 12} d 
d = {(1,2,3) : 12} d 

Another way to add items to a dictionary is with the update() method:

d = {} d 
d.update( {10 : 'newvalue', 20: 'newervalue', 3: 14, 'a':[1,2,3], 15:[]} ) d 
d[15].append(7) d 

We can iterate through the keys, or values, or both, of a dictionary.

d = {10 : 'newvalue', 20: 'newervalue', 3: 14, 'a':[1,2,3]} 
[key for key in d] 
[key for key in d.iterkeys()] 
[value for value in d.itervalues()] 
[(key, value) for key, value in d.iteritems()] 

Exercise: Consider the following directed graph. 

Create a dictionary whose keys are the vertices of the above directed graph, and whose values are the lists of the vertices that it points to. For instance, the vertex 1 points to the vertices 2 and 3, so the dictionary will look like:

d = { ..., 1:[2,3], ... }

Then try

g = DiGraph(d)



Using sage types: The srange command

Example:  Construct a $3 \times 3$ matrix whose $(i,j)$ entry is the rational number $\frac{i}{j}$.

matrix([[ i/j for j in range(1,4)] for i in range(1,4)]) 
matrix([[ i/j for j in srange(1,4)] for i in srange(1,4)]) 

Modifying lists has consequences!

Try to predict the results of the following commands.

a = [1, 2, 3] L = [a, a, a] L 
a.append(4) L 

Now try these:

a = [1, 2, 3] L = [a, a, a] L 
a = [1, 2, 3, 5] L 
L[0].append(4) L 

You can use the command deepcopy to avoid these issues.

a = [1,2,3] L = [deepcopy(a), deepcopy(a)] L 
a.append(4) L 

The same issues occur with dictionaries.

d = {1:'a', 2:'b', 3:'c'} dd = d d.update( { 4:'d' } ) dd 

Some Sums

Let's now explore binomial expansions, Pascal's triangle, and more interesting sequences of sums.

var('x,y') p = [1] n = 7 for k in range(1,n): p.append(expand((x+y)*p[k-1])) print p[k] 
# That leads directly to Pascal's triangle n = 13 Pt = matrix(ZZ,n) Pt[0,0] = 1 for i in range(1,n): Pt[i,0] = 1 for j in range(i): Pt[i,j+1] = Pt[i-1,j] + Pt[i-1,j+1] print Pt 
pc = Pt.columns() # print pc s0 = copy(pc[0]) for k in range(1,len(s0)): s0[k] += s0[k-1] print s0 s1 = copy(pc[1]) for k in range(1,len(s1)): s1[k] += s1[k-1] print s1 s2 = copy(pc[2]) for k in range(1,len(s2)): s2[k] += s2[k-1] print s2 s3 = copy(pc[3]) for k in range(1,len(s3)): s3[k] += s3[k-1] print s3 

Observation 1:

If we sum the first $N$ elements in column $k$ that sum will be the value of the element in the next row and column, namely, row $N+1$ and column $k+1$.

$\sum\limits_{n=k}^N P_{n,k} = P_{N+1,k+1}$

Since $P_{n,k} = \left( {\begin{array}{*{20}c} n \\ k \\ \end{array}} \right)$ this is the same as the statement

$\sum\limits_{n=k}^N \left( {\begin{array}{*{20}c} n \\ k \\ \end{array}} \right) = \left( {\begin{array}{*{20}c} N+1 \\ k+1 \\ \end{array}} \right)$

or, if we use the factorial definition of $n$ choose $k$

$\sum\limits_{n=k}^N \frac{n!}{(n-k)! k!} = \frac{(N+1)!}{(N-k)!(k+1)!}$

Actually, this observation holds for all cases, and the Theorem is easy to prove by mathematical induction. 


Observation 2:

Look at each row of Pascal's Triangle.  In all the cases in the Table, if $n$ is a prime (2, 3, 5, 7, 11), then $n$ divides all the entries in the row between the first and last, which are ones.  On the other hand, when $n$ is not a prime (4, 6, 8, 9, 10), then $n$ does not divide every entry in the row.  In particular, 4 doesn't divide 6; 6 doesn't divide 15 or 20; 8 doesn't divide 28; 9 doesn't divide 84; 10 doesn't divide 45; and 12 doesn't divide 66.

This observation also holds for all cases, and the result along with Eisenstein's criterion can be used to show that polynomials of the form

$x^{p-1} + x^{p-2} + \cdots + x + 1$

cannot be factored over the rational numbers when $p$ is a prime.

# Now let's explore some partial row sums nc = 6 Pr = matrix(ZZ,n,nc) for j in range(nc): for i in range(n): if j == 0: Pr[i,j] = Pt[i,j] else: Pr[i,j] = Pt[i,j] + Pr[i,j-1] print Pr 
pc = Pr.columns() 
# A000124 numbers -- Check with the Online Encyclopedia of Integer Sequences print pc[2] 
# ?? Numbers g = Graphics() g += circle((2,2),2) g += line([(0,0),(4,4)]) g += line([(0,2),(4,2)]) g += line([(0,4),(4,0)]) # g += line([(0,3.5),(3.5,0)], rgbcolor=(1,0,0)) show(g,aspect_ratio=1) 
# A000125 numbers -- Check with the Online Encyclopedia of Integer Sequences print pc[3] 
# A000127 numbers -- Check with the Online Encyclopedia of Integer Sequences print pc[4] 
np = 3 g = Graphics() g += circle((2,2),2) ck = cos(2.0*pi/np) sk = sin(2.0*pi/np) xk = 0 yk = 2 points = [(xk+2,yk+2)] for k in range(1,np): xt = ck*xk - sk*yk yk = sk*xk + ck*yk xk = xt points.append((xk+2,yk+2)) g += point(points,pointsize=20) for i in range(np): for j in range(i,np): g += line([points[i],points[j]]) show(g,aspect_ratio=1) 
# A006261 numbers -- Check with the Online Encyclopedia of Integer Sequences print pc[5] 
# A000045 numbers diag_sums = [] for i in range(n): s = 0 for j in range(i+1): s += Pt[i-j,j] diag_sums.append(s) print diag_sums 
F = transpose(Matrix(diag_sums)) print F 
R = Matrix(ZZ,n) for i in range(n): R[n-1-i,i] = 1 
F2 = Pt*R*F print F2 
F3 = Pt*R*F2 print F3 
print (F[n-1]/F[n-2]).n(digits=30) print (F2[n-1]/F2[n-2]).n(digits=30) print (F3[n-1]/F3[n-2]).n(digits=30)