We will introduce some of the fundamental concepts in Python and 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:
For a full list see: https://www.sagemath.org/links-components.html
It can be integrated with existing commercial packages like
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
Some packages that are not shipped with SageMath are available to install afterwards.
sage -i <package_name>
- to install optional Sage packagesage --pip install <py_package_name>
- to install extra Python package using pip
sage: installed_packages()
- to list all installed packagessage: optional_packages()
- to list installed and available optional packagesFor 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`
On your own machine
For the Sage REPL (read-eval-print loop) or command-line interactive use:
sage
to run Sagesage --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
<function_name>?
- to read documentation (docstring)<function_name>??
- to read the source (if available)<object_name>.<tab>
- use tab completion to see available propertiesExample: Try using these features for yourself below. Run a code cell by first selecting it and then either
Run
button in the toolbar abovegcd?
x=1
x.
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.
int
, float
, complex
str
list
tuple
dict
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))
In SageMath defaults for many numeric types are overwritten. For example:
Integer
type, which has different behaviour from Python's int
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 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))
SageMath stores variables to higher levels of precision than Python
RR
keyword denotes the Real Field with 53 bits of precision (double precision)CC
keyword denotes the Complex Field with 53 bits of precisionRealField(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
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)
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:
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
Ei
for the exponential integral function# 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)
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")
# 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
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:
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 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))
Standard (used in most programming languages):
for
and while
loopsMore Python specific:
# 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
for
loops and list comprehensionsWe 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
plot
functionThe 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)
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:
When the notebook is started we have $x$ as a variable:
x
We can define additional variables
var('x, y, z')
y*z
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)
f = x * y / (x^2 + y^2 ) z = diff(f, x, y) z
z.simplify_full()
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)
m=Matrix(ZZ,[[1,1],[2,3]]); m
m.characteristic_polynomial()
m.eigenvalues() ## Root finding in the above
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()