2. Introduction to SageMath: A Computer Algebra System for All

2.1. 1. What is SageMath?

  1. An open source system for advanced mathematics.

  2. An open source mathematics distribution (like Linux) with Python as the glue.

  3. A tool for learning and teaching mathematics.

  4. A tool for mathematics research.

image info

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

2.1.1. 1.1. A Brief Overview

  • Created in 2005 by William Stein.

  • Free and open, GPL license.

  • Includes about 100 open source packages.

  • Now has around 500,000+ lines of new code, by several hundred mathematician-programmers.

Some of the 100 packages includes:

Click here to see all software packages that constitute SageMath.

2.1.2. 1.2 Using Other Non-Commercial Softwares

As mentioned earlier, SageMath includes about 100 open source packages and these packages can be run separately in either a SageMath worksheet or in a Jupyter notebook. There are several ways one can use this. I will only demonstrate a few.

S_8 = gap('Group( (1,2), (1,2,3,4,5,6,7,8) )'); S_8
%%gap

S := Group( (1,2), (1,2,3,4,5,6,7,8) )
maxima('lsum(x^i, i, [1, 2, 7])')
%%maxima

lsum (x^i, i, [1, 2, 7]);
%%python

print(map(lambda x : x**2, [1,3,5])) # note the operation ** and not the same as ^ in Python but in Sage they are.
%%singular

ring R = 0,(x,y,z), lp; R;
%%r
library("MASS")
data(Cars93)
mean(Cars93$MPG.city)
xtabs( ~ Origin + MPG.city, data=Cars93)
%%r

d = c(1,2,3,4)
%%octave
rand (3, 2) # This is will output a 3x2 matrix whose entries are randomly generated between 0 and 1.

Click here if you want to go back to the table of content of this notebook. Otherwise, continue to the next Section.

2.2. 2 Help Inside SageMath

There are various ways to get help in SageMath. Here are several common ways to get help as you are working in the Jupyter notebook.

2.2.1. 2.1 Documentation

SageMath includes extensive documentation covering thousands of functions, with many examples, tutorials, and other helps.

  • One way to access these is to click the “Help” link at the top right of any worksheet. This page has lots of useful commands for the notebook. At the top it has a list of documents; you can click your preferred option at the top of the help page.

  • They are also available any time online at the SageMath website, which has many other links, like video introductions.

  • The Quick Reference Cards are another useful tool once you get more familiar with SageMath.

2.2.2. 2.2 Tab Completion

The most useful help available in the notebook is “tab completion”. The idea is that even if you are not one hundred percent sure of the name of a command, the first few letters should still be enough to help find it. Here’s an example.

  • Suppose you want to do a specific type of plot - maybe a slope field plot - but aren’t quite sure what will do it.

  • Still, it seems reasonable that the command might start with pl .

  • Then one can type pl in an input cell, and then press the Tab key to see all the commands that start with the letters pl .

Try tabbing after the pl in the following cell to see all the commands that start with the letters pl. You should see that plot_slope_field is one of them.

pl

Click here if you want to go back to the table of content of this notebook. Otherwise, continue to the next Section.

2.3. 3 Complete Access to Source Code

Unlike other commercial softwares, SageMath gives its users a complete access to the source code. For example let’s see the source code for the SageMath function is_square and is_squarefree. Both are in the same file. All we need to do is to type the name of the function and after that a double question mark (??) and press the Enter key.

Exercise 1:

Find the source code for the SageMath functions: derivative and factor.

Click here if you want to go back to the table of content of this notebook. Otherwise, continue to the next Section.

2.4. 4 Just the Basics

Why do we need computers for math, anyway? Hopefully by now you and I don’t need help with the following computations:

##%display latex # this is just for since print out.
4 + 2
-2^2
(-2)^2

2.4.1. 4.1 Sage as a Calculator

Like every scientific calculator, SageMath adheres to the standard order of operations, PEMDAS (parenthesis, exponents, multiplication, division, addition, subtraction).

2+3*4^6-45/21
numerical_approx(86015/7)

