Notebook 1 - An Introduction to Sage.ipynb

Kernel: SageMath 9.1

A Brief Introduction to Sage Math

This document aims to give a crash-course to Sage. There are many additional resources for help, including the built-in documentation (discussed below), the official Sage tutorial, and the (highly recommended) open textbook Computational Mathematics with SageMath.

Sage is free and open source. Information on running a local installation can be found on the Sage installation guide. Alternatively, Sage can be run "in the cloud" by making a (free) account on the CoCalc website.

This document is written as a Jupyer notebook, the most common (and convenient) way to write and execute Sage code. A notebook is composed of cells. Most of the cells in this notebook consist of an Input section (containing Sage code) and (potentially) an output section (containing the result of evaluating that Sage code) $-$ some code cells simply perform a computation without returning anything (for instance, updating the values of variables). A few cells (including the current one) consist of formatted text and LaTeX equations, written using the Markdown markup language. A third type of cell contains plain, unformatted text.

To execute a piece of Sage code, click on the Input section of the corresponding code cell and hit Shift + Enter (only hitting Enter simply adds a new line). The reader should execute each statement as they work through the notebook, and is encouraged to modify the code and play around as they go. Note that skipping a cell may result in errors when later cells are executed (for instance, if one skips a code block defining a variable and later tries to run code calling that variable). There are a selection of short exercises throughout, and a few larger exercises in the final section. To add a new cell, click to the left of any cell and press the "a" key. To delete a cell, click to the left of a cell and press the "d" key. These (and other) tasks can also be accomplished through the menu bars at the top of the page.

Additional details on the topics most closely related to combinatorics are covered in a follow-up notebook, available by clicking here.

Part 1: The Basics of Sage

We begin with an explanation of arithmetic in Sage, if statements, for and while loops, and Sage functions (in the programming sense; symbolic mathematical functions are described in Part 2 below).

# Any text on a line after a '#' symbol is considered a comment, and is not evaluated by Sage
# Sage can be used as a calculator using usual math notation (recall you need it hit Shift + Enter to evaluate a cell)
1 + 3*4

13

# Code in a Sage cell should be entered in sequence, with one "instruction" per line
# (multiple instructions can be put on the same line using a semicolon -- see below)
# The result of previous executions of Sage cells (both the current cell and other cells) is stored by Sage
# All (non-commented) code is executed, however only the output of the final line is printed by default
1 + 2
2 ^ 3

8

# To print other lines, use the print command
print(2^3) # This is an integer
print(5/10) # This is an exact rational number
print(5.0/10) # This is a floating point number
print(11 % 3) # This is 11 mod 3

8
1/2
0.500000000000000
2

# Multiple instructions can be put on the same line using a semicolon
# Again, only the output of the final line is displayed by default
2+2
print(1+2); 3*5

3

15

# The "show" command outputs a latex-like pretty output
# Sage knows common math constants such as pi (lower case), e, and I (or the alternative i)
# Most common mathematical functions are supported by default
show(5/10 + 1/pi)
show(sin(pi/3))
show(log(2))
show(log(2).n()) # Adding ".n()" gives a numerical approximation (can use ".n(k)" to get k bits)
show(exp(i*2*pi/3))
show(exp(I*2*pi/101))

# Mathematical objects and variables can be inserted in print commands using {} and .format(math1,math2,...)
print("The area of a circle of radius 1 is approximately",pi.n())
print("The square-root of {} is approximately {}".format(log(2),sqrt(log(2)).n()))

The area of a circle of radius 1 is approximately 3.14159265358979
The square-root of log(2) is approximately 0.832554611157698

# To access the help page for a function, type the name of the function and then add a question mark
# For instance, evaluating the expression "sin?" (without quotes) gives the help page for sine
# THIS WILL OPEN A NEW WINDOW AREA
# Using two question marks (e.g. "sin??") shows the source code of the function
# (Note that for many low level functions like sin Sage relies on outside packages and the source code is not very informative)
sin?

# Variables are defined using the "=" operator
# Note that some operations (such as variable assignment) do not have an output
# So we add another line to print the value of our variable
mult = 11 * 12
mult

132

# The values "_", "__", and "___" (using one, two, and three underscores) are the last three output values
# Printing displays content but does *not* return a value, and a block ending with a print statement has no output
print(_)

132

# The output of previous cells can be directly accessed using the syntax Out[k]
# Note cells are numbered by execution order: running cells repeatedly or in different orders will change their numbering
Out[8]

132

# Sage is built on Python, and uses much of the same syntax. In particular, indentation is extremely important.
# If statements are defined with indentation
a = 2

if a<0:
    print("Negative")
elif a==0:
    print("Zero")
else:
    print("Positive")

Positive

# In addition, an if statement can be written in one line
a = 2
if a>0: print("Positive")

Positive

# Functions are also defined with indentation
def fact(n):
    if n==0:
        return 1;
    else:
        return n*fact(n-1)

print(fact(4))
print(fact(10))
print(fact) # Prints reference to the function
# fact?? # uncommenting and running this code displays the source code for the function fact, just as for built-in functions

24
3628800
<function fact at 0x7f2e610a9950>

################################################
# EXERCISE: Make a function to compute the nth Fibonacci number (defined by fib(n+2) = fib(n+1) + fib(n)
# and fib(0)=fib(1)=1). What is the highest number you can compute in 5 seconds?
#
# Adding the line "@cached_function" directly before the function definition tells Sage to store the 
# result of each function return, which will speed up the computation. Add this line then see how high
# you can go in 5 seconds.
################################################

# Lists in Sage are defined with square brackets
LST = [1,2,3,4]

print(LST)
print(len(LST))
print(LST[0]) # Access elements from the left (starting at index 0)
print(LST[1])
print(LST[-2]) # Access elements from the right
print(LST[1:3]) # Define sublist
print(LST[1:-1])
print(LST[1:]) # Define sublist from start to end of list

[1, 2, 3, 4]
4
1
2
3
[2, 3]
[2, 3]
[2, 3, 4]

