CoCalc Public FilesSD109 - Tutorial.ipynbOpen with one click!
Author: Fredrik Strömberg
Views : 200
Description: Tutorial made for the Virtual Global Sage Days 109

Introduction to Python and SageMath

We will introduce some of the fundamental concepts in Python and SageMath.

Learning outcomes

  1. Introduction to SageMath
  2. Learn how to run Sage command line or notebook
  3. Basic expressions and variables in Python and SageMath
  4. Introduce some basic packages for plotting etc.
  5. For the full lessons see

What is SageMath?

SageMath was originally developed by William Stein and first released in 2005. It was initially called SAGE, then Sage, and now SageMath, or Sage for short. The intention behind SageMath is to provide a free, open source alternative to Magma, Maple, Mathematica, MATLAB etc. with an easy to use interface (both for end-users and developers) built on Python (Cython) and which can integrate with other free or commercial packages (installed separately).

The default installation guide (wiki) or official page: installation contains many (over 150) different free packages, including:

  • Pari/GP
  • GAP
  • Singular
  • R
  • Numpy
  • matplotlib
  • ...

For a full list see:

It can be integrated with existing commercial packages like

  • Magma
  • Mathematica
  • Maple

Full documentation is included in the source of SageMath and is available online, together with a multitude of tutorials etc.
See for instance

Optional packages

Some packages that are not shipped with SageMath are available to install afterwards.

  • In the terminal (on Windows, use the "SageMath Shell"):
    • sage -i <package_name> - to install optional Sage package
    • sage --pip install <py_package_name> - to install extra Python package using pip
  • In a SageMath session (on Windows, use "SageMath" or "SageMath Notebook"):
    • sage: installed_packages() - to list all installed packages
    • sage: optional_packages() - to list installed and available optional packages
In [ ]:
For example, - install additional GAP packages by running `sage -i gap_packages` - install JupyterLab by running `sage --pip install jupyterlab` - install Pandas by running `sage --pip install pandas`

Running SageMath

On your own machine

For the Sage REPL (read-eval-print loop) or command-line interactive use:

  • sage to run Sage
  • sage --help to see all options, e.g. --python, --pip, etc.

For the Jupyter Notebook interface, run one of:

  • sage -n jupyter --notebook-dir=<dir>
  • sage --notebook=jupyter --notebook-dir=<dir>
  • sage --jupyter notebook --notebook-dir=<dir>

For the JupyterLab interface (if installed -- see above), run one of:

  • sage -n jupyterlab --notebook-dir=<dir>
  • sage --notebook=jupyterlab --notebook-dir=<dir>
  • sage --jupyter lab --notebook-dir=<dir>


Getting help with SageMath

  • <function_name>? - to read documentation (docstring)
  • <function_name>?? - to read the source (if available)
  • <object_name>.<tab> - use tab completion to see available properties
  • Help button in the toolbar links to the online documentation

Example: Try using these features for yourself below. Run a code cell by first selecting it and then either

  • Press the Run button in the toolbar above
  • Use the keyboard shortcut Ctrl+Enter or Shift+Enter
In [ ]:
In [ ]:
In [ ]:

Programming in SageMath

Check that your kernel is set to SageMath (e.g. SageMath 9.0) in the top right of the toolbar.
You can change your kernel using the Kernel button on the toolbar to set the default to Python (2.7 or 3.x depending on your version of SageMath).
Since version 9.0 released in January 2020, SageMath is using Python 3.

Since the main functionality of SageMath is written in Python it means that most of the programming can be done using Python. There are however some subtle differences.

Variable types in Python

  • Numerical types: int, float, complex
  • Strings: str
  • Lists: list
  • Tuple: tuple
  • Dictionary: dict
  • Set: set

Example: Run the cell below to explore basic data types in Python

In [ ]:
%%python3 # We run Jupyter cells in SageMath mode by default but can run individual cells in Python using Jupyter magics # As Python 2 is still the default for older versions of SageMath it is better to be explicit about Python version. x = 2 print(x, type(x)) # int # Combine multiple statements with ';' -- avoid this in practice since it reduces readability y = 2.0; print(y, type(y)) # float r = 3/4; print(r, type(r)) # float z = 1 + 1j; print(z, type(z)) # complex s = "hello"; print(s, type(s)) # string l = [1, 2, 2]; print(l, type(l)) # list t = (1, 2, 2); print(t, type(t)) # tuple d = {1: 'a', 'b': 2}; print(d, type(d)) # dictionary A = set([1, 2, 3, 2]); print(A, type(A)) # set

