SD109 - Tutorial.ipynb

Kernel: SageMath (system-wide)

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 https://github.com/fredstro/sage-lesson-nt

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: https://www.sagemath.org/links-components.html

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 https://www.sagemath.org/tour.html

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

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>

Online

• SageMathCell: https://sagecell.sagemath.org/ - write SageMath code directly in the browser

• Cocalc: https://cocalc.com - full interactive notebook environment (compare with Google Colab)

• Binder: https://opendreamkit.org/try/ - temporary instance

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

gcd?

x=1

x.

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

%%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"

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

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

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

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.

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

Sage has some additional useful builtin-functions:

list(cartesian_product([A,B]))

list(powerset(B))

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)

for field in [RR, RealField(106), CC, ComplexField(106)]:
    print(field)

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

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

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

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

x.parent()  # Returns a copy of F

x.category()  # The category of x, a category of elements

F.category()  # The category of F

In SageMath there is something called "coercion": natural operations are made using common parents.

x=RR(1); y=CC(5)
type(x+y)

CC.coerce_map_from(RR)

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.

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^{\pi}$ as a floating point number to 225 decimal places.

Python solution: Use the decimal library (documentation)

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

%%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 $\log_{2}(10) \approx 3.32$

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

Functions

• 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: http://doc.sagemath.org/html/en/developer/coding_basics.html

# A first function
def is_zero_mod_3(x):
    return x % 3 == 0  # The % is the Python modulo operation: remainder in integer division

is_zero_mod_3(5)

# Does this make sense?
is_zero_mod_3(5.4)

is_zero_mod_3(1+1j)  # Raises type error

Sage has an alternative modulo reduction function

Mod(5, 3)

type(Mod(5, 3))

# 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

# Better and more informative error is raised
is_zero_mod_3(1 + 1j) 

is_zero_mod_3(int(6))

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

# 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

is_zero_mod_3(6)

Docstrings

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

is_zero_mod_3?

# 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

# Now a potential user knows what input and output to expect.
is_zero_mod_3?

is_zero_mod_3("test")

Variables in functions:

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

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

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

# 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 $x$ and a positive integer $n$

2. Returns True if $x$ is congruent to a square modulo $n$ and otherwise False.

3. Handles errors in input with sensible error messages.

4. Includes a docstring which describes the function.

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

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)

category(is_zero_mod_3)

# some property...
is_zero_mod_3.__code__.co_varnames

Example of monkey patching

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

print(TestA().foo(1))

%%python3
import six
class TestA(object):
    def foo(self,x):
        return x
if six.PY3:
    TestA.foo = 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

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

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!

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.

for i in range(10^(10^10)):
    pass

Example: for loops and list comprehensions

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

# 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

sum([n**(-2) for n in range(1, 100)])  # List comprehension

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.

plot?

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

p

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

latex(p)

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 $x$ as a variable:

x

We can define additional variables

var('x, y, z')

y*z

Differentiation

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

g.differentiate(x, 2)

Simplification

f = x * y / (x^2 + y^2 )
z = diff(f, x, y)
z

z.simplify_full()

Substitution

z1=z.substitute(x=2*y^2+1); z

z1.simplify_full()

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

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

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

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

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

• (e) None of the above?

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 $$\frac{x}{e^x - 1} = \sum_{m \ge 0} \frac{B_m x^m}{m!}$$

F = x / (e^x - 1); F  # The generating function

Try to find the first Taylor coefficient:

g = derivative(F, x, 1); g  # another name for .diff()

print(g.simplify_full())
# yes - there are other types of simplifications.
g.substitute(x=0)

We can't just divide 0 / 0. We need L'Hopital's rule!

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

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

bernoulli(1)

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

F.taylor(x, 0, 10)

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

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

q = p.add_bigoh(12); q

print(q.parent())
type(q)

x = q.parent().gen()
x

q + (x + 1).add_bigoh(8)

We can get coefficients of certain terms in Taylor expansions

F.taylor(x, 0, 10).coefficient(x^4)

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

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

F=GF(3)['x, y']
F

x,y=F.gens(); x+4*y

f=x*y+2*x**2
g=x**2+y*x**3

f.gcd(g)

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

m=Matrix(ZZ,[[1,1],[2,3]]); m

m.characteristic_polynomial()

m.eigenvalues() ## Root finding in the above

Algebraic numbers

What is the "?" ?

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

Can be evaluated to arbitrary precision

m.eigenvalues()[0].n(1000)

Can use Numpy

m.numpy()

import numpy
numpy.linalg.eigvals(m)

Can work over finite fields

m=Matrix(GF(5),[[1,1],[2,3]]); m

m.eigenvalues() # "exact"

m.eigenvalues()[0].parent()

Does it work with interval arithmetic?

RIF=RealIntervalField()
h=1e-10
x=(1-h,1+h)
x

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

m.eigenvalues() # Not implemented!

Can use Newton - Raphson to find the roots:

p=m.charpoly()
p.newton_raphson(10,0.2)

p.newton_raphson(10,3)

_[-1].lower(),_[-1].upper()