Instead of using the numerical approximation function numerical_approx, we could also use its alias n for short.

(86015/7).n()

Observe that the above division, after simplification, is the rational number \frac{86015}{7} and NOT an approximation like 12287.8571428571.

Continuing, there is no limit to the size of integers or rational numbers, unless otherwise due to the available memory of the computer in use.

Exercise 2:

Find the numerical approximation to 200 digits of 2^{300}

Sage can compute numerical approximation of numbers to any large precision. Let us increase the precision of the rational number \frac{11}{7} to show the periodicity of its digit expansion.

Can we confirm if the number of digits is really 100? Yes we can, see the code below.

sum(map(len, str(numerical_approx(11/7, digits=100)).split(".")))

Very impressive! Do not panic if you do not understand the above code. Keep calm and relax, hopefully by the end of the session, you should be able to understand it.

Continuing, the operators // and % returns the quotient and the remainder of the division of two integers respectively.

40 // 12
40 % 12

The output of the two operators // and % can also be obtained directly with the divmod function.

divmod(40, 12)

There are several in-built functions that takes integers as inputs. Two of such functions are the factorial and the binomial.

factorial(1000)

Wow! That was very fast! You may not be impressed by this, but this computation took less than a second. To check it, we can either use timeit or time functions.

timeit("factorial(1000)")
time(factorial(1000))

Very impressive! If you are still not impressed, then in Exercise 15, I ask you to write a recursive Python function that returns the factorial of any positive integer and compare the time it takes your function to compute 100! with that of SageMath.

We can also decompose an integer into prime factors. Let us factor the integer 1729.

factor(1729)

Do you know anything about this number 1729? It is a taxicab number and there is a story behind it which goes like this:

Hardy said that it was just a boring number: 1729. Ramanujan replied that 1729 was not a boring number at all: it was a very interesting one. He explained that it was the smallest number that could be expressed by the sum of two cubes in two different ways.

I also observed that the consecutive prime factors of the taxicab number 1729 differ by 6.

Click here if you want to go back to the table of content of this notebook. Otherwise, continue to the next Subsection.

2.4.2. 4.2 Elementary Functions and Usual Constants

Some of the usual mathematical constant in sage includes pi, e, golden_ratio, euler_gamma, catalan. In other to see their real values, we will have to use Sage functions numerical_approx or its short form n.

pi.n()
numerical_approx(e, digits=10)
golden_ratio.n(digits=15)
euler_gamma.n()
catalan.n(digits=5)

There are also elementary functions like: - the exponential function exp; - the logarithmic function log; - trigonometric functions and their inverses: sin, cos, tan, and arcsin, arccos, arctan respectively; - inverse trigonometric functions and their inverses sec, csc, cot and sec, csc, cot respectively; - hyperbolic functions and their inverses: sinh, cosh, tanh, and arcsinh, arccosh, arctanh respectively; - inverse hyperbolic functions and their inverses: sech, csch, coth and arcsech, arccsch, arccoth respectively;

Let’s look at a few of them in action. Note that the computations here are again exact. Recall that you can use the SageMath function numerical_approx or its alias n to get their numerical approximation.

cos(pi)
sin(pi/3)
arccos(0)
exp(3 * I * pi)
exp(oo)
exp(-oo)
arccos(sin(pi/3))
sqrt(2)

As you may have observed in the previous cell, one does not always get the expected results. Indeed, only few simplifications are done by SageMath automatically. If needed, it is possible to explicitly call the simplification function simplify.

simplify(arccos(sin(pi/3)))

There is more we can explore about simplification of symbolic expression. We can talk about this later if there is time.

Click here if you want to go back to the table of content of this notebook. Otherwise, continue to the next Subsection.

2.4.3. 4.3 Python Variables

So far, none of the results of our computations have been saved. So if we want to reuse any result of the previous computations, we will have to either copy and paste it or type them again. This can be very tedious. However, we can avoid all this by using the assignment operator, = to assign the result to a variable. You should already be familiar with this from your Python class.