Note that in Python 2, dividing two int variables uses "integer division with floor floor function"

In [ ]:
%%python2 r = 3/4; print(r, type(r)) r = -3/4; print(r, type(r))

Variable types in SageMath

In SageMath defaults for many numeric types are overwritten. For example:

  • Integers are by default cast to SageMath's Integer type, which has different behaviour from Python's int
  • Instead of float, SageMath uses RealLiteral and Rational (where appropraite)
In [ ]:
x = 2 # Integer y = 2.0 # Real/floating point numbers with `RealLiteral` r = 3/4 # Rational z_num = 1 + 1j # Numerical complex numbers with `ComplexNumber` z_sym = 1 + 1*I # Symbolic complex numbers with `Expression` # Reduce writing and improve readability by iterating over a list using a `for` loop for variable in [x, y, r, z_num, z_sym]: print(variable, type(variable))


Sets are useful for making lists with unique elements, are order ambivalent, and can be defined using either braces {} or the set() constructor function.

In [ ]:
X = set([1, 2, 3, 4, 4, 4, 5]) Y = {5, 3, 2, 4, 2} print(X) print(Y) print(list(X))

The set type features many useful methods such as unions and intersections. See the documentation here.

In [ ]:
A = set([1, 2, 3, 4]) B = set(['a', 'b', 'c', 1, 2, 8]) A.intersection(B)

Sage has some additional useful builtin-functions:

In [ ]:
In [ ]:

Precision of variables in SageMath

SageMath stores variables to higher levels of precision than Python

  • The RR keyword denotes the Real Field with 53 bits of precision (double precision)
  • The CC keyword denotes the Complex Field with 53 bits of precision
  • We can construct fields with higher precision using RealField(prec) and ComplexField(prec)
In [ ]:
for field in [RR, RealField(106), CC, ComplexField(106)]: print(field)

We can contruct high precision complex numbers as field elements as follows:

In [ ]:
my_complex_field = ComplexField(106) z1 = my_complex_field(1, 10) z2 = my_complex_field(5, 5) for num in [z1, z2, z1 + z2]: print(num)

There are also builtin constants that can be given to arbitrary precision

In [ ]:
print(RR.pi()) print(RealField(106).pi()) high_prec_one = RealField(1000)(1) print(high_prec_one) print(exp(high_prec_one)) # the builtin `exp` function print(high_prec_one.exp()) # the `exp` method of the `RealNumber` class

Objects, Parents, Elements and categories

In Python and even more so in SageMath everything is an object and has properties and methods.
In SageMath most (every?) object is either an Element or a Parent

In [ ]:
F = RealField(53) # Parent x = F(2) # Element print(type(F)) print(type(x)) F.is_parent_of(x) # True print(x in F) # True
In [ ]:
x.parent() # Returns a copy of F
In [ ]:
x.category() # The category of x, a category of elements
In [ ]:
F.category() # The category of F
In [ ]:
In SageMath there is something called "coercion": natural operations are made using common parents.
In [ ]:
x=RR(1); y=CC(5) type(x+y)


String formatting

In previous versions of Python there were more types of strings but in current Python there are basically only strings str and bytestrings bytes (which we do not discuss here)

We can control how strings are formatted when printing by using either the format method or by using formatted strings, which are strings with an f character before the first quotation mark. Please see the examples below, and follow this link for a reference on string formatting.

In [ ]:
x = RR(2) # The real number 2, up to 53 bits of precision print("x = {0}".format(x)) # format method print(f"x = {x}") # f-strings print(f"exp(x) = {x.exp()}") # Can call functions in f-strings print(f"π is about {float(RR.pi()):0.5f}") # As floating point number with 5 decimals

Exercise 1. Print eπe^{\pi} as a floating point number to 225 decimal places.

Python solution: Use the decimal library (documentation)

In [ ]:
%%python3 import decimal, math import numpy as np decimal.getcontext().prec = 225 print(f"{decimal.Decimal(math.pi)}\n{decimal.Decimal(np.pi)}")

The cell above shows that there is not enough precision in the stored values of pi in the math and numpy libraries for our application, so we must use an algorithm to calculate it to a higher precision. We use the Gauss-Legendre algorithm for its fast convergence rate.