# Lists can be concatenated with '+'
# Strings can be concatenated similarly (they are lists of characters)
print([1,2,3] + ["a","b","c"])
print("hello" + " " + "world")

[1, 2, 3, 'a', 'b', 'c']
hello world

# For loops work over lists or generators, and are indented similarly to if statements and functions
LST = [0,1,2,3,4,5]
for k in LST:
    print(k^2)

0
1
4
9
16
25

# The notation a..b can be used to build the list of numbers between integers a and b
for k in [0..5]:
    print(k^2)

0
1
4
9
16
25

# As in Python, Sage contains the function range(k), which encodes the numbers 0,1,...,k-1
# In Sage versions 8.x and lower, which rely on Python 2, range(k) returns the list [0,1,...,k-1]
# In Sage 9.0 and higher, which use Python 3, range(k) returns an *iterator* that can be converted to a list
# (think of an iterator as a method to construct elements of a list, one by one, which can be more efficient)
# More details on this change can be found here: https://docs.python.org/3.0/whatsnew/3.0.html#views-and-iterators-instead-of-lists

print(range(5)) # printing an iterator just returns the iterator
print(list(range(5))) # the iterator can be converted to a regular list
for k in range(5): # loops can range directly over iterators, which can be more efficient
    print(k^2)

range(0, 5)
[0, 1, 2, 3, 4]
0
1
4
9
16

# For loops can also go over one line
for k in range(5): print(k^2)

0
1
4
9
16

# There is a powerful alternate syntax for building lists using a function f(x) on the elements of a list LST -- [f(k) for k in LST]
print([k^2 for k in range(5)])
print([cos(k*pi) for k in [1..5]])

[0, 1, 4, 9, 16]
[-1, 1, -1, 1, -1]

################################################
# EXERCISE: Compute the sum of the first 100 perfect squares. (Use the add function -- type "add?" for its documentation)
################################################

# While loops are defined similarly
k = 0
while k<5:
    print(k^2)
    k = k+1

0
1
4
9
16

# While loops can be broken using 'break'
k = 0
while True:
    if k>= 5:
        break
    print(k^2)
    k = k+1

0
1
4
9
16

# The map operator applies a function to each element of a list 
# Similar to range, in Sage 9.0+ map returns a "map object" / iterator
# The list function can be applied to obtain an honest list from a map object
list(map(abs,[-1,2,-3,4]))

[1, 2, 3, 4]

# User defined functions can be mapped, where appropriate
def myfun(k):
    return 2*k

list(map(myfun,[-1,2,-3,4]))

[-2, 4, -6, 8]

# Can also use map with 'lambda expressions' to define a function in place
print(list(map(lambda t: t^2, [-1,2,-3,4])))
print(lambda t: t^2) # Defines the function f(x) = x^2 in place

[1, 4, 9, 16]
<function <lambda> at 0x7f2e610a97b8>

# Can filter a list using 'filter'
# Similar to range and map, in Sage 9.0+ filter returns a "filter object" / iterator
# The list function can be applied to obtain an honest list from a filter object

print(filter(lambda t: t>0, [-1,2,-3,4]))
print(list(filter(lambda t: t>0, [-1,2,-3,4])))

<filter object at 0x7f2e60fa94a8>
[2, 4]

# Can also use the list 'comprehension form' to filter elements when building a list
[p^2 for p in [1..10] if is_prime(p)]

[4, 9, 25, 49]

# Can sort lists, when appropriate -- sort is a function of a list (see Part 2 below for additional details)

L = [1,4,3,2,7,6]
L.sort() # Modifies the list L in place
print(L)
L = ["a","acb","abc","ab"]
L.sort() # Sort in lexicographical order
print(L)

[1, 2, 3, 4, 6, 7]
['a', 'ab', 'abc', 'acb']

# Lists are stored by reference. Simply writing L2 = L1, makes L1 and L2 point to the *same* list
# Use L2 = copy(L1) to make an independent copy
# Note copy only works at one level of lists (use deepcopy to copy a list of lists, and make everything independent)

L1 = [1,2,3]
L2 = L1 # L2 points to the same list as L1
L3 = copy(L1) # L3 points to a new list, initialized with the same values as L1
L2[0]=-1
print(L1); print(L2); print(L3)

[-1, 2, 3]
[-1, 2, 3]
[1, 2, 3]

# Sage also supports dictionaries, which for the user behave as lists indexed by strings
L = {}; L['apple'] = 1; L['pear'] = pi; L['banana'] = sin(2)
show(L)
show(L['apple'] + L['pear'])
del(L['apple'])
show(L)

###########################################################################
# EXERCISE: Create a list containing the first 100 primes of the form 4*k+1
###########################################################################

Part 2: The Symbolic Ring and Sage Types

We now see how to manipulate symbolic variables and abstract functions, including basic calculus operations and plotting, and how to determine the type and parent of an object.

# Before running this section, we reset all variables
reset()

# By default, when opening a notebook the variable "x" can be used to define a symbolic function / expression
poly = x^2 - 1
print(poly)

x^2 - 1

# Using any other (undeclared) variable will give an error
poly2 = y^2 - 1

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
/ext/sage/sage-9.1_1804/local/lib/python3.7/site-packages/sage/all_cmdline.py in <module>()
      1 # Using any other (undeclared) variable will give an error
----> 2 poly2 = y**Integer(2) - Integer(1)

NameError: name 'y' is not defined

# If the variable x is assigned a different value, this does not change the value of our symbolic expression!
# Of course, any new symbolic expressions will use this updated value of x.
x = 2
print(poly)
poly2 = x^2 - 1
print(poly2)

x^2 - 1
3

MAKE SURE YOU UNDERSTAND THIS CRUCIAL POINT: This behaviour occurs because Sage variables (for instance, x on the left hand side of x = 2) are distinct from the underlying symbolic variables used to define symbolic expressions. By default the Sage variable x is initialized to a symbolic variable "x", and the expression poly above is defined in terms of this symbolic variable. Changing the Sage variable x to a new value does nothing to the underlying symbolic variable "x", which is why the value of poly does not change after setting x = 2.