y = 2 * sin(pi/3)

y is a Python variable. Note that in general, the result of a computation is not automatically printed when it is assigned to a variable. In other to print the content of a variable in the same cell it is assigned, we will have to do the following.

y = 2 * sin(pi/3); y

To reuse later, we only have to type y.

(y + sqrt(2)) * 1/y

Additionally, SageMath saves the last three results in the special variables _, __ and ___:

2 * sin(pi/3)
(_ + sqrt(2)) * 1/_
_ * __
___ - _ * __

You already know from the Programming with Python lectures that you should avoid using keywords as variables. It is a bad practise to redefine predefined constants and functions in SageMath . Although doing this does not influence the internal behaviour of SageMath, it could yield surprising results.

len = 3; len
len([3,3])
pi = -I/3; pi
exp(3*I*pi)

To restore the original value, one can type for example:

from sage.all import pi

Let’s check if this corrects our previous computation.

exp(3*I*pi)

Alternatively, you can also call the restore() function. Doing this will restore all predefined variables and functions to their default values. Let us try this on the previous example.

pi = -I/3; pi
exp(3*I*pi)
restore()
exp(3*I*pi)

But this does not correct the Python function len. Let us try to find the length of the list object [1,3,5].

len([1,3,5])

To correct this you will have to call the function reset(). However it is importatnt to note that calling reset() does a complete reset. That is, it clears all user-defined variables. Let’s call reset() to reset the Python function len to its default definition and test it by finding the length of the list object [1,3,5].

reset()
len([1,3,5])

Click here if you want to go back to the table of content of this notebook. Otherwise, continue to the next Subsection.

2.4.4. 4.4 Symbolic Variables

SageMath is especially useful in dealing with expressions containing variables like x^{2} - y^{2}\,z^{2} + z^{3} or \sin^{3}(x) - 2\,\cos^{2}(x) + \sin(x) - 1. These are the kind of expressions that mathematicians and physicists are used to. The variables x, y, and z that appear in those expressions are symbolic variables. In general, they differ from the Python variables discussed in the previous subsection. In SageMath, the symbolic variables should be explicitly declared before being used except for the symbolic variable x which is already defined. To declare a symbolic variable, we do the following:

y = SR.var('y'); y

We can also use double quote instead of single quote. That is

y = SR.var("y"); y

The \texttt{SR} here stands for \texttt{Symbolic Ring}. Let me briefly explain what exactly is happening in the commands just above. The command SR.var('y') or SR.var("y") builds and returns a symbolic variable which is then assigned to the Python variable y on the left hand side. Note that we could have use any other name for the Python variable instead of y. That is, we could also have done

z = SR.var("y"); z

Thus, assigning the symbolic variable y to the Python variable y is just a convention, which is however recommended to avoid confusion. Let us see if the Python variables y and z are the same in memory.

id(y) == id(z)

Now that we have the symbolic variable y saved in the Python variables y or z, we can use it to build more complex symbolic expressions.

y^(2/3)+ 4*y^2 + sin(y*pi)
z^(3/2)+ 6*z^(1/2) + exp(z*pi)

Here is something important to note: The Python variables y or z does not interact with the symbolic variable y. Here is what I mean: Let’s build a complex expression with the symbolic variable y and assign it to a Python variable. Here, I choose to call my Python variable my_expr.

my_expr = y + 3; my_expr

Now let’s assign the Python variable z to 1.

z = 1; z

What do you think the value of the expression my_expr is?

my_expr

You may be surprised at the result of the last cell. Perhaps you taught it should return 4 right? But if this is the case, how can we evaluate any symbolic expression for example the symbolic expression assigned to the Python variable my_expr? In other to do this, we use the substitution operation.

my_expr(1)
my_expr(y=1)

There is more I can say about the substitution operation but I will leave that one for you to explore.

Continuing, what if you would like to create a large number of symbolic variables. This can be tedious as you may think. However, you can use the shortcut x = SR.var(x, n) where n a positive number which is the number of symbolic variables you would like to create. Note that the indexing starts from 0. Let us create the symbolic variables t_{0}, t_{1}, \dots, t_{199}.