In [ ]:
%%python3 from decimal import Decimal, getcontext import numpy as np getcontext().prec = 227 num_iter = 8 one = Decimal(1) two = Decimal(2) four = Decimal(4) def gauss_legendre_algorithm(num_iter): r"""Implement the Gauss-Legendre algorithm for the specified number of iterations Return an `np.array` of approximations to pi of increasing accuracy""" a, b, t, = [one], [two.sqrt()/two], [one/four] approx = list() for i in range(num_iter): a.append((a[i] + b[i]) / two) b.append((a[i] * b[i]).sqrt()) t.append(t[i] - getcontext().power(two, Decimal(i)) * (a[i] - a[i+1]) * (a[i] - a[i+1])) approx.append((a[i+1] + b[i+1]) * (a[i+1] + b[i+1]) / (four * t[i+1])) return np.array(approx) approx = gauss_legendre_algorithm(num_iter) print(approx == approx[-1]) # Stabilised to 225 decimal places precision after 7 iterations precise_pi = approx[-1] result = precise_pi.exp() print(result) print(f"Number of decimal places = {len(str(result).split('.')[-1])}") # Count digits after decimal point

SageMath solution:

  • Note that this can be calculated using many fewer lines of code since there are builtin SageMath algorithms running behind the scenes.
  • Note also that we must convert from bit (binary) precision to precision in decimal places by multiplying by log2(10)3.32\log_{2}(10) \approx 3.32
In [ ]:
required_bit_prec = 228 * RR(10).log2() # Calculate required binary precision from decimal exp_pi = RealField(required_bit_prec).pi().exp() print(exp_pi) print(f"Number of decimal places = {len(str(exp_pi).split('.')[-1])}") # Count digits after decimal point


  • Function names should be descriptive and use "snake_case" (this is a Python convention)
  • Variable names in Python should also be lower, snake_case
  • In SageMath we also use mathematical conventions even if they break this, e.g. Ei for the exponential integral function
  • For full document about coding conventions see:
In [ ]:
# A first function def is_zero_mod_3(x): return x % 3 == 0 # The % is the Python modulo operation: remainder in integer division
In [ ]:
In [ ]:
# Does this make sense? is_zero_mod_3(5.4)
In [ ]:
is_zero_mod_3(1+1j) # Raises type error

Sage has an alternative modulo reduction function

In [ ]:
Mod(5, 3)
In [ ]:
type(Mod(5, 3))
In [ ]:
# Add handling of input def is_zero_mod_3(x): if not isinstance(x, int): raise ValueError(f"Received input of type {type(x)}. This function needs an integer!") return x % 3 == 0
In [ ]:
# Better and more informative error is raised is_zero_mod_3(1 + 1j)
In [ ]:
In [ ]:
# However, this doesn't work as intended is_zero_mod_3(6)

The issue is that we checked for type int but in a SageMath environment the default type for integers is Integer. Need to check for both types!

In [ ]:
# Add handling of input def is_zero_mod_3(x): if not isinstance(x, (int, Integer)): # Multiple types should be given as a tuple raise ValueError(f"Received input of type {type(x)}. This function needs an integer!") return x % 3 == 0
In [ ]:


When writing our own functions we can help future users check the expected inputs in advance by writing a docstring (documentation string).

In [ ]:
In [ ]:
# Let's add a docstring def is_zero_mod_3(x): r""" Return True if input is congruent to 0 mod 3, otherwise return False INPUT: - ``x`` -- integer OUTPUT: A boolean describing whether the input is congruent to 0 mod 3 or not. """ if not isinstance(x, (Integer, int)): raise ValueError(f"Received input of type {type(x)}. This function needs an integer!") return x % 3 == 0
In [ ]:
# Now a potential user knows what input and output to expect. is_zero_mod_3?
In [ ]:

Variables in functions:

In [ ]:
# Scalar values get passed by value into functions, i.e. what happens in the function stays in the function: x = 1 def add_one(x): x = x + 1 print(f"x = {x}") y = 1 print(y) add_one(y) print(y) print(x)
In [ ]:
# Dictionaries get passed by reference, i.e. what happens in the function also happens outside. d = {'a': 1} print(d['a']) # Get d['a'] if d has 'a' as a key, raise an error if not print(d.get('b', 'default')) # Get d['b'] if 'b' is a key, otherwise return 'default' (or None if no default value is given) def add_one(d): d['a'] = d.get('a') + 1 print(f"d = {d}") print(d) add_one(d) print(d)
In [ ]:
[1] + 1 # Addition of integers and lists in this way is not supported

Type hints

In Python 3 it is also possible to add type hints (not enforced by Python).

