2. Introduction to SageMath: A Computer Algebra System for All¶
2.1. 1. What is SageMath?¶
An open source system for advanced mathematics.
An open source mathematics distribution (like Linux) with Python as the glue.
A tool for learning and teaching mathematics.
A tool for mathematics research.
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:
Groups, Algorithms, Programming (GAP) - group theory
PARI - rings, finite fields, field extensions, number theory
Singular - commutative algebra
NumPy/SciPy - numerical linear algebra, scientific computing
Integer Matrix Library (IML) - integer, rational matrices
CVXOPT - linear programming, optimization
NetworkX - graph theory
Pynac - symbolic manipulation
Maxima - calculus, differential equations
R - for statistical computing
Ore Algebra - computations with Ore operators
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
andfactor
.
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 and NOT an approximation like .
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 digits of
Sage can compute numerical approximation of numbers to any large precision. Let us increase the precision of the rational number 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 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 differ by .
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 or . These are the kind of expressions that mathematicians and physicists are used to. The variables , , and 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 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 here stands for .
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 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 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 . Here is what I
mean: Let’s build a complex expression with the symbolic variable
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 . Let us create the symbolic variables
.
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 algebraMaple
,Mathematica
orMaxima
. 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 is defined.
parent(M)
I can compute the determinant of .
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 .
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.
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.
SageMath can also handle definite integral. For example:
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 defined in the previous cell. As an example, let’s count the number of spanning trees of the complete graph .
K_3.spanning_trees_count()
SageMath cannot only count the number of spanning trees of the complete graph 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 >