t = SR.var("t", 200); # This creates a tuple of 200 symbolic variables.
(t[0]^3 + t[4]) * t[199]

Use the show function to get a pretty display of the output.

show(_)

In other computer algebra systems like Maple, Mathematica or Maxima, one does not declare symbolic variables before using them. As you have seen already in this Subsection, this is not so with SageMath. However, when one uses the old SageNB notebook interface (which is obsolete), and typing the command, automatic_names(True), one can make SageMath behave like the computer algebra systems I mentioned earlier. That is, one can decide not to declare symbolic variables before using them and SageMath will not complain. However, this do not work when one uses the Jupyter notebook. But it does not mean that it can not be done. In fact, if you are able to this, you will be improving on SageMath and you will become famous. Therefore, I propose this exercise which is a project in itself.

Exercise 3: Hard

Build a SageMath package so that if one types in the command automatic_names(True) when one is running SageMath in the Jupyter notebook, SageMath works like the old SageNB notebook or like the computer algebra Maple, Mathematica or Maxima. In other words, one does not have to declare symbolic variables before using them.

Click here if you want to go back to the table of content of this notebook. Otherwise, continue to the next Subsection.

2.4.5. 4.5 Linear Algebra — Some Basics

What about matrix algebra? Can we perform computations with matrices in SageMath?

M = matrix([[1, 2, 3], [4, 5, 7], [6, 8, 9]]); M

SageMath can tell me the domain over which the entries of the matrix M is defined.

parent(M)

I can compute the determinant of M.

det(M)

I can also compute the inverse.

M^-1

I can also specify the domain over which the entries of a matrix should be defined. For example I want to define a matrix whose elements are from the finite field \mathbb{Z}_{5}.

A = matrix(GF(5), 3, [1,3,5,2,4,9,-1,-3,7]); A
det(A)

Click here if you want to go back to the table of content of this notebook. Otherwise, continue to the next Subsection.

2.4.6. 4.6 Calculus — Some Basics

Many of you are in the applied sciences and will need integrals. But even look up tables can have errors (and often had many more before they were checked against computer integration).  We can use SageMath to evaluate some complicated integrals. Here are a few examples.

\int_{1}^{2} \frac{3 \, x}{\sqrt{25\,x^{2}-3}} \, {\rm d} x, \int_{1}^{e} x \, \log(x) \, {\rm d} x, \int e^{x} \, \log(x) \, {\rm d} x.

integrate( log(x)*x, x, 1, e)
integrate( log(x)*e^x, x)
integral( 3*x^2/sqrt(25*x^2-3), x)

Hmm, the computation is okay, but nowadays computers look nicer than that!  This antiderivative looks kind of ugly.  Let me try to simplify it and typeset it nicely.

show(expand(integral(3*x^2/sqrt(25*x^2-3),x)))

That’s better.  Well, maybe I should just type it in my notes, so I don’t forget.

\frac{3}{50}\sqrt{25x^2-3x}+\frac{9}{250}\log\left(50x+10\sqrt{25x^2-3}\right)

SageMath can also handle definite integral. For example:

\int_{1}^{2} \frac{3 \, x}{\sqrt{25\,x^{2}-3}} {\rm d}x.

integral(3*x^2/sqrt(25*x^2-3),x,1,2)
show(integral(3*x^2/sqrt(25*x^2-3),x,1,2))

And maybe I needed it numerically approximated, with a built-in error bound computed.

numerical_integral( 3*x^2/sqrt(25*x^2-3), 1,2)

And maybe I need to visualize that as well…