In [ ]:
# Let's add a docstring def is_zero_mod_3(x: int) -> bool: r""" Return True if input is congruent to 0 mod 3, otherwise return False INPUT: - ``x`` -- integer OUTPUT: A boolean describing whether the input is congruent to 0 mod 3 or not. """ if not isinstance(x, (Integer, int)): raise ValueError(f"Received input of type {type(x)}. This function needs an integer!") return x % 3 == 0

Exercise Write a function is_square_residue with the following specifications:

  1. Takes as input an integer xx and a positive integer nn
  2. Returns True if xx is congruent to a square modulo nn and otherwise False.
  3. Handles errors in input with sensible error messages.
  4. Includes a docstring which describes the function.
In [ ]:
def is_square_residue(x, n: int) -> bool: r""" Return `True` if $x$ is congruent to a squared integer modulo $n$, otherwise return `False` INPUTS: - ``x`` -- an integer - ``n`` -- a positive integer OUTPUT: A boolean describing whether there are any solutions $a$ to $x \equiv a^{2} \mod n$ for inputs `x` and `n` """ if not (isinstance(x, (Integer, int)) and isinstance(n, (Integer, int)) and n > 0): raise ValueError(f"Received inputs (x, n) = ({x}, {n}) of type ({type(x)}, {type(n)}). Inputs should be integers, with n positive!") return (x % n) in set(a**2 for a in range(n))
In [ ]:
n = 8 for x in range(n): print(x, is_square_residue(x, n))

Functions and objects

Functions are also objects and can be introspected via their properties. Class methods and imports etc. can be modified at runtime (monkey-patching)

In [ ]:
In [ ]:
# some property... is_zero_mod_3.__code__.co_varnames

Example of monkey patching

In [ ]:
%%python2 import six class TestA(object): def foo(self,x): return x if six.PY3: = lambda self,x : x+1 print(TestA().foo(1))
In [ ]:
%%python3 import six class TestA(object): def foo(self,x): return x if six.PY3: = lambda self,x : x+1 print(TestA().foo(1))

Python Control Structures

Standard (used in most programming languages):

  • Iteration with for and while loops
  • If-then-else statements

More Python specific:

  • Generator expressions
  • List comprehensions
In [ ]:
# Iterate over a range of integers using the `range` function. Note that range starts from 0. for i in range(5): print(i)

The output of the range function is an example of a generator expression. It does not actually allocate all elements until needed.

In [ ]:
range(5, 12) # Specify start and end points

We can cast a range to a list to see all of its elements. Note that range starts at the left endpoint and stops at the right endpoint without including it!

In [ ]:
list(range(5, 12))