# The easiest way to define a new symbolic variable, having the same name as a Sage variable, is with the "var" command
x = 2
print(x) # Prints the value of the Sage variable x
var('x') # Makes the Sage variable x point to the symbolic variable "x"
print(x) # Prints the value of the Sage variable x, which is now the symbolic variable "x"

2
x

# Multiple variables can be defined at the same time
var('a b c') # Initialize Sage variables a, b, and c to symbolic variables "a", "b", and "c"
poly2 = (a + b*c)^2
print(poly2)

(b*c + a)^2

##############################################################################################################
# EXERCISE: Guess the result of uncommenting and running the following code. Explain the output that does occur.
##############################################################################################################
# var('p q r'); r = q; q = p; p = r
# print(p); print(q); print(r)

# The commands "type" and "parent" illustrate the domains in which Sage objects live
# Symbolic expressions, defined with symbolic variables, live in the Symbolic Ring
# Sage automatically determines where objects defined with "=" should live
# Part 3 below gives many more details about the types of objects in Sage, and how to specify and manipulate them
var('a')
poly = a
print(type(poly))
print(parent(poly))

<class 'sage.symbolic.expression.Expression'>
Symbolic Ring

# Some additional examples
poly = 10
print(type(poly))
print(parent(poly))
print(type(1/2))
print(parent(1/2))
print(type(0.5))
print(parent(0.5))

<class 'sage.rings.integer.Integer'>
Integer Ring
<class 'sage.rings.rational.Rational'>
Rational Field
<class 'sage.rings.real_mpfr.RealLiteral'>
Real Field with 53 bits of precision

# In Sage (as in Python) objects can have their own functions
# Type "poly2." (without quotes) then hit the *Tab* key to see the available functions for poly2
# Run the command "poly2.subs?" to see the help page for the function poly2.subs
# Run the command "poly2.subs??" to see the source code for the function poly2.subs
var('a b c'); poly2 = (a + b*c)^2
print(poly2.subs(a=1,b=2,c=3))
print(poly2.subs)
print(parent(poly2.subs))

49
<built-in method substitute of sage.symbolic.expression.Expression object at 0x7f2e6106dc88>
<class 'builtin_function_or_method'>

##################################################################################################
# EXERCISE: Uncomment and hit <Tab> with the cursor on the far right of the following line see the
# available functions for poly2 which begin with "s". Select one of the functions, look up its
# documentation, then run the function (with appropriate arguments if necessary)
##################################################################################################
#poly2.s

# Sage has many commands for manipulating expressions, such as simplify and expand.
# Be careful -- in the Symbolic Ring, Sage simplifies without checking restrictions on variables
var('x')
print(((x-1)*(x+1)/(x-1)).simplify())
print(simplify((x-1)*(x+1)/(x-1))) # simplify(p) is a built-in shortcut for type(p).simplify()
print(expand((x-1)*(x+1)))

pol = x^2-1
print(pol.factor())

x + 1
x + 1
x^2 - 1
(x + 1)*(x - 1)

# Equations are also symbolic expressions, defined using "=="
eq1 = x^3 == 1
print(eq1)
print(eq1.lhs())
print(type(eq1))

x^3 == 1
x^3
<class 'sage.symbolic.expression.Expression'>

# The solve command works with equations
show(solve(eq1,x))
show(eq1.solve(x))

# Symbolic functions can be defined with symbolic variables
var('x')
f = sin(x)+2*cos(x)
print(f)
print(parent(f))

2*cos(x) + sin(x)
Symbolic Ring

# The syntax for a "callable" symbolic function is analogous
var('x')
f(x) = sin(x)+2*cos(x)
print(f)
print(parent(f))

x |--> 2*cos(x) + sin(x)
Callable function ring with argument x

# The find_root command can be used to approximate roots numerically
# (additional details on finding roots of polynomials are given in the next section)
f(x).find_root(-10, 10)

8.31762924297529

# Symbolic functions can be plotted with the plot command
plot(f(x),x,0,10) # Syntax is plot(function, variable, xmin value, xmax value)

# Plots can be "added" to overlap -- run "plot?" for plotting options
p1 = plot(sin(x),x,0,pi)
p2 = plot(cos(x),x,0,pi,color="red")
p1+p2

# Common plot options include
# plot_points (default 200)
# xmin and xmax
# color
# detect_poles (vertical asymptotes)
# alpha (line transparency)
# thickness (line thickness)
# linestype (dotted with ':', dashdot with '-.', solid with '-')
# Use commands list p.set_aspect_ratio, etc.

pp = plot([]) # Define empty plot
cols = rainbow(20) # Define 20 evenly spaced colors
for k in range(20):
    pp = pp + plot(log(x+k),1,10,color=cols[k],thickness=2)

pp # Print superimposed graphs

# 3d plots are similar (user may need to load the 3D viewer after running this code)
var('x y')
f(x,y) = x^2 + sin(x)*cos(y)
plot3d(f, (x,-1,1), (y,-1,1))

# The series command computes power series of symbolic expressions representing functions
show(log(1/(1-x)).series(x,10))

# The series command can also compute Laurent series
show(tan(x).series(x==pi/2,10))

# Series are not symbolic expressions, but symbolic series (note the Symbolic Ring is the parent of both)
# This means some methods available to symbolic expressions may not be available for symbolic series
ser = arctan(x).series(x,10)
print(type(ser))

<class 'sage.symbolic.series.SymbolicSeries'>

# To go from a series expression to the polynomial defined by the terms, use the truncate command
pol = arctan(x).series(x,10).truncate()
print(pol)
print(type(pol))

1/9*x^9 - 1/7*x^7 + 1/5*x^5 - 1/3*x^3 + x
<class 'sage.symbolic.expression.Expression'>

# The diff/differentiate command computes derivatives
# The integral/integrate command is similar
print(diff(sin(x),x))
print(integrate(sin(x),x))
print(integrate(exp(-x^2),x,-infinity,infinity))

cos(x)
-cos(x)
sqrt(pi)

# This also works with abstract symbolic functions
var('x y')
function('f')(x)
print(type(f))
print(diff(sin(f(x)),x,x))
function('g')(x,y)
print(diff(g(sin(x),f(x)),x))

<class 'sage.symbolic.function_factory.function_factory.<locals>.NewSymbolicFunction'>
-sin(f(x))*diff(f(x), x)^2 + cos(f(x))*diff(f(x), x, x)
cos(x)*D[0](g)(sin(x), f(x)) + diff(f(x), x)*D[1](g)(sin(x), f(x))

# Here we set up and solve a differential equation
x = var('x'); y = function('y')(x)
desolve(y*diff(y,x) == x, y, show_method=True)

[1/2*y(x)^2 == 1/2*x^2 + _C, 'separable']

# Now compute a Laplace transform
x, s = var('x, s');
print(sin(x).laplace(x,s))
print((1/(s^2 + 1)).inverse_laplace(s,x))

1/(s^2 + 1)
sin(x)

# Sage can calculate some symbolic sums
var('a k n')
show(sum(a^k,k,0,n))
show(sum(binomial(n,k), k, 0, n))

# You can use the assume to command to put restrictions on variables
# This can be useful for infinite series
# Just running sum(a^k,k,0,infinity) throws a "ValueError" -- need to *assume* |a|<1
assume(abs(a)<1)
show(sum(a^k,k,0,infinity))

# Assumptions stay with variables until cleared using forget
# Run this cell directly after the previous cell
print(bool(a<1)) # Test if a < 1, under assumption abs(a) < 1
forget(abs(a)<1)
print(bool(a<1)) # Test if a < 1, under no assumptions

True
False

Part 3: Mathematical and Algebraic Objects (rings and fields, polynomials, and power series)

Working over the symbolic ring should be familiar anyone with experience in the Maple and Mathematica computer algebra packages. However, the separation of Sage variables from symbolic variables allows Sage great flexibility for defining and working with mathematical objects. In this section we discuss the different mathematical domains supported in Sage and how to define and work with mathematical objects.

# Reset all variables
reset()

# We will define mathematical objects using various algebraic domains (rings, fields, vector spaces, etc.)
# Here are a few common domains and their default shortcuts (note that it is possible to accidentically overwrite the default shortcuts!)

print(SR) # Symbolic Ring
print(ZZ) # Integers (ZZ is default shortcut for IntegerRing())
print(QQ) # Rational numbers (QQ is default shortcut for RationalField())
print(QQbar) # Algebraic numbers (QQbar is default shortcut for AlgebraicField())
print(GF(9)) # GF(p^k) is the finite field with p^k elements (GF(p^k) is default shortcut for FiniteField(p^k))
print(AA) # Real algebraic numbers (AA is default shortcut for AlgebraicRealField())
print(IntegerModRing(10)) # IntegerModRing(n) = ring Z/nZ
print(RealField(10)) # RealField(k) is real floating points with k bits of precision (RR is a shortcut for RealField(53))
print(ComplexField(10)) # ComplexField(k) is complex floating points with k bits of precision  (CC is a shortcut for ComplexField(53))

Symbolic Ring
Integer Ring
Rational Field
Algebraic Field
Finite Field in z2 of size 3^2
Algebraic Real Field
Ring of integers modulo 10
Real Field with 10 bits of precision
Complex Field with 10 bits of precision

# Domains can be defined recursively (we do not get into user defined domains here)
# Here we define the Sage variable t to point to the symbolic variable "t" generating a polynomial ring with rational coefficients
# The notation "A.<t>" tells Sage to make the Sage variable t equal the symbolic variable "t"
# Running A = QQ['t'] would make A the same field, but would not update the Sage variable t
A.<t> = QQ['t']
print(A)
print(t^2+1)

# The Sage variable representing the symbolic variable doesn't need to match the symbolic variable
B.<v> = QQ['VAR']
print(B)
print(v^2+1)
# print(VAR) # running this would give an error as the Sage variable VAR is undefined

Univariate Polynomial Ring in t over Rational Field
t^2 + 1
Univariate Polynomial Ring in VAR over Rational Field
VAR^2 + 1

# Two more examples, defining power series and matrix rings
R = IntegerModRing(10)
print(R[['w']])
print(MatrixSpace(R,3,2))

Power Series Ring in w over Ring of integers modulo 10
Full MatrixSpace of 3 by 2 dense matrices over Ring of integers modulo 10

# These domains are themselves objects in Sage, and have associated functions
# If no symbolic variable name is specified, the Sage variable is used for the symbolic variable by default
A.<t> = QQ[]
B.<x> = GF(2)[]
print(A)
print(A.is_commutative())
print(A.is_euclidean_domain())
print(A.gen())
print(B)
print(B.random_element(5))
print(list(B.polynomials(of_degree=3)))

Univariate Polynomial Ring in t over Rational Field
True
True
t
Univariate Polynomial Ring in x over Finite Field of size 2 (using GF2X)
x^5 + x^3 + x^2 + 1
[x^3, x^3 + 1, x^3 + x, x^3 + x + 1, x^3 + x^2, x^3 + x^2 + 1, x^3 + x^2 + x, x^3 + x^2 + x + 1]

# We make the polynomial expressions in t live in the field QQ[t]
# Functions in these expressions, such as roots, factor, etc. can now work over the appropriate domain (or give an error when unsupported)
# By default, the roots command returns the a list of (root, multiplicity)
A.<t> = QQ[]
pol = (t^2 - 1)*(t^2-2)
print(parent(pol))
print(pol.factor())
print(pol.roots())
print(pol.factor) # We see the factor method uses the FLINT package

Univariate Polynomial Ring in t over Rational Field
(t - 1) * (t + 1) * (t^2 - 2)
[(1, 1), (-1, 1)]
<built-in method factor of sage.rings.polynomial.polynomial_rational_flint.Polynomial_rational_flint object at 0x7f2e50c330b8>

# Here we accidentally set the Sage variable t to a constant, instead of the Symbolic variable "t"
# This does not affect our polynomial pol, but hinders defining new polynomials
# We can fix this by setting the Sage variable back to the generator of our field A = QQ[t]
A.<t> = QQ[]; pol = (t^2 - 1)*(t^2-2) # Define t symbolically as above

t = 2
print("The polynomial {} is no longer entered as (t^2 - 1)*(t^2-2), which now evaluates to {}".format(pol,(t^2 - 1)*(t^2-2)))
print("The problem is that the variable t is a constant, now lying in {}".format(parent(t)))

# Set Sage variable t back to symbolic variable "t" generating the ring A
t = A.gen()
print("Now the variable t is back to lying in {}".format(parent(t)))
print(A == parent(t))

The polynomial t^4 - 3*t^2 + 2 is no longer entered as (t^2 - 1)*(t^2-2), which now evaluates to 6
The problem is that the variable t is a constant, now lying in Integer Ring
Now the variable t is back to lying in Univariate Polynomial Ring in t over Rational Field
True

# The domain of an object can be explicitly defined by casting into the domain (when possible)
a = QQ(1/2) # Cast 1/2 into the field of rational numbers
b = QQbar(1/2) # Cast 1/2 into the field of real algebraic numbers
c = RR(1/2) # Cast 1/2 into the field of real floating point numbers to 53 bits of precision
print("We have the constants a = {} (over {}), b = {} (over {}), and c = {} (over {})".format(a,parent(a),b,parent(b),c,parent(c)))

A = ZZ['t']
pol = A('t^2-1') # Define a polynomial by casting from a string to ZZ[t]
pol2 = QQ['t'](pol) # Cast pol from ZZ[t] into QQ[t]
print("Although they look the same, pol = {} lives in {} while pol2 = {} lives in {}".format(pol,parent(pol),pol2,parent(pol2)))

We have the constants a = 1/2 (over Rational Field), b = 1/2 (over Algebraic Field), and c = 0.500000000000000 (over Real Field with 53 bits of precision)
Although they look the same, pol = t^2 - 1 lives in Univariate Polynomial Ring in t over Integer Ring while pol2 = t^2 - 1 lives in Univariate Polynomial Ring in t over Rational Field

print(GF(3)(1/2)) # This conversion is allowed, because 2 has an inverse in Z/3Z
#print(GF(3)(1/3)) # Uncomment and run this code to get an error, because 3 is not invertible in Z/3Z

2

# change_ring and change_variable_name can also be used to change parameters of a polynomial ring
R = ZZ['t']
p = R.random_element(5)
print("p = {} is an element of {}".format(p, parent(p)))
q = p.change_ring(GF(3))
print("q = {} is an element of {}".format(q, parent(q)))
r = p.change_variable_name('T')
print("r = {} is an element of {}".format(r, parent(r)))

p = t^5 + 3*t^4 - t^3 - 1 is an element of Univariate Polynomial Ring in t over Integer Ring
q = t^5 + 2*t^3 + 2 is an element of Univariate Polynomial Ring in t over Finite Field of size 3
r = T^5 + 3*T^4 - T^3 - 1 is an element of Univariate Polynomial Ring in T over Integer Ring

# Can substitute values that can be interpreted as elements of an A-algebra into polys of A[x]
sub1 = p.subs(3)
sub2 = q.subs(t=3) # Variable can be explicitly specified (optional for univariate polys)
sub3 = q.subs(matrix([[1,2],[3,4]]))

print("The value {} obtained from the polynomial {} lies in {}".format(sub1,p,parent(sub1)))
print("The value {} obtained from the polynomial {} lies in {}".format(sub2,p,parent(sub2)))
print("The following matrix, obtained from the polynomial {}, lies in {}".format(p,parent(sub3)))
show(sub3)

The value 458 obtained from the polynomial t^5 + 3*t^4 - t^3 - 1 lies in Integer Ring
The value 2 obtained from the polynomial t^5 + 3*t^4 - t^3 - 1 lies in Finite Field of size 3
The following matrix, obtained from the polynomial t^5 + 3*t^4 - t^3 - 1, lies in Full MatrixSpace of 2 by 2 dense matrices over Finite Field of size 3

# Can apply functions to all coefficients
p.map_coefficients(lambda z: z^2)

t^5 + 9*t^4 + t^3 + 1

# Factorization works over different coefficient rings
A.<x> = QQ[]
p = QQ['x'](x^2-2)
print(p.factor()) # Factor over the rationals
print(p.change_ring(GF(7)).factor()) # Factor in finite field of order 7
print(p.change_ring(QQ[sqrt(2)]).factor()) # Factor over Q adjoin the square-root of 2
print(parent(p.factor())) # Factorizations lie in a Factorization class, which allows one to access the elements of the factorization, the powers, etc.

x^2 - 2
(x + 3) * (x + 4)
(x - sqrt2) * (x + sqrt2)
<class 'sage.structure.factorization.Factorization'>

# The roots method works the same way
# Over the field of algebraic numbers, the solutions are stored exactly.
# Algebraic numbers (among other special classes of numbers) are denoted by a "?" at the end of a decimal
p = QQ['x']((x-1)*(x^3-x+1))
print("The roots of {} in QQ are {} \n".format(p,p.roots(multiplicities=False)))

rts = p.roots(QQbar,multiplicities=False)
print("The roots of {} in QQbar are {} \n".format(p,rts))

r1 = rts[0]
print("For instance, the algebraic number {} is a root of {}, encoded as an element of {}".format(r1,p,parent(r1)))
print("The algebraic number {} is stored exactly. Its minimal polynomial is {} and a numeric approximation to 100 bits is {}".format(r1,r1.minpoly(),r1.n(100)))

The roots of x^4 - x^3 - x^2 + 2*x - 1 in QQ are [1] 

The roots of x^4 - x^3 - x^2 + 2*x - 1 in QQbar are [-1.324717957244746?, 1, 0.6623589786223730? - 0.5622795120623013?*I, 0.6623589786223730? + 0.5622795120623013?*I] 

For instance, the algebraic number -1.324717957244746? is a root of x^4 - x^3 - x^2 + 2*x - 1, encoded as an element of Algebraic Field
The algebraic number -1.324717957244746? is stored exactly. Its minimal polynomial is x^3 - x + 1 and a numeric approximation to 100 bits is -1.3247179572447460259609088545

############################
# Short Exercise: Type "r1." and press <Tab> to the methods available for an algebraic number.
# Run a few of these methods (look at the corresponding help pages, if necessary)
############################

# The structure of the roots can be explored further using Sage's methods
p = QQ['x'](x^3-x+1)
gal = p.galois_group()
print("The Galois group of {} is the group {} \n".format(p,gal.as_finitely_presented_group()))

K.<alpha> = NumberField(x^3-x+1) # Number field where alpha is a root of x^3-x+1
print("{} has the factorization {} in {} \n".format(p,p.change_ring(K).factor(),K))

L.<beta> = p.splitting_field() # Splitting field of p
print("A splitting field of {} is {} \n".format(p,L))
print("Over this splitting field, {} has factorization {}".format(p,p.change_ring(L).factor()))

The Galois group of x^3 - x + 1 is the group Finitely presented group < a, b | a^2, b^3, (b*a)^2 > 

x^3 - x + 1 has the factorization (x - alpha) * (x^2 + alpha*x + alpha^2 - 1) in Number Field in alpha with defining polynomial x^3 - x + 1 

A splitting field of x^3 - x + 1 is Number Field in beta with defining polynomial x^6 + 3*x^5 + 19*x^4 + 35*x^3 + 127*x^2 + 73*x + 271 

Over this splitting field, x^3 - x + 1 has factorization (x - 1/69*beta^5 - 2/69*beta^4 - 13/69*beta^3 - 2/69*beta^2 - 12/23*beta + 100/69) * (x + 3/575*beta^5 + 3/575*beta^4 + 1/575*beta^3 - 147/575*beta^2 - 6/23*beta - 906/575) * (x + 16/1725*beta^5 + 41/1725*beta^4 + 14/75*beta^3 + 491/1725*beta^2 + 18/23*beta + 218/1725)

# The ring of power series works up to some precision, and automatically keeps track of precision
R.<x> = QQ[[]]
f = 1 + 2*x + 3*x^2 + O(x^3)
g = 1 - x + 3*x^2 + O(x^4)
print(f+g)
print(f*g)
print("The precision of [{}] + [{}] is {}".format(f,g,(f+g).prec()))

2 + x + 6*x^2 + O(x^3)
1 + x + 4*x^2 + O(x^3)
The precision of [1 + 2*x + 3*x^2 + O(x^3)] + [1 - x + 3*x^2 + O(x^4)] is 3

# Be careful with equality of formal series
# Equality only tested up to lowest precision
print(1+2*x == 1+x+O(x^2))
print(1+x == 1+x+O(x^2))
print(1+x+x^3 == 1+x+O(x^2))
print(0 == O(x^2))

False
True
True
True

# Sometimes series can be calculated by a direct conversion from the symbolic ring
show(sqrt(1+x) + O(x^10))
show(sin(x+2*x) + O(x^10))

# Ideals and quotient rings can also be defined
R.<x> = QQ[]
J = R.ideal(x^2+1)
print(J)

K = R.quo(J)
xb = K(x) # Make xb the image of x in the residue ring K
print("The generator of the ring [{}] is {}".format(K,xb))
print("In K, xbar^2 = {}".format(xb^2))

Principal ideal (x^2 + 1) of Univariate Polynomial Ring in x over Rational Field
The generator of the ring [Univariate Quotient Polynomial Ring in xbar over Rational Field with modulus x^2 + 1] is xbar
In K, xbar^2 = -1

# The lift command gives an element back in the original ring
# (the unique equivalent polynomial of degree less than the modulus)
print("With the lift command, we can recover the variable {} in {} (which maps to xbar in the residue ring)".format(xb.lift(),parent(xb.lift())))

With the lift command, we can recover the variable x in Univariate Polynomial Ring in x over Rational Field (which maps to xbar in the residue ring)

# Sage also supports interval arithmetic (RealIntervalField(n) = field with 100 bits of precision)
R = RealIntervalField(100)
a = R(sqrt(2))
print(a)
print(a.str(style='brackets'))
print(a.center())
print(a.precision())
print(a.endpoints())

1.414213562373095048801688724210?
[1.4142135623730950488016887242091 .. 1.4142135623730950488016887242108]
1.4142135623730950488016887242
100
(1.4142135623730950488016887242, 1.4142135623730950488016887243)

# Define an algebraic number, then get an interval of precision 100 bits containing the algebraic number
r1 = QQbar['x'](x^3-x+1).roots(multiplicities=False)[0] # Select the first root over the algebraic numbers
print(r1)
print(r1.interval(RealIntervalField(100)).str(style='brackets'))

-1.324717957244746?
[-1.3247179572447460259609088544787 .. -1.3247179572447460259609088544770]

# RIF = RealIntervalField(53)
print(RIF(sqrt(2),sqrt(3)).str(style='brackets')) # Defines interval
print(RIF)
si = sin(RIF(pi))
print(si.str(style='brackets'))
print(si.contains_zero())

[1.4142135623730949 .. 1.7320508075688775]
Real Interval Field with 53 bits of precision
[-3.2162452993532733e-16 .. 1.2246467991473533e-16]
True

# An alternative to interval arithmetic is ball arithmetic
# ComplexIntervalField(n) and ComplexBallField(n) are the analogous domains for complex numbers
RealBallField(100)(pi)

[3.14159265358979323846264338328 +/- 2.25e-30]

Part 4: Linear Algebra