show(plot(3*x^2/sqrt(25*x^2-3),(x,1,2),fill=True)+plot(3*x^2/sqrt(25*x^2-3),(x,.5,3)),figsize=3)
html("This shows $\int_1^2 \\frac{3\,x^2}{\sqrt{25\,x^2-3}}\; \mathrm{d}x$")
show(plot(3*x^2/sqrt(25*x^2-3),(x,1,2),fill=True)+plot(3*x^2/sqrt(25*x^2-3),(x,.5,3)),figsize=3, title="This shows $\int_1^2 \\frac{3x^2}{\sqrt{25x^2-3}}\; \mathrm{d}x$")

Even if you haven’t done these sorts of things yourself, you know you can do so with the help of a computer, and you possibly have something in the computer lab or on your laptop which can do them - or you use some web applications or other softwares.

Click here if you want to go back to the table of content of this notebook. Otherwise, continue to the next Subsection.

2.4.7. 4.7 Plotting in Brief

SageMath is not only useful in dealing with symbolic expressions but also with plotting. There are lots of plotting commands in SageMath. In the previous Subsection, we saw some examples of the SageMath plot function. In addition, to the plot function, I will give examples of the SageMath functions: plot3d and parametric_plot3d. As you may have have observed in the previous Subsection, the plot command makes it easy to draw the curve of a real function on a given interval. On the other hand, the plot3d command is its counterpart for three-dimensional graphics, or for the graph of a real function of two variables. A similar description applies to the parametric_plot3d. Here are examples of those three commands:

plot(cos(sqrt(2) * x), x, -pi, pi)
plot3d(cos(pi*sqrt(x^2 + y^2))/sqrt(x^2+y^2),(x,-5,5), (y,-5,5))
u, v = var('u, v')
f1 = (4+(3+cos(v))*sin(u), 4+(3+cos(v))*cos(u), 4+sin(v))
f2 = (8+(3+cos(v))*cos(u), 3+sin(v), 4+(3+cos(v))*sin(u))
p1 = parametric_plot3d(f1, (u,0,2*pi), (v,0,2*pi), texture="red")
p2 = parametric_plot3d(f2, (u,0,2*pi), (v,0,2*pi), texture="blue")
show(p1 + p2)

Click here if you want to go back to the table of content of this notebook. Otherwise, continue to the next Subsection.

2.4.8. 5 Demonstration of Some SageMath Components

One reason SageMath is so powerful is that it uses the “best of” other softwares. Click on this link to see the complete list of all software packages that constitutes SageMath.  We’ve already seen it contains R, and it also includes Numpy/SciPy, which I believe you were introduced to in your Python course.

Let us see a more theoretical example of this.  It isn’t important whether you have seen any of the things here before, just that you can see us using many subsystems.

K_3 = graphs.CompleteGraph(3); K_3.plot(figsize=2) # Takes NetworkX graph

I can use SageMath to compute some properties of the complete graph K_{3} defined in the previous cell. As an example, let’s count the number of spanning trees of the complete graph K_{3}.

K_3.spanning_trees_count()

SageMath cannot only count the number of spanning trees of the complete graph K_{3} but it can plot them as well.

M = Matroid(K_3);
M_bases = M.bases();
[K_3.subgraph(edges = list(M_bases[basis])).show(figsize=1.5) for basis in range(len(M_bases))]
H = K_3.cartesian_product(K_3) # Uses SageMath's functionality
show(H) # Uses matplotlib
Sym = H.automorphism_group() # Native Sage Code using C backend for graphs
len( Sym.list() )  # Needs GAP to calculate
Sym._gap_().ComputedPCentralSeriess()  # Directly within GAP
chrom_poly = H.chromatic_polynomial(); chrom_poly # Uses compiled C code from Cython in SageMath
deriv = derivative(chrom_poly, x); deriv # Uses Ginac/Pynac for symbolics

Believe it or not, this has meaning! It tells us the graph is connected.

deriv.subs(x=0)
H.is_connected()
integral(chrom_poly,x) # Uses Maxima for integration, symbolic summation, etc.

No one knows if this has any meaning.  The book “Quo Vadis, Graph Theory?: A Source Book for Challenges and Directions” asks for interpretations of such things.

Click here if you want to go back to the table of content of this notebook.


< Table of contents | 2. Review of Python programming >