Calling list(range(10^(10^10))) would run out of memory, but iterating over it is fine (although probably won't finish).
If evaluating the cell below, please call keyboard interrupt or select Kernel -> Interrupt from the toolbar above.

In [ ]:
for i in range(10^(10^10)): pass

Example: for loops and list comprehensions

We can evaluate ζ(2)\zeta(2), the zeta function at s=2, using a loop, and compare to the known value π6\frac{\pi}{6}.

In [ ]:
# Calculate the partial sum of the first 100 terms of `zeta(2)` result = 0 for n in range(1, 100): result += n**(-2) print(result)

We may also use a list comprehension, or directly use the generator expression to compute the partial sum

In [ ]:
sum([n**(-2) for n in range(1, 100)]) # List comprehension
In [ ]:
sum(n**(-2) for n in range(1, 100)) # Generator expression - most memory efficient

Plotting in Sagemath with the plot function

The plot() function in SageMath takes as its first argument a function or list of functions. In the case that the passed functions take a single argument, the next two arguments, start and end, can be used to define the range over which to evaluate and plot them.

In [ ]:
In [ ]:
# A first plot of zeta along the horizontal line Im(z) = 1, 2 < Re(z) < 10 p = plot(lambda x: zeta(CC(x, 1)).real(), 2, 10)
In [ ]:

Most objects have a latex property which "pretty-prints" them in a format suitable for inclusion in a LaTeX document.

In [ ]:

Symbolic expressions in SageMath

When starting Sage the name 'x' is defined as a symbolic variable and other variables have to be declared using the var('x, y, z, ...') command. Symbolic expressions can be treated in various ways:

  • differentiation
  • simplifications

When the notebook is started we have xx as a variable:

In [ ]:

We can define additional variables

In [ ]:
var('x, y, z')
In [ ]:


In [ ]:
g = 1 / (x^2 + y^2) print(f"{'g':20s} = {g}") print(f"{'dg/dx':20s} = {diff(g, x)}") print(f"{'d^2g/dxdy':20s} = {diff(g, x, y)}")
In [ ]:
g.differentiate(x, 2)


In [ ]:
f = x * y / (x^2 + y^2 ) z = diff(f, x, y) z
In [ ]:


In [ ]:
z1=z.substitute(x=2*y^2+1); z
In [ ]:

Exercise Determine the value of fu(0,1,1)\frac{\partial f}{\partial u}(0,1,1) for f(s,t,u)=s+tus2+t2+u2 f(s,t,u)=\frac{s+tu}{\sqrt{s^2+t^2+u^2}} Is it

  • (a) 2{\sqrt{2}}

  • (b) 12\frac{1}{\sqrt{2}}

  • (c) 122\frac{1}{2\sqrt{2}}

  • (d) 152\frac{1}{5\sqrt{2}}

  • (e) None of the above?

In [ ]:
var('s,t,u') f = (s+t*u)/sqrt(s**2+t**2+u**2) diff(f,u).substitute(s=0,t=1,u=1)

We can now combine everything and compute the Bernoulli numbers. Recall that they can be defined by the generating series xex1=m0Bmxmm! \frac{x}{e^x - 1} = \sum_{m \ge 0} \frac{B_m x^m}{m!}

In [ ]:
F = x / (e^x - 1); F # The generating function

Try to find the first Taylor coefficient:

In [ ]:
g = derivative(F, x, 1); g # another name for .diff()
In [ ]:
print(g.simplify_full()) # yes - there are other types of simplifications. g.substitute(x=0)
In [ ]:
We can't just divide 0 / 0. We need L'Hopital's rule!
In [ ]:
# Differentiate the numerator and denominator and divide again: g.numerator().diff(x) / g.denominator().diff(x)

Still of the form 0 /0. Need one more derivative!

In [ ]:
# The second parameter gives the number of times we differentiate p = g.numerator().diff(x, 2) / g.denominator().diff(x, 2) print(p) p = p.simplify_full() print(p) p.substitute(x=0)
In [ ]:

So the first Bernoulli number B1=12B_1=-\frac{1}{2}. This method is a bit cumbersome but fortunately there is a builtin command in Sage for Taylor expansions

In [ ]:
F.taylor(x, 0, 10)

We can convert this to a polynomial over Q\mathbb{Q}:

In [ ]:
p = F.taylor(x, 0, 10).polynomial(QQ) print(type(p)) print(p) print(latex(p))

For a polynomial we can add a big-Oh

In [ ]:
q = p.add_bigoh(12); q
In [ ]:
print(q.parent()) type(q)
In [ ]:
x = q.parent().gen() x
In [ ]:
q + (x + 1).add_bigoh(8)

We can get coefficients of certain terms in Taylor expansions

In [ ]:
F.taylor(x, 0, 10).coefficient(x^4)

We can now write a function that returns the j-th Bernoulli number

In [ ]:
def B(j): F = x / (e^x - 1) return F.taylor(x, 0, j).coefficient(x^j)*factorial(j) [B(j) for j in range(1, 10)]

We can also work with polynomials in many variables

In [ ]:
F=GF(3)['x, y'] F
In [ ]:
x,y=F.gens(); x+4*y
In [ ]:
f=x*y+2*x**2 g=x**2+y*x**3
In [ ]:

Linear algebra can be done using "native" Sage or numpy

In [ ]:
m=Matrix(ZZ,[[1,1],[2,3]]); m
In [ ]:
In [ ]:
m.eigenvalues() ## Root finding in the above

Algebraic numbers

What is the "?" ?

In [ ]:
print(m.eigenvalues()[0].parent()) print(type(m.eigenvalues()[0]))

Can be evaluated to arbitrary precision

In [ ]:

Can use Numpy

In [ ]:
In [ ]:
import numpy numpy.linalg.eigvals(m)

Can work over finite fields

In [ ]:
m=Matrix(GF(5),[[1,1],[2,3]]); m
In [ ]:
m.eigenvalues() # "exact"
In [ ]:
In [ ]:
Does it work with interval arithmetic?
In [ ]:
RIF=RealIntervalField() h=1e-10 x=(1-h,1+h) x
In [ ]:
m=Matrix(RealIntervalField(53),[[RIF(1-h,1+h),RIF(1-h,1+h)],[RIF(2-h,2+h),RIF(3-h,3+h)]]); m
In [ ]:
m.eigenvalues() # Not implemented!

Can use Newton - Raphson to find the roots:

In [ ]:
p=m.charpoly() p.newton_raphson(10,0.2)
In [ ]:
In [ ]:
In [ ]: