CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutSign UpSign In

Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place. Commercial Alternative to JupyterHub.

| Download
Project: admcycles
Views: 68
Visibility: Unlisted (only visible to those who know the link)
Image: ubuntu2004
Kernel: SageMath 9.1

Lecture 3 : Numbers and symbolic expressions

References:

Summary:
In todays lecture, we'll start with a reminder of conditional execution (if) and loop (for) in Python, and the datatype of dictionaries. Then we discuss various types of numbers in SageMath, such as integers, rationals, algebraic numbers and more. In the last part, we discuss symbolic expressions and how to use them to solve equations.

Python warmup : conditional execution (if/else)

Sometimes we want to perform some operations in our code (like computing the inverse of a matrix AA) only if certain conditions are satisfied (the matrix AA is invertible), and maybe do something else otherwise (print a warning). This can be done in Python using the if/else construction:

A = matrix([[1,2],[3,5]]); A
if A.is_invertible(): print('The inverse of A is:') print(A.inverse()) else: print('A is not invertible!')

Here the general pattern is:

if CONDITION: do_something() do_more_of_something() else: do_something_else() maybe_more_of_something_else()

Note that the indentation (using a Tab = 4 spaces) is important!

Here CONDITION is an expressiong evaluating to a bool, which is either True or False. Basic ingredients for such conditions:

  • A == B, A != B for checking equality (or not-equality)

  • A <= B, A > B etc for comparisons of two values

  • Logical operations: Cond1 and Cond2, Cond1 or Cond2, not Cond1

It is possible to omit the else part, in which case nothing is done if CONDITION is not satisfied.

Exercise

Given real parameters p,qp,q of the quadratic equation x2+px+q=0 x^2 + p x + q = 0 write some code which prints the set of solutions xRx \in \mathbb{R} (which you have to compute yourself, since solving equations only comes later in the lecture).

Hint: It's never too late in life to look at the formula for the solutions of a general quadratic equation.

Solution: (uncomment to see)

Note that we can have an if-block inside another if- or else-block, using double indentation (and so on):

A = -6 if A.is_prime(): if A == 2: print('A is an even prime') else: print('A is an odd prime') else: if A < 0 and (-A).is_prime(): print('A is negative, but -A is a prime.') else: print('A is not a prime, since it factors as:') print(factor(A))

Python warmup: for-loops

Sometimes we want to perform an operation for all nn in a range of values. This can be accomplished by a for-loop. As an example, below we compute the sum of the numbers n=0,1,,100n=0,1, \ldots, 100:

count = 0 for n in range(101): count = count + n print(count)

The general pattern is:

for n in ITERABLE: do_something()

Above, we have used range(N) as the ITERABLE, which outputs the numbers 0,1,,N10,1, \ldots, N-1, which is why we had to take range(101) to sum from 00 to 100100.

Exercise

Count for how many of the numbers n=0,1,,2022n=0,1, \ldots, 2022 the value sin(n)\sin(n) is negative.

Solution (uncomment to see)

Similarly, we have the variations:

  • range(N0, N) for the numbers N0,N0+1,,N1N_0, N_0+1, \ldots, N-1,

  • range(N0, N, d) for the numbers N0,N0+d,N0+2d,,N0+kdN_0, N_0+d, N_0+2d, \ldots, N_0+k\cdot d, where kk is maximal such that N0+kd<NN_0+k\cdot d<N

Some more interesting examples:

for a in [2,4,'cat']: print(a)

Since we are not just using Python, but SageMath, we also have some cool iterables, like the symmetric group SnS_n given by SymmetricGroup(n):

for a in SymmetricGroup(3): print('The element {} of S_3 has order {}'.format(a, a.order()))

For more information about iterables, see this tutorial.

Python warmup: dictionaries

Dictionaries are a useful data structure in Python, which encodes a map from a finite set of objects to some other objects:

# The dictionary d sending 1 -> 2, 2 -> 2, 3 -> 'cat' d = {1:2, 2:2, 3: 'cat'}
d[1]

We can easily add new assignments to the dictionary (or change old ones):

d[4] = 3 d[3] = 'dog' d

There is also a connection to the for-loops we saw above: on the one hand, dictionaries are themselves Iterables:

for a in d: print(a,d[a])

On the other hand, similar to the constructions of lists like [k^2 for k in range(10)] that we already saw, we can also create dictionaries in this way:

s = {i : i+1 for i in range(10)} s

Exercise

  • Create dictionaries f,g implementing the maps f:{1,2,3}{5,6},15,26,35,g:{4,5,6}{7,8,9},49,58,69. f : \{1,2,3\} \to \{5,6\}, 1 \mapsto 5, 2 \mapsto 6, 3 \mapsto 5,\\ g : \{4,5,6\} \to \{7,8,9\}, 4 \mapsto 9, 5 \mapsto 8, 6 \mapsto 9.

  • Create the dictionary h implementing the composition gfg \circ f.

Solution (uncomment to see)

Numbers

Last time, we already saw some interesting rings implemented in SageMath:

  • ZZ the integers

  • QQ the rational numbers

  • RR the real numbers (floating point operations!)

  • CC the complex numbers (floating point operations!)

  • GF(q) the finite field Fq\mathbb{F}_q with qq elements

In this section, we will learn more about elements of these rings (commonly known as "numbers") and what we can do with them.

Integers

When we type an expression like 42 in SageMath, this will produce an Integer:

type(42)

Integer division and remainder

Apart from the usual division a/b of two Integers (which mostly produces a Rational), we can also use a//b to compute the largest integer which is at most a/ba/b:

5/3
5//3
6//3

The remainder (or modulus) of this division can be computed with a%b:

5 % 3

Exercise

Check whether the number 1717117^{17}-1 is divisible by 33.

Solution (uncomment to see)

Exercise

Compute the last 33 digits of the integer a=(((22)2))2 a = (((2^2)^2)^\ldots)^2 where the number 22 appears

  • 55 times, or

  • 500500 times

in the formula for aa.

Hint 1: Since the number of digits of an integer roughly doubles when taking its square, we can estimate that aa has more decimal digits than the number of atoms in the universe. Keep in mind that the code you use for solving the exercise must run on a computer which is part of the universe!
Hint 2: If you have ignored Hint 1, you probably at some point want to click on Kernel -> Interrupt in the menu to stop the computation, and maybe also on Kernel -> Restart to clear your memory.

Solution (uncomment to see)

Factorization

We can also compute the prime factorization of an integer:

u = factor(99); u
type(u)

To actually extract the information from this factorization, to use it in further computations, the following conversions are useful:

list(u)
dict(u)

Exercise

Find the largest m1m \geq 1 such that the number 100!100! is divisible by 2m2^m.

Solution (uncomment to see)

Conversely, we can look up the primes in the interval [a,b1][a,b-1] using primes(a,b):

primes(4,19)

Oops, this function returns a generator, i.e. an iterable that can be used in a for loop, without having to first compute all the numbers it will put out:

for p in primes(1,10^100): print(p) sleep(1)

To just look at the numbers, you can convert this generator object into a list:

list(primes(4,19))

Python int vs. SageMath Integer

Remember the exercise from the first lecture to create a multiplication table of the integers from 11 to 1010:

matrix([[a*b for a in range(1,11)] for b in range(1,11)])

Let's see what happens instead if we wanted to create a division table:

matrix([[a/b for a in range(1,11)] for b in range(1,11)])

By golly, what happened? It turns out that the integers that range spits out are actually of type int in Python, and not of type Integer:

for i in range(10): print(type(i))

Therefore, the division above is a division of int values, which returns a float:

a = int(3)/int(7); a
type(a)

To repair this, you can convert either numerator or denominator of the fraction into an Integer:

b = ZZ(int(3))/int(7); b
type(b)

Thus a nice and mathematically accurate division table can be obtained as follows:

matrix([[ZZ(a)/b for a in range(1,11)] for b in range(1,11)])

This type of conversion works more generally:

c = 6/2 (c,type(c))
d = ZZ(c) (d,type(d))

Exercise

An integer aa is a perfect power if there exist integers b,eb,e with e>1e>1 such that a=bea=b^e.

[a.is_perfect_power() for a in [3, 9, 6]]

For the integers 1,,201, \ldots, 20 print line by line whether they are perfect powers or not, e.g. the output should start like this:

1 True 2 False 3 False 4 True 5 False ...

Solution (uncomment to see)

Rational numbers

If you want to do an exact computation with real numbers (e.g. in linear algebra), it is best when you can stay inside the rational numbers, since they have an exact representation as a fraction. Many operations we saw for Integers above work analogously for Rationals:

a = 24/9; a
print(a.numerator(), a.denominator())
b = QQ(1.45); b
factor(b)

Some more useful functions (which also work for integers or real numbers):

a = - 5/3 floor(a) # rounding down
ceil(a) # rounding up
sign(a) # the sign
abs(a) # absolute value

Exercise

Among the numbers of the form p/qp/q for p,q{1,,20}p,q \in \{1, \ldots, 20\}, which is closest to π\pi?
Hint: Have a variable storing the best approximation you have so far, and when going through all numbers p/qp/q replace the current optimum if one of them comes closer.

Solution (uncomment to see)

Real and complex numbers

We already saw that RR and CC are the fields of real and complex numbers, implemented with floating-point arithmetic. In particular, the computations here are no longer exact:

piRR = RR(pi); piRR
sin(2*piRR)
a = CC(- pi * I); a
exp(a)

Since we mostly talk about symbolic computations in this course, we won't go deeper into this topic for now, but here are some links where you can find more information:

Algebraic numbers

One nice field in which exact computations are possible is the field Q\overline{\mathbb{Q}} of algebraic numbers. These are complex numbers which are a zero of some polynomial with rational coefficients. The unique such polynomial of minimal degree and with leading coefficient 11 is called the minimal polynomial.

z = QQbar(3/2 + I); z
z.minpoly()
z^2 - 3*z + 13/4

Then we can do things like taking roots of an algebraic number.

u = sqrt(z); u

Even though below it is displayed as a decimal number, behind the scenes it still remembers its exact representation:

u.minpoly() # for some reason this produces an error on the SageMath installation at the University

Then we can do fun things, like compute an expansion of a real algebraic number as a continued fraction:

w = sqrt(QQbar(2)); w
continued_fraction(w)

This says that 2=1+12+12+12+ \sqrt{2} = 1 + \frac{1}{2 + \frac{1}{2 + \frac{1}{2+ \ldots} }} If you haven't seen it, it's a fun exercise to prove this!

Exercise

Compute the minimal polynomial and continued fraction of the golden ratio ϕ=1+52. \phi = \frac{1 + \sqrt{5}}{2}\,.

Solution (uncomment to see)

Finite fields

Given a prime number pp, there exists the finite field Fp=Z/pZ\mathbb{F}_p = \mathbb{Z}/p\mathbb{Z} of integers modulo pp, which has precisely pp elements. More generally, for a prime power q=prq=p^r, (r1)(r \geq 1), there exists a unique field Fq\mathbb{F}_q with qq elements. It is not given by Z/qZ\mathbb{Z}/q \mathbb{Z} though! Instead it is given as an extension Fq=Fp[z]/(f(z)), \mathbb{F}_q = \mathbb{F}_p[z]/(f(z)), where zz is a variable and fFp[z]f \in \mathbb{F}_p[z] is an irreducible polynomial of degree rr.

This finite field can be created in SageMath using FiniteField(q) (or GF(q) for short).

GF(9)

If we actually need to write elements of Fq\mathbb{F}_q, we can use the following line to specify the variable used in the definition:

R.<z> = GF(9) R

We can check e.g. that z+z+z=0z+z+z=0, since z+z+z=3z=0z=0z+z+z = 3 \cdot z = 0 \cdot z = 0, as 3=0F323 = 0 \in \mathbb{F}_{3^2}.

print(z+z) print(z+z+z)

We can also type formulas using standard integers like 17. Of course by default they give elements of type Integer, but SageMath is smart enough to convert them to elements of FpFq\mathbb{F}_p \subseteq \mathbb{F}_q. We are later going to see more about how this works.

u = z + 17; u
u+1

And we can check which irreducible polynomial ff was used to define Fq=Fp[z]/(f(z))\mathbb{F}_q = \mathbb{F}_p[z]/(f(z)); SageMath automatically picks one for us:

z.minimal_polynomial()

Here is a cool exercise using finite fields, which hints at some quite deep mathematics:

Exercise

Given a prime power qq, let Lq(f)\mathbb{L}_q(f) be the set of solutions (x,y,z)Fq3(x,y,z) \in \mathbb{F}_q^3 of the equation f=x3y3x2z2yz+4=!0Fq. f = x^3 y - 3x^2 z^2 - yz + 4 \stackrel{\text{!}}{=} 0 \in \mathbb{F}_q.

  • Compute the cardinality L2(f)|\mathbb{L}_2(f)| by going through all points in F23\mathbb{F}_2^3 and counting how many of them satisfy the equation.

  • Compute the cardinality Lq(f)|\mathbb{L}_q(f)| for a few more values of qq and then make a guess how it grows with qq.
    Hint: There is a function in SageMath giving you a list of all prime powers qq up to n1n-1. Can you guess its name and confirm with Tab-completion?

Solution: (uncomment to see)

Symbolic expressions

We have already seen in some examples that SageMath can do computations involving formal variables. These happen in the so-called Symbolic Ring, which is denoted by SR in SageMath.

Creating and comparing symbolic expressions

To get started, we declare some variable names using the function var.

var('a,b') f = (a+b)^3; f
type(f)
f in SR

Then we can ask to perform operations, like expanding the formula we entered above:

expand(f)

We can also try to check whether two formulas are equivalent:

t1 = (a^2+2*a*b+b^2)/(a+b); t1
t2 = a+b; t2
t1 == t2
bool(t1 == t2)

This is not restricted to algebraic equations:

bool(sin(a)^2 + cos(a)^2 == 1)

Exercise

Try to check the following identities, some of which are true and some of which are false. Which one does SageMath get wrong? amax(a,b)exp(a)exp(b)=exp(a+b)11a1b=abab12(a+b)max(a,b)(a+b+c)2=a2+b2+c2+2(ab+ac+bc)sin(π4)=12 a \leq \max(a,b) \\ \exp(a) \cdot \exp(b) = \exp(a+b) \\ \frac{1}{\frac{1}{a}-\frac{1}{b}} = \frac{a \cdot b}{a - b} \\ \frac{1}{2}(a+b) \leq \max(a,b) \\ (a+b+c)^2 = a^2 + b^2 + c^2 + 2 \cdot(ab + ac + bc) \\ \sin\left(\frac{\pi}{4}\right) = \frac{1}{\sqrt{2}}

Solution (uncomment to see)

Substitution

We can make variable substitutions as follows:

var('a,b,c,d') f = a^2 + b*c; f
f.subs({a:b, c:b+1})

Note that the argument of f.subs is a dictionary, sending variables appearing in f to their new values. Alternatively, we can treat f as if it was a function, and call it (but specifying the names of the variables as follows):

f(a=2,b=3,c=4)
f(a=d,b=d,c=d)

Formal sums

We can also evaluate some formal summations, e.g. the famous identity 1+2++n=n(n+1)2 1 +2 + \ldots + n = \frac{n(n+1)}{2} Indeed, to compute a=a0a1f(a) \sum_{a=a_0}^{a_1} f(a) we use sum(f,a,a0,a1), where a is a variable and f is again a symbolic expression.

var('n') g = sum(a,a,1,n); g

The answer is again a symbolic expression, so we can e.g. substitute some value for n or compute a further symbolic sum:

g.subs({n:100}) # The famous sum that Gauss allegedly computed in elementary school
var('m') sum(g,n,1,m)

Solving equations

Maybe you have wondered why writing f == g for two symbolic expressions does not immediately give either True or False, but again a symbolic expression:

g == 5050
type(g == 5050)

The reason is that we can now use this expression to state an equation (or a system of equations) that we want to solve:

solve(g == 5050, n)

The second argument is the variable for which we want to solve. To specify multiple equations (or inequalities) we give a list of these as the first argument of solve:

solve([g == 5050, n>=0], n)

To get an output that can be used in further computations, we can specify the optional parameter solution_dict=True to obtain a list of dictionaries giving our solutions:

solve([g == 5050, n>=0], n, solution_dict=True)
var('a,b') L = solve(a^2 - b^2 == 0, b, solution_dict=True); L
for l in L: print((a^2 - b^2).subs(l))

Finally, we can also solve for multiple variables:

var('y') solve([x*y==0, (x-1)*(y-1)==0],[x,y])

Exercise

Solve the following (systems of) equalities for the variables x,y,n,x,y,n,\ldots: x3+2x221x+18=0,x2+px+q=0,k=0n2k=4095x2+y2=1,x+y=1 x^3 + 2 x^2 - 21 x + 18 = 0, \\ x^2 + p x + q = 0,\\ \sum_{k=0}^n 2^k = 4095\\ x^2 + y^2 = 1, x + y = 1

Solution (uncomment to see)

We can also add some assumptions about the variables:

solve(x^3==1, x, solution_dict=True)
assume(x, 'real') solve(x^3==1, x, solution_dict=True)

Here is how we get rid of the assumptions again:

forget(assumptions()) solve(x^3==1, x, solution_dict=True)

Assignments

Exercise

Given a number d{2,3,,9}d \in \{2,3, \ldots, 9\} here are the rules of a game named dd-Blup:

The participants stand in a circle and count upwards starting from 11, except that for every number either divisible by dd or ending on the digit dd, they have to say "Blup" instead. If you say something wrong or take too long, you are out, and the others start again.

A correct 33-Blup would start as follows:

Player 1 : "1" Player 2 : "2" Player 3 : "Blup" Player 4 : "4" ... Player 11 : "11" Player 12 : "Blup" Player 13 : "Blup" ...

Given d, write a program that prints the first 50 things the players need to say, i.e.

1 2 Blup 4 ...

Hint: Recall that given integers a,b with b nonzero, you can compute the remainder of the division of a by b using a % b.

Solution: (uncomment to see)

Exercise

The Collatz conjecture states that for the function $$ f : \mathbb{Z}_{\geq 1} \to \mathbb{Z}_{\geq 1}, n \mapsto \begin{cases} 3n+1 & \text{for $nParseError: KaTeX parse error: Expected 'EOF', got '}' at position 5: odd}̲,\\ n/2 & \text…nParseError: KaTeX parse error: Expected 'EOF', got '}' at position 6: even}̲ \end{cases} andanypositiveinteger and any positive integer n,thesequence, the sequence n,, f(n),, f(f(n)),, \ldots$ eventually hits the number 11. Test the Collatz conjecture for the first 100100 integers.
Hint: The while-loop

while CONDITION: do_something()

repeats the following code block until the CONDITION becomes False.

Solution (uncomment to see)

Exercise

The Pólya conjecture states that given a natural number n>1n>1, there are among the set {1,,n}\{1, \ldots, n\} at least as many elements with an odd number of prime factors as elements with an even number of prime factors. Here prime factors are counted with multiplicity, i.e. 12=22312 = 2 \cdot 2 \cdot 3 has 33 prime factors according to this definition.

  • Test the Pólya conjecture for the integers 2nn02 \leq n \leq n_0, for an n0n_0 of your choice.

  • Google the Pólya conjecture.

Solution: (uncomment to see)

Exercise

Given an integer c0c \geq 0 let TcZ2T_c \subseteq \mathbb{Z}^2 be the set of points in Z2\mathbb{Z}^2 contained in the triangle spanned by (0,0),(0,c),(c,0)(0,0), (0,c), (c,0). Or in other words Tc={(x,y)Z2:x,y0,x+yc}. T_c = \{(x,y) \in \mathbb{Z}^2 : x,y \geq 0, x+y \leq c\}. We claim that the function g:Z0Z,c(x,y)Tcxy g : \mathbb{Z}_{\geq 0} \to \mathbb{Z}, c \mapsto \sum_{(x,y) \in T_c} x \cdot y is given by a polynomial in cc. Compute this polynomial and check your result in a few cases.

Solution (uncomment to see)