# MatrixSpace defines the spaces of matrices
# Here we define some matrices and access their entries
MS = MatrixSpace(ZZ,3,4)
print(MS)
show(MS()) # Zero matrix
A = MS([1,2,3,4,5,6,7,8,9,10,11,12]) # Since the dimension is fixed, a matrix can be entered as a vector
show(A)

Full MatrixSpace of 3 by 4 dense matrices over Integer Ring

# A matrix can also be defined as a list of rows
show(MS([[1,2,3,4],[5,6,7,8],[9,10,11,12]]))

# Elements are accessed similar to arrays
show(A[1,2]) # Access single entry
show(A[-1,:]) # Acess a row (-1 = final row/column index and ":" = all rows/columns)
show(A[0:2,2:]) # Acess a submatrix (notation "k:" goes from index k to final index)
show(MS.change_ring(GF(3)).random_element())

# Matrix groups can also be defined
A = matrix(GF(11), 2, 2, [1,2,3,4])
B = matrix(GF(11), [[0,1],[1,0]])
MG = MatrixGroup([A,B])
show("The matrices ", A, "and ", B, " generate the matrix group ", MG, " over ", MG.base_ring());
if identity_matrix(GF(11),2) in MG:
    print("This group contains the identity")
else:
    print("This group does not contain the identity")

This group contains the identity

# Just like lists, matrices are stored as pointers -- need to use B = copy(A) to have independent copy
A = Matrix(ZZ,[[1,2],[3,4]])
B = A
C = copy(A)
A[0,0]=-1
show(A,"and ",B,"are both changed, but ",C," is not")

# One can also construct the image and kernel of a matrix
A = matrix(QQ,3,5,list(range(15)))
im = A.image()
ker = A.right_kernel()
show("We define the matrix A = ",A, " (over the rational numbers)")
show("The image of A is: ",im)
show("The right kernel of A is: ",ker)

# Same thing over a different ring
B = A.change_ring(GF(3))
show(B.image())
show(B.right_kernel())

# We can compute the eigenvalues of a matrix (over the appropriate field)
A = matrix(QQ,[[1,2],[3,4]])
B = matrix(GF(3),[[1,2],[3,4]])
print(A.eigenvalues())
print(B.eigenvalues())

[-0.3722813232690144?, 5.372281323269015?]
[1, 1]

# And get the eigenvector corresponding to each eigenvector
# (we get a list of the form (e,V,n) where e is the eigenvalue, V is a list of eigenvectors forming a
# basis for the corresponding right eigenspace, and n is the algebraic multiplicity of the eigenvalue)
print(A.eigenvectors_right())

[(-0.3722813232690144?, [(1, -0.6861406616345072?)], 1), (5.372281323269015?, [(1, 2.186140661634508?)], 1)]

# As noted above, algebraic numbers are stored exactly, even if only the decimal expansion is printed
eig1 = A.eigenvalues()[0]
print(eig1)
print(eig1.minpoly())

-0.3722813232690144?
x^2 - 5*x - 2

# Over appropriate rings, one can compute the Diagonalization, Jordan Form, etc. of a matrix
# A.jordan_form() # running this results in "RuntimeError: Some eigenvalue does not exist in Rational Field.""
R = QQ[eig1] # Define the extension of QQ by adding an eigenvalue of A
print(R)

A = A.change_ring(R)
M1,M2 = A.jordan_form(transformation=True, subdivide=False) # Now we can diagonalize (using the Jordan form command)

show(M1, M2)
show(M2,"times ", M1," times ",M2.inverse(),"equals A=",M2*M1*M2.inverse(), "where a is the algebraic number defined by ", R.gen().minpoly(), "=0, where x is approximately",  R.gen().n(20))

Number Field in a with defining polynomial x^2 - 5*x - 2 with a = -0.3722813232690144?

Part 5: Polynomial System Solving

# Reset all variables
reset()

# Multivariate polynomial rings are defined similarly to univariate rings
R.<x,y,z> = QQ['x,y,z']
p = x^2+y+z
show(p," is a polynomial in the ring ", parent(p))

# The order of variables matters when defining a ring
A = QQ['x,y']
B = QQ['y,x']
print(A==B)
print(A(x*y^2+y*x^2))
print(B(x*y^2+y*x^2))
print(A(x*y^2+y*x^2) == B(x*y^2+y*x^2))

False
x^2*y + x*y^2
y^2*x + y*x^2
True

#################################################################
# EXERCISE: Define a polynomial ring R over the rational numbers whose variables are xk for all primes k less than 100 (so the variables are x2,x3,x5,...)
# HINT: First make the string str = "x2,x3,x5,..." then use R = QQ[str]
#################################################################

# One can use an arbitrarily large number of variables, extending the number of variables as necessary
# (note that running this cell a second time will give a different behaviour than the first run)
R.<x,y> = InfinitePolynomialRing(ZZ, order='lex')
p = x[0]+y[0]-x[1]+3*y[1]
show(p)
show(parent(p))
show(parent(p.polynomial()))
show(parent((p + x[5]).polynomial()))
show(parent(p.polynomial()))

# Coefficients of a polynomial can be specified by a monomial or a list of exponents
R.<x,y> = QQ[]
p = (x+y)^2
print(p[x*y])
print(p[(1,1)])

2
2

# Ideals can be defined using the ideal method of a polynomial ring
R.<x,y,z> = QQ[]
J = R.ideal(x+y,x*z^2+y*x^3,z+x^2*y)

# The typical algebraic geometry methods can be applied to ideals and the corresponding varieties
# Can find solutions for zero-dimensional ideals
print(J.dimension())
print(J.variety()) # Solutions over the defining field QQ
print(J.variety(QQbar)) # Solutions over the algebraic closure QQbar

0
[{z: 0, y: 0, x: 0}, {z: 1, y: -1, x: 1}]
[{z: 0, y: 0, x: 0}, {z: 1, y: -1, x: 1}, {z: 1, y: 0.50000000000000000? - 0.866025403784439?*I, x: -0.50000000000000000? + 0.866025403784439?*I}, {z: 1, y: 0.50000000000000000? + 0.866025403784439?*I, x: -0.50000000000000000? - 0.866025403784439?*I}]

# To access elements of a variety, the user needs to define variables for poly ring over the correct field
(xx, yy, zz) = QQbar['x,y,z'].gens()
# print(x.subs(J.variety(QQbar)[-1])) # Gives error as "x" not variable for QQbar -- it is used above for QQ[x,y,z]
print(xx.subs(J.variety(QQbar)[-1]))
print([ pt[xx].degree() for pt in J.variety(QQbar) ])
print(set(pt[xx].minpoly() for pt in J.variety(QQbar)))

-0.50000000000000000? - 0.866025403784439?*I
[1, 1, 2, 2]
{x^2 + x + 1, x - 1, x}

# We can also define quotient rings
xbar = R.quo(J)(x)
print(parent(xbar))
print([xbar^k for k in range(5)])

Quotient of Multivariate Polynomial Ring in x, y, z over Rational Field by the ideal (x + y, x^3*y + x*z^2, x^2*y + z)
[1, -ybar, ybar^2, zbar, -ybar*zbar]

# q.lift(J) writes q as a polynomial combination of the elements of J
# Note we use R(1) to specify the constant 1 in R (instead of the integer 1, as is the default)
R(1).lift(ideal(1-x*y+y^2, x^3-y^2, x+2*y))

[-72/67*x^2 - 72/67*x*y - 48/67*y^2 - 18/67*x - 36/67*y + 1, -72/67*y + 9/67, 24/67*y^3 - 9/67*x^2 - 18/67*y^2 + 72/67*x - 5/67*y + 18/67]

# Can find radicals, primary decompositions, etc.
show(J.radical())
show(J.primary_decomposition())

# Can get elimination ideals
J.elimination_ideal(x)

Ideal (z^3 - z^2, y*z^2 - y*z, y^3 + z) of Multivariate Polynomial Ring in x, y, z over Rational Field

# Can reduce via repeated division algorithm (output depends on the order of the factors unless J is a Grobner Basis)
p = x^5-y^2+1
print(p.reduce(J))

y^2*z - y^2 + 1

# Can compute Grobner Bases
# Monomial term orders include lex, invlex, deglex, and degrevlex
GB = J.groebner_basis()
print(GB)

[y^3 + z, y*z^2 - y*z, z^3 - z^2, x + y]

Exercises

(Exercise 1)

An efficient method of computing high powers $x^N$ for an element $x$ of a ring is to use binary powering, through the formula $$\begin{equation} x^N = \begin{cases} \left(x^{N/2}\right)^2 &: \text{ $N$ even} \\ x \cdot \left(x^{(N-1)/2}\right)^2 &: \text{ $N$ odd} \end{cases} \end{equation}$$ Write a recursive function bin_pow(x,n) which takes an element $x$ from a ring (your implementation should not care which) and a natural number $n$ and returns $x^n$ in approximately $\log_2 n$ multiplications using binary powering.

(Exercise 2)

Recall again the Fibonacci numbers defined by the recursive formula $f_n = f_{n-1} + f_{n-2}$ for $n \geq 2$ and $f_0=f_1=1$. Show that the $N$th Fibonacci number $f_N$ can be recovered from the matrix product $$\begin{equation} \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^N \; \begin{pmatrix} 1 \\ 1 \end{pmatrix} \end{equation}$$ Using the bin_pow function from Exercise 1, find the one millionth Fibonacci number. If $M$ is the matrix being raised to the $N$th power, the command print(timeit('bin_pow(M,1000000)')) will measure the time it takes Sage to run this command. How does the time change if $M$ is defined over different rings (for instance, as a matrix over the integers, over the rationals, and over algebraic numbers)?

Note: Because Sage has built-in optimizations for powering a matrix, using your binary powering function may not be faster than running M^N. However, your function should only be slower by a small constant amount of time (not changing much, if at all, with $N$).

(Exercise 3)

Create the field extension of the rational numbers obtained by adjoining an eigenvalue of the matrix $M$ from Exercise 2. Working over this field extension, diagonalize $M$. Deduce an explicit expression for the $N$th Fibonacci number.

(Exercise 4)

The fastest method of determining terms in a linearly recurrent sequence, such as the Fibonacci numbers, is not to use matrix powers (although this is much faster than naively using the recurrence). Instead, the idea is to compute powers of elements in a residue class ring defined by such a recurrence.

Skipping the details of why this works, define the residue class ring of polynomials in $x$ with rational coefficients modulo the polynomial $x^2-x-1$. Let $\overline{x}$ be the image of $x$ in this ring. Use Sage to compute the powers $\overline{x},\overline{x}^2,\dots,\overline{x}^5$. Can you see how to recover the $N$th Fibonacci number from $\overline{x}^N$?

Use your bin_pow method from Exercise 1 to compute the one millionth Fibonacci number using this approach. Use the timeit command described above to compare the amount of time it takes to run bin_pow(xbar,1000000), bin_pow(M,1000000), and M^1000000, where xbar denotes $\overline{x}$ and $M$ is the matrix from Exercise 2. Which is fastest?

(Exercise 5)

A simple random walk of length $n$ on the lattice $\mathbb{Z}^2$ is a sequence of $n$ steps from the set $\{(-1,0),(1,0),(0,-1),(0,1)\}$ picked uniformly at random. Create a function which takes a natural number $n$ and returns a plot displaying a simple random walk of length $n$. Create a plot overlaying the result of 100 random walks of length 100, each of a distinct colour.

# The following functions will be useful.

# The choice function picks a random element of a list. For instance,
print([choice([1,2,3,4]) for k in range(10)])

# The line command plots a line from a list of points. For instance,
show(line([[0,0],[0,1],[2,2]], color = "black", thickness=2))

[1, 2, 3, 4, 3, 4, 2, 2, 4, 2]