3. Number Theory#
3.1. What is Number Theory?#
Number Theory, to put it very simply, is the study of integers, i.e., the set
It is one of the oldest branches of mathematics. In particular, prime numbers (as we will discuss below) have a central role in Number Theory.
For a long time it was an area of with little concrete applications. G. H. Hardy, a very famous number theorist, once called Number Theory “useless”.
I have never done anything ‘useful’. No discovery of mine has made, or is likely to make, directly or indirectly, for good or ill, the least difference to the amenity of the world.
(Emphasis mine.)
But with the advent of computers, Number Theory has found many applications, and one the most notable is cryptography.
We will then spend sometime learning some of the basics of Number Theory (with a very concrete and computational approach) before we start with the applications.
3.2. Bibliography#
Although we will introduce all the necessary topics in number theory, there are many good resources for the reader interested in learning more about this beautiful and important subject. Here are some good introductory texts:
William Stein’s Elementary Number Theory: besides providing an excellent introduction to the theory, topics are often illustrated with computations in Sage. (Stein was the creator of Sage!) It makes a very good “companion” to this text. (It is also freely available.)
Niven, Zuckerman, and Montgomery’s An Introduction to the Theory of Numbers: a classic and also a great introduction.
Kenneth Rosen’s Elementary Number Theory and It’s Applications: another good introduction. It also contains applications to cryptography.
Kenneth Ireland and Michael Rosen’s A Classical Introduction to Modern Number Theory: a more advanced book, but it does start with the basics. It is very well written and provided many good exercises.
3.3. Long Division#
You probably remember long division from school. You can break an integer as a multiple of another plus some remainder smaller than what you are dividing by. More precisely:
Theorem 3.1 (Long Division)
Given \(a,b \in \mathbb{Z}\), with \(b \neq 0\), there is a unique pair of integers \((q,r)\) such that \(a=bq+r\) and \(r \in \{0, 1, 2, \ldots, |b|-1\}\).
Note
Remember that the symbol \(\in\) means “belongs to”. So \(a, b \in \mathbb{Z}\) is saying that \(a\) and \(b\) belong to the set of integers, i.e., \(a\) and \(b\) are integers.
Note that \(|b|\) above is the absolute value of \(b\), meaning that if \(b\) is non-negative, then \(|b|=b\), and if it is negative, then \(|b| = -b\), i.e., the “positive version” of \(b\).)
For positive integers, the process is fairly simple: we start subtracting \(b\) from \(a\) until what we have is less than \(b\). The number of times we subtract \(b\) from \(a\) is the quotient \(q\) and what is left (the number less than \(b\), perhaps \(0\)) is the remainder \(r\).
So, \(b\) goes into \(a\) exactly \(q\) times, with a left-over of \(r\).
Example (Long Division)
What are the quotient and remainder of \(183\) when divided by \(37\)?
a, b = 182, 37
a - b
145
a - 2*b
108
a - 3*b
71
a - 4*b
34
Ah-ha! Since \(34\) is less than \(37\), we have that the remainder is \(34\) and the quotient is \(4\). So, we have
Let’s automate this process:
def positive_long_div(a, b):
"""
Given two positive integers a and b, return the quotient and remainder
of the long division of a by b.
INPUTS:
a: the dividend (positive integer)
b: the divisor (positive integer)
OUTPUT:
(q, r): q the quotient and r the remainder of the long division of a by b.
"""
q = 0 # initialize the quotient as 0
r = a # initialize the remainder as a
while r >= b:
q += 1
r -= b
return (q, r)
Let’s test it!
positive_long_div(182, 37)
(4, 34)
If the dividend \(a\) is already less than the divisor \(b\), the quotient should be \(0\) and the remainder \(a\):
positive_long_div(5, 7)
(0, 5)
We should also be able to do the same when the dividend \(a\) is negative. Note that we still want the remainder to be between \(0\) and \(b-1\) (inclusive in both ends). In this case, instead of subtracting, we need to add \(b\) until we get a non-negative number.
So, if we were to divide \(-182\) by \(37\):
a, b = -182, 37
a + b
-145
a + 2*b
-108
a + 3*b
-71
a + 5*b
3
Since \(3\) is our first non-negative number, it is the remainder, and \(-5\) (note the negative) is the quotient:
In principle \(b\) could also be negative, but let’s not bother with that case, as it is rather unusual. So, lets automate the process for negative and positive dividends in a single function:
def long_div(a, b):
"""
Given an integer a and a positive integer b, return the quotient and remainder
of the long division of a by b.
INPUTS:
a: the dividend (integer)
b: the divisor (positive integer)
OUTPUT:
(q, r): q the quotient and r the remainder of the long division of a by b.
"""
if a >= 0:
q = 0 # initialize the quotient as 0
r = a # initialize the remainder as a
while r >= b:
q += 1
r -= b
return (q, r)
# no need for else, as we only get here when a < 0:
q = -1 # initialize the quotient as -1
r = a + b # initialize the remainder as a + b
while r < 0:
q -= 1
r += b
return (q, r)
Testing:
long_div(182, 37)
(4, 34)
long_div(-182, 37)
(-5, 3)
The work above is a good exercise in thinking algorithmically and writing Sage/Python code, but, as one would expect, we could also have done this directly, as // gives the quotient and % gives the remainder!
182 // 37, 182 % 37
(4, 34)
-182 // 37, -182 % 37
(-5, 3)
In fact, Sage (but not Python) has the function divmod:
help(divmod)
Help on built-in function divmod in module builtins:
divmod(x, y, /)
Return the tuple (x//y, x%y). Invariant: div*y + mod == x.
So, to make sure our code is OK, let’s test our function against Sage’s:
number_of_tests = 1_000
max_int = 1_000_000
for _ in range(number_of_tests):
a = randint(-max_int, max_int)
b = randint(1, max_int) # it cannot be zero!
result = long_div(a, b) # my result
expected = divmod(a, b) # Sage's
if result != expected: # check for problems and print data!
print(f"Failed for {a = } and {b = }.")
print(f"Sage's result: {expected}.")
print(f"My result: {result}")
break # stop!
else: # runs when for loop ended without break!!!!!
print("It worked (for all tests)!")
It worked (for all tests)!
Important
Note that for-loops in Python/Sage can have an else statement. It runs when the loop runs reaches its end without being manually stopped by the break keyword. (A better name would have been nobreak instead of else.)
Definition 3.1 (Divisibility)
We say that \(b\) divides \(a\) if the remainder of the long division of \(a\) by \(b\) is zero, or, equivalently, if there is an integer \(q\) such that \(a = bq\).
If “\(b\) divides \(a\)”, we can also say:
\(b\) is a factor of \(a\), or
\(a\) is a multiple of \(b\).
(They all mean the same thing.)
Notation: We write \(b \mid a\) for “\(b\) divides \(a\)”.
Attention
Note that \(b / a\) is a number, meaning \(b\) divided by \(a\). On the other hand \(b \mid a\) is a boolean (i.e. True or False), meaning \(b\) divides \(a\).
If we only want to test for divisibility, we can simply use a % b in Sage/Python. If it is zero, then b does divide a, and if not zero, then b does not divide a.
3.4. Primes#
Definition (Prime and Composite)
A positive integer \(p\) is prime if it is not \(1\) and the only positive divisors are \(1\) and \(p\) itself. Integers greater than \(1\) that are not prime are called composite.
Note
Note that \(1\) is neither prime nor composite!
Examples are \(2\), \(3\), \(5\), \(7\), and \(11\), the first five primes. The number \(4\) is composite, since \(2\) is a positive divisor different from \(1\) and \(4\).
Testing if a positive integer is prime is a very important problem in Number Theory. Let’s write a very naive function to test if a number is prime, by looking for all possible divisors:
def prime_test(n):
"""
Given an integer n >= 1, returns True if n is prime and False
if it is not.
INPUT:
n: number to be test for primality (integer >= 1).
OUTPUT:
boolean for whether or not n is prime.
"""
# we do n = 2 separate
if n == 2:
return True
# now, test for divisors
for div in xsrange(2, n): # stops at n-1
if n % div == 0: # found divisor
return False
# if we get here, no divisor was found and the number is prime!
return True
Let’s test for some small numbers that we know to be prime:
for n in xsrange(2, 20):
if prime_test(n):
print(f"{n = } is prime!")
else:
print(f"{n = } is composite.")
n = 2 is prime!
n = 3 is prime!
n = 4 is composite.
n = 5 is prime!
n = 6 is composite.
n = 7 is prime!
n = 8 is composite.
n = 9 is composite.
n = 10 is composite.
n = 11 is prime!
n = 12 is composite.
n = 13 is prime!
n = 14 is composite.
n = 15 is composite.
n = 16 is composite.
n = 17 is prime!
n = 18 is composite.
n = 19 is prime!
It seems to work, but it is relatively very slow. Before we time it, here is some information about the computer running these computations:
!inxi --system --cpu --memory
System:
Host: dell7010 Kernel: 6.17.8-1-siduction-amd64 arch: x86_64 bits: 64
Desktop: KDE Plasma v: 6.5.3 Distro: siduction 22.1.2 Masters_of_War -
kde - (202303220911)
Memory:
System RAM: total: 128 GiB available: 125.51 GiB used: 29.15 GiB (23.2%)
Array-1: capacity: 128 GiB slots: 4 modules: 4 EC: None
Device-1: DIMM1 type: DDR5 size: 32 GiB speed: spec: 4800 MT/s
actual: 3600 MT/s
Device-2: DIMM2 type: DDR5 size: 32 GiB speed: spec: 4800 MT/s
actual: 3600 MT/s
Device-3: DIMM3 type: DDR5 size: 32 GiB speed: spec: 4800 MT/s
actual: 3600 MT/s
Device-4: DIMM4 type: DDR5 size: 32 GiB speed: spec: 4800 MT/s
actual: 3600 MT/s
CPU:
Info: 24-core model: 13th Gen Intel Core i9-13900 bits: 64 type: MCP cache:
L2: 32 MiB
Speed (MHz): avg: 1011 min/max: 800/5300:5600:4200 cores: 1: 1011 2: 1011
3: 1011 4: 1011 5: 1011 6: 1011 7: 1011 8: 1011 9: 1011 10: 1011 11: 1011
12: 1011 13: 1011 14: 1011 15: 1011 16: 1011 17: 1011 18: 1011 19: 1011
20: 1011 21: 1011 22: 1011 23: 1011 24: 1011
Now, for instance, \(30{,}000{,}001\) is prime, but it takes us a little while to get there:
%%time
prime_test(3 * 10^7 + 1)
CPU times: user 2.19 s, sys: 1.52 ms, total: 2.2 s
Wall time: 2.2 s
True
But here is a simple mathematical idea that can speed up this quite a bit!
If \(n = a \cdot b\), with \(a\) and \(b\) integers greater than \(1\) (and hence also smaller than \(n\)), what can we say about the size of the smallest factor, say \(a\)? (Note that we allow \(a = b\).)
If \(a > \sqrt{n}\), then also \(b > \sqrt{n}\), and so
i.e., \(n > n\), which is impossible! Therefore, we must have that, if \(n\) is composite, then the smallest divisor must be less than or equal to \(\sqrt{n}\)!
Therefore, we only need to test for divisors up to \(\sqrt{n}\). If none is found up to there, then \(n\) must be prime.
So, in our original function, the worst case scenario was that we performed \(n-2\) divisions (looking for factors), but with this idea we perform at most \(\sqrt{n} - 2\) divisions.
That is a huge difference for large numbers. For instance, for \(n = 30{,}000{,}001\), we had to perform \(20{,}999{,}999\) divisions, since it was prime. With this idea, we perform only \(5{,}475\)!
Let’s implement and compare!
def prime_test_2(n):
"""
Given an integer n >= 1, returns True if n is prime and False
if it is not.
INPUT:
n: number to be test for primality (integer >= 1).
OUTPUT:
boolean for whether or not n is prime.
"""
# we do n = 2 separate
if n == 2:
return True
# now, test for divisors
for div in xsrange(2, floor(sqrt(n)) + 1): # stops at sqrt(n) (NOTE the "+1")
if n % div == 0: # found divisor
return False
# if we get here, no divisor was found and the number is prime!
return True
%%time
prime_test_2(3 * 10^7 + 1)
CPU times: user 2.5 ms, sys: 16 µs, total: 2.51 ms
Wall time: 2.46 ms
True
Important
This very simple idea gave a huge improvement, but required some mathematical knowledge (and maturity). That is often the case: no matter how good a coder you are, or what language and computing tricks you use, knowing some math can have a much greater impact in your code!
3.5. The Sieve of Eratosthenes#
Let’s see if we can improve the previous algorithm. A very natural idea, which we would likely do when checking for primality by hand, is to only check divisibility by \(2\) and then skip all other even numbers. After all if \(2\) does not divide a number, no other even number will. This would cut the time (in worst case scenarios, i.e., when the number is prime) in half!
But them, similarly, after we check if the number is divisible by \(3\), we do not need to check any other multiples of \(3\), further cutting the time in one third!
And then we do the same for \(5\). (No need for \(4\), since it is even.) Then \(7\). (Six is also even.) Then \(11\). (Eight and \(10\) are even and \(9\) is divisible by \(3\).) And so on…
We then see that we only need to check divisibility by primes!
But wait, we don’t know which numbers are prime yet. (That’s exactly what we are trying to do here.)
The Sieve of Eratosthenes offers a solution that uses this natural idea! It is a bit more memory intensive, since we have store (more or less) all numbers up to the one we want to test in memory. But, at the end, it does not only check if the number itself is prime: it gives us all primes up to that number!
Let’s see how it works.
Suppose we want to know if \(101\) is prime. We start with a list of all integer between \(2\) and \(101\) both inclusive:
We know that \(2\) is prime and we can then cross all multiples of \(2\) from the list, which we can do without division, by just crossing every other element:
Then, we are sure that the next non-crossed out element is also prime, in this case, \(3\). And again, we can cross out all its multiples, by jumping with a step size of \(3\):
Now we repeat! The next non-crossed out element, namely \(5\) is prime, and we cross out all its multiples:
Then repeat for \(7\):
But now, the next non-crossed out number is \(11\) which is larger than \(\sqrt{101}\) and we can stop! Every non-crossed out element is actually prime:
So, the list of primes less than or equal to \(101\) is: \(2\), \(3\), \(5\), \(7\), \(11\), \(13\), \(17\), … Let’s check with Sage’s prime_range:
prime_range(102)
[2,
3,
5,
7,
11,
13,
17,
19,
23,
29,
31,
37,
41,
43,
47,
53,
59,
61,
67,
71,
73,
79,
83,
89,
97,
101]
Note that if our interest is only whether or not a number is prime (and not finding all primes up to that number), we can stop if it is crossed out, as then we know that it is not prime.
Now, let’s implement it:
def prime_sieve(maxn):
"""
Given an integer maxn, returns a list of all primes up to (and possibly including) maxn.
INPUT:
maxn: the uppder bound for our list of primes.
OUTPUT:
A list of all primes up to (and possibly including) maxn.
"""
# dictionary for primes: key = number, value: is it prime?
prime_dict = {i: True for i in range(2, maxn + 1)}
last_to_try = floor(sqrt(maxn)) # we can stop here
p = 2 # we start with 2
while p <= last_to_try:
# mark multiples of p as not prime
for j in range(2 * p, maxn + 1, p):
prime_dict[j] = False
# move on to the next prime
p += 1
while (p <= last_to_try) and (not prime_dict[p]):
p += 1
return [x for x in prime_dict if prime_dict[x]]
Let’s test it:
prime_sieve(101) == prime_range(102)
True
To be sure, let’s test for a large integer:
prime_sieve(10^6) == prime_range(10^6 + 1)
True
Let’s see how efficient it is:
%%time
n = 3*10^7 + 1
n in prime_sieve(n)
CPU times: user 5.49 s, sys: 376 ms, total: 5.86 s
Wall time: 5.88 s
True
A lot slower, but gives a list of primes. Let’s see how long our previous method would take to give the same list:
%%time
n = 10^5
test_1 = [p for p in range(2, n + 1) if prime_test_2(p)]
CPU times: user 2.88 s, sys: 67 µs, total: 2.88 s
Wall time: 2.89 s
%%time
test_2 = prime_sieve(n)
CPU times: user 9.72 ms, sys: 0 ns, total: 9.72 ms
Wall time: 9.75 ms
So, for lists of primes, it is a lot faster!
Let’s check that the results match:
test_1 == test_2
True
3.6. The Greatest Common Divisor#
Given two integers \(a\), \(b\), the greatest common divisor (GCD) of \(a\) and \(b\), denoted by either \(\gcd(a, b)\) or simply \((a, b)\) (when it is clear that is not simply an ordered pair) is the greatest positive integer that divides both \(a\) and \(b\).
Property 3.1 (Properties of the GCD)
Since \(1\) divides any integer \(a\) (since \(a = 1 \cdot a\), and \(a\) is an integer), we have that \(\boxed{\gcd(1, a) = 1}\) for any integer \(a\).
Any integer \(a\) divides \(0\), since \(0 = a \cdot 0\) (and \(0\) is an integer). So, if \(a\) is a positive integer, then \(\boxed{\gcd(a, 0)=a}\).
Signs do not affect divisibility: if \(a \mid b\), then \(-a \mid b\), \(a \mid -b\) and \(-a \mid -b\), and so \(\boxed{\gcd(a, b) = \gcd(|a|, |b|)}\).
If \(a\) and \(b\) are positive, and \(a \mid b\), then \(\gcd(a, b) = a\).
If \(p\) is prime, then
\[\begin{split} \gcd(p, a) = \begin{cases} 1,& \text{if $p \nmid a$,}\\ p,& \text{if $p \mid a$.} \end{cases}\end{split}\]
Definition (Relatively Prime)
We say that two integers are relatively prime if their GCD is \(1\).
The following result about GCD is quite useful, but not as simple as the ones above. One can prove it using Bezout's Lemma or Prime Factorization below.
Proposition 3.1
Suppose that \(a \mid c\) and \(b \mid c\). If \(a\) and \(b\) are relatively prime, i.e., if \(\gcd(a, b)=1\), then we also have that \(ab \mid c\).
This is certainly not necessarily the case if \(a\) and \(b\) are not relatively prime. For instance \(6\) and \(9\) both divide \(18\), but \(6 \cdot 9 = 36\) clearly does not.
3.6.1. Computing the GCD#
A first naive method to compute the GCD is to find the divisors of both numbers and take the largest common one. For instance, say we want \(\gcd(12, 15)\). We have
But the process of finding all these divisors might take a while. Here is a second (still naive) approach. We can simply start testing:
Does \(1\) divide both \(12\) and \(15\)? Yes! So, so far the GCD is \(1\).
Does \(2\) divide both \(12\) and \(15\)? No. So, so far the GCD is still \(1\).
Does \(3\) divide both \(12\) and \(15\)? Yes! So, so far the GCD is \(3\).
Does \(4\) divide both \(12\) and \(15\)? No. So, so far the GCD is still \(3\).
When do we stop? We stop when we reach the smallest of the two numbers, namely, \(12\). The last common divisor will be the GCD. (In this case, it is \(3\).)
But we can greatly improve it! Since we are interested in the largest common divisor, we start with the largest it could possibly be, namely, the smallest of the two numbers (\(12\) in our case). So, we do:
Does \(12\) divide both \(12\) and \(15\)? No, so we go to the next number.
Does \(11\) divide both \(12\) and \(15\)? No, so we go to the next number.
Does \(10\) divide both \(12\) and \(15\)? No, so we go to the next number.
When do we stop? When we find a common divisor! This will automatically be the GCD. So, in this case, we perform \(18\) divisions instead of \(22\). In some other cases, we can perform a single one. For instance, in computing \(\gcd(1{,}000{,}000, 2{,}000{,}000)\) we perform a single division, instead of \(999{,}999\)!
But, there is a much more efficient way of doing it!
3.7. The Euclidean Algorithm#
We have the following:
Theorem 3.2
We have that \(d\) divides both \(a\) and \(b\) if, and only if, \(d\) also divides both \(a\) and \(b - a\). This implies that
and hence \(\boxed{\gcd(a, b) = \gcd(a, b-a)}\).
Here is a quick idea of why that is true: say \(d \mid a, b\). Then \(a/d\) and \(b/d\) are integers. So,
a difference of integers, and so an integers as well, and hence \(d \mid (b -a)\).
And, if \(d \mid a, (b-a)\), then \(a/d\) and \((b-a)/d\) are both integers, and thus
a sum of integers, and so an integer too, and therefore \(d \mid b\).
If you are not convinced by the argument above, let’s do some tests using random numbers and Sage’s gcd:
number_of_tests = 100
maxn = 1_000_000
for i in range(number_of_tests):
# a, b as two random integers from 1 to maxn
a = randint(1, maxn)
b = randint(1, maxn)
if gcd(a, b) != gcd(a, b - a):
print(f"Failed for {a = } and {b = }!")
break # get out of the loop!
else: # no break
print("It worked for all tests!")
It worked for all tests!
But how does that help? Say that we need to compute the \(\gcd(12, 15)\) again. Well, we have
Here, since we know \(3 \mid 12\), we could already see that \(\gcd(12, 15) = \gcd(12, 3) = 3\)! But, let’s pretend we did not notice that. We can repeat the idea! Replace the larger number by its difference with the smaller.
So, we can do:
This last line follows from \(\gcd(a, 0) = a\) when \(a>0\). We did not perform a single division, only subtractions.
Let’s do a larger example, say \(\gcd(159, 48)\):
This already works well, but note that by repeating the idea:
Hence, when we had \(\gcd(12, 3)\) we could subract \(4 \cdot 3 = 12\) form \(12\) to get \(\gcd(0, 3) = 3\) above.
More generally, if we do the long division \(a = bq + r\) (with \(0 \leq r \leq b\)), then
Let’s apply this same idea to the steps above:
This is more efficient! This process above it called the Euclidean Algorithm. Here are the steps.
Algorithm 3.1 (Euclidean Algorithm)
To compute \(\gcd(a, b)\):
Perform the long division of the larger by the smaller and replace the larger by the remainder. (This remainder now is smaller than the previous smallest element!)
Repeat until you get a remainder of \(0\).
The GCD is the last divisor (meaning, the last non-zero remainder/smallest element).
So, we have a series of long divisions, until we get a remainder of \(0\):
So, let’s do another example: \(\gcd(1{,}027, 349)\):
1027 % 349 # remember that '%' gives the remainder
329
So, the remainder is \(329\). The next divide the smaller number (\(349\)) by this remainder:
349 % 329
20
Now, repeat: get the remainder of the smaller of the two numbers by the new remainder:
329 % 20
9
And keep going until we get a remainder equal to \(0\):
20 % 9
2
9 % 2
1
2 % 1
0
The last divisor (or the last non-zero remainder) is the GCD, so \(\gcd(1{,}027, 349) = 1\). We can check with Sage’s own gcd function:
gcd(1027, 349)
1
Homework
In your homework you will automate this process, i.e., implement the Euclidean Algorithm as a function.
3.8. How Fast Is the Euclidean Algorithm#
With some (relatively simple) math, one can prove that if the sequence of remainder is \(r_1\), \(r_2\), \(r_3\), …, then
This means that after two long divisions, the remainder is at most half the size.
Hence, we get:
(Here we are using \(b\) as the “zero-th remainder”, since we use it as the divisor for the first long division, so our bound above still works.)
This means that \(r_{2k} < \dfrac{b}{2^k}\). Since \(2^k\) grows very fast as \(k\) increases, \(r_{2k}\) gets small very fast. With a little more math, one can show that
where, again, \(b\) is the smallest between \(a\) and \(b\).
Note
Remember that \(\log_2\) denotes the logarithm (or simply log) base \(2\).
For instance, \(\log_3(9) = 2\), since \(2\) is the power of the base, i.e., \(3\), that gives the argument, i.e., \(9\): \(3^2 = 9\).
Similarly:
\(\log_2(16) = 4\), since \(2^4 = 16\) and
\(\log_5(1/25) = -2\), since \(5^{-2} = 1/25\).
To illustrate that this is not much at all, if \(b = 1{,}000{,}000\), we need to perform at most \(2 \log_2(1{,}000{,}000) + 2\) long divisions:
2 * floor(log(10^6, 2)) + 2
40
So, at most \(40\)! And, most of the time, it will be considerably less! To illustrate, using Sage’s own function, which does use the Euclidean Algorithm:
%%time
gcd(76497326597236475295264957264957246971641974697326592743619743, 89435098282058743250982430584059850982437095214350820)
CPU times: user 10 µs, sys: 0 ns, total: 10 µs
Wall time: 10.5 µs
3
3.9. The Extended Euclidean Algorithm#
Here is a result that will be quite important to us:
Lemma 3.1 (Bezout’s Lemma)
Let \(a\) and \(b\) be positive integers. Then, there are integers \(u\) and \(v\) (maybe negative), such that:
Let’s illustrate it with same examples.
We had the \(\gcd(12, 15) = 3\), and
We had the \(\gcd(159, 48) = 3\), and
Let’s check this second one:
159 * (-3) + 48 * 10
3
We have that \(\gcd(1{,}027, 349) = 1\), and
Again, let’s check:
1027 * 157 + 349 * (-462)
1
Note
Note that these \(u\) and \(v\) are not unique!
For example:
In fact, there are infinitely many possibilities: let \(k\) be any integer. Then
So,
work (for \(a=15\) and \(b=12\)) for any integer \(k\)!
But let’s now focus on how to find one possible pair.
Suppose we want to find (some pair of) integers \(u\) and \(v\) such that
We will not even assume that we know \(\gcd(159, 48)\): we will compute it at the same time we find \(u\) and \(v\).
Here is how we can do it: start with the trivial equalities:
Now, we will perform the Euclidean Algorithm described above: we do the long division \(149 = 48 \cdot 3 + 15\). Then we subtract \(149 - 3 \cdot 48\), but we will do it with the whole equations above:
So, now we have the two equations:
Repeat, as with the Euclidean Algorithm: we divide \(48 = 15 \cdot 3 + 3\):
Note
We are using the distributive property here. More precisely, we are factoring out the common terms \({\color{green} 159}\) and \({\color{green} 48}\):
and
So, now our two equations are:
Next \(15 = 3 \cdot 5 + 0\), so we know that \(\gcd(149, 49) = 3\) and our last equation gives us \(u\) and \(v\):
These are the same \(u\) and \(v\) we’ve checked above!
Let’s do this now with the help of Sage for our computations. Let’s try to find \(u\) and \(v\) so that \(\gcd(1{,}027, 349) = 1{,}027 \cdot u + 349 \cdot v\).
Now, we always deal with two equations, so we will use x and y for the numbers on the left side, and u1 and v1 for the first equation and u2 and v2 for the second. When we start, we always have u1 = 1 and v1 = 0 and u2 = 0 and v2 = 1:
a, b = 1037, 349
# initial values
x, y = a, b
u1, v1 = 1, 0
u2, v2 = 0, 1
Now, long division:
q, r = divmod(x, y) # quotient and remainder
Now, the “new” x is the old y and the new y is the remainder r:
x, y = y, r
x, y
(349, 339)
Important
Here is very important that we use the “double assignment”, or we would need an intermediate variable!
We also need the new u’s’ and v’s. But the new u1 and v1 are just the old u2 and v2 and the new u2 and v2 are simply the u1 - q * u_2 and v1 - q * v_2:
u1, v1, u2, v2 = u2, v2, (u1 - q * u2), (v1 - q * v2)
u1, v1, u2, v2
(0, 1, 1, -2)
(Again, it is very important to keep this as a multiple assignment!)
So, we have:
Now, we repeat until we get a remainder of zero:
q, r = divmod(x, y)
x, y = y, r
u1, v1, u2, v2 = u2, v2, (u1 - q * u2), (v1 - q * v2)
print(f"{r = }, {x = }, {u1 = }, {v1 = }, {y = }, {u2 = }, {v2 = }")
r = 10, x = 339, u1 = 1, v1 = -2, y = 10, u2 = -1, v2 = 3
So:
Repeat, since the remainder was not zero:
q, r = divmod(x, y)
x, y = y, r
u1, v1, u2, v2 = u2, v2, (u1 - q * u2), (v1 - q * v2)
print(f"{r = }, {x = }, {u1 = }, {v1 = }, {y = }, {u2 = }, {v2 = }")
r = 9, x = 10, u1 = -1, v1 = 3, y = 9, u2 = 34, v2 = -101
Now:
Repeat, since the remainder was not zero:
q, r = divmod(x, y)
x, y = y, r
u1, v1, u2, v2 = u2, v2, (u1 - q * u2), (v1 - q * v2)
print(f"{r = }, {x = }, {u1 = }, {v1 = }, {y = }, {u2 = }, {v2 = }")
r = 1, x = 9, u1 = 34, v1 = -101, y = 1, u2 = -35, v2 = 104
We have:
Repeat:
q, r = divmod(x, y)
x, y = y, r
u1, v1, u2, v2 = u2, v2, (u1 - q * u2), (v1 - q * v2)
print(f"{r = }, {x = }, {u1 = }, {v1 = }, {y = }, {u2 = }, {v2 = }")
r = 0, x = 1, u1 = -35, v1 = 104, y = 0, u2 = 349, v2 = -1037
Ah-ha! Now the remainder is zero and the \(\gcd(1{,}037, 349)\) is x and \(u\) and \(v\) are u1 and v1 respectively:
x, u1, v1
(1, -35, 104)
Let’s check:
1 == 1037 * (-35) + 349 * 104
True
Success!
There is one small improvement we can make, to avoid some unnecessary computations. (This will only make any noticeable difference for very large numbers, but we will deal with very large numbers soon, so we might as well think about it.)
We can skip computing the v’s and only focus on the u’s. At the end, we can compute v from what we have.
Here is how it would go:
a, b = 1037, 349
# initial values
x, y = a, b
u1 = 1
u2 = 0
q, r = divmod(x, y)
x, y = y, r
u1, u2 = u2, (u1 - q * u2)
print(f"{r = }, {x = }, {u1 = }, {y = }, {u2 = }")
r = 339, x = 349, u1 = 0, y = 339, u2 = 1
Remainder is not zero, so continue:
q, r = divmod(x, y)
x, y = y, r
u1, u2 = u2, (u1 - q * u2)
print(f"{r = }, {x = }, {u1 = }, {y = }, {u2 = }")
r = 10, x = 339, u1 = 1, y = 10, u2 = -1
Repeat:
q, r = divmod(x, y)
x, y = y, r
u1, u2 = u2, (u1 - q * u2)
print(f"{r = }, {x = }, {u1 = }, {y = }, {u2 = }")
r = 9, x = 10, u1 = -1, y = 9, u2 = 34
Repeat:
q, r = divmod(x, y)
x, y = y, r
u1, u2 = u2, (u1 - q * u2)
print(f"{r = }, {x = }, {u1 = }, {y = }, {u2 = }")
r = 1, x = 9, u1 = 34, y = 1, u2 = -35
Repeat:
q, r = divmod(x, y)
x, y = y, r
u1, u2 = u2, (u1 - q * u2)
print(f"{r = }, {x = }, {u1 = }, {y = }, {u2 = }")
r = 0, x = 1, u1 = -35, y = 0, u2 = 349
Remainder is \(0\), so \(\gcd(1{,}037, 349) = 1\) and \(u = -35\). But what is \(v\)?
Well
so, solving for \(v\):
So, we get:
v = (x - a * u1) // b
v
104
This indeed gives us the same \(u\) and \(v\) as above!
3.9.1. Summary (Algorithm)#
Here is the actual extended Euclidean Algorithm:
Algorithm 3.2 (Extended Euclidean Algorithm)
To compute the GCD of a and b and u and v such that gcd(a, b) = a * u + b * v:
Initialize
x, y, u1, u2 = a, b, 1, 0 r = x % y
While
ris not zero, repeat the stepsq, r = divmod(x, y) x, y, u1, u2 = y, r, u2, (u1 - q * u2)
When
ris zero, we have that:xis the GCD ofaandb;u1is the requiredu;(x - a * u1) // bis the requiredv.
Homework
You will implement this algorithm in your homework.
Of course, Sage has its own function for this, xgcd:
xgcd(1037, 349)
(1, -35, 104)
This tells us that the GCD is \(1\), \(u=-35\), and \(v=104\), the same we’ve found above.
3.10. Other Solutions#
As we’ve seen above, the \(u\) and \(v\) from Bezout’s Lemma are not unique. But, we actually know how to find all solutions!
Theorem 3.3
Let \(a\) and \(b\) be positive integers, \(d\) be their GCD, and \(u_0\) and \(v_0\) be some pair of integers such that
Then, \(u\) and \(v\) are integers such that
if, and only if,
for some integer \(k\).
Note
The same integer \(k\) must be use for both \(u\) and \(v\).
Note that \(b/d\) and \(a/d\) are both integers since \(d\) is a common divisor.
Note that the signs are different for \(u\) and \(v\).
Let’s understand what this is saying. Remember that we’ve found:
According to the theorem we can produce infinitely many other solutions. For instance, using \(k=1\), we get
Let’s check that they also work:
159 * 13 + 48 * (-43)
3
Let’s make Sage check for us:
bound = 100
a, b, u0, v0 = 159, 48, -3, 10
d = gcd(a, b)
for k in range(-bound, bound + 1):
u = u0 + k * (b // d)
v = v0 - k * (a // d)
if a * u + b * v != d: # fail!
print(f"Formula fails for {k = }.")
break
else: # no break
print("It seems to work!")
It seems to work!
It is actually easy to see why. If
and
then
On the other hand, the theorems tells us more! We’ve just checked that if \(u_0\) and \(v_0\) give us the GCD, then so do
But the theorem tells us that those are all the solutions, i.e., every single solution can be obtained by choosing some integer \(k\) for the formula.
It is not too hard the check that this is true with some algebra, but let’s just do a computer check instead. Note that if \(au + bv = d\) with \(u\) and \(v\) integers, then
so \(b \mid (d - au)\). So, to search for solutions, we try integers \(u\) in some range, and if \(b \mid (d-au)\), then the formula above would gives us \(v\).
max_number = 10^4 # max size of a and b
repetitions = 100 # number of tries
max_k = 30 # maximum k from the formula; go from -max_k to max_k
for _ in range(repetitions):
# random a and b in the range
a = randint(2, max_number)
b = randint(2, max_number)
# get GCD and initial u and v
d, u0, v0 = xgcd(a, b)
# solutions from the formula
formula_solutions = {(u0 + k * (b // d), v0 - k * (a // d)) for k in range(-max_k, max_k + 1)}
# initialize found solutions as empty set
found_solutions = set()
# find solutions with u in the same range used by the formula
for u in range(u0 + (-max_k) * (b // d), u0 + max_k * (b // d) + 1):
# if the division is exact, then the quotient is v
v, r = divmod(d - a * u, b)
if r == 0:
# add solution found ones!
found_solutions.add((u, v))
# check if the sets match!
if found_solutions != formula_solutions:
print("Failed!")
# give information of example for which it failed
print(f"{a = }, {b = }, {d = }, {u0 = }, {v0 = }.")
# does the formula miss solutions?
print(f"Solutions found, but not with the in formula:\n {found_solutions - formula_solutions}\n\n")
# does the formula gives solutions we did not found?
print(f"Solutions from the formula, but not found:\n {formula_solutions - found_solutions}")
break
else: # no break
print("It seems to work!")
It seems to work!
3.10.1. Finding a Specific Solution#
Often we want a specific solution \(u\) and \(v\). Say that want the smallest positive \(u\).
For instance, \(\gcd(293, 58) = 1\) and
a, b = 293, 58
u0, v0 = -135, 682
d = gcd(a, b)
d
1
d == a * u0 + b * v0
True
The previous theorem tells us that we if add (since \(u_0\) is negative) \(b/d = 58\) (an integer!) to \(u_0\) and subtract \(a/d = 293\) to \(v_0\), we will have another solution! (Here we are assuming that \(a\) and \(b\) are positive!)
# initialize
u, v = u0, v0
# replace
u = u + b // d
v = v - a // d
print(u)
-77
But it is still negative. So, we repeat until u becomes positive. At that point it will be the smallest one!
u = u + b // d
v = v - a // d
print(u)
-19
Still negative…
u = u + b // d
v = v - a // d
print(u)
39
Ah-ha! It is positive! The corresponding v is:
v
-197
So, \(\boxed{u=39, \; v = -197}\) is the solution (of \(\gcd(293, 58) = 293 u + 58 v\)) with the least positive \(u\).
What if we start with an already positive \(u_0\), say \(u_0 = 155\) and \(v_0 = -783\):
u0, v0 = 155, -783
# check:
d == a * u0 + b * v0
True
In this case we subtract \(b / d\) from \(u_0\) and add \(a / d\) to \(v_0\) until we get a negative \(u\). At that point we went one step to far. So, we add \(b / d\) from \(u_0\) and subtract \(a / d\) to \(v_0\) to get our solution. Let’s see it in practice:
# initialize
u, v = u0, v0
# replace: note the signs!
u = u - b // d
v = v + a // d
print(u)
97
Still positive, so repeat:
u = u - b // d
v = v + a // d
print(u)
39
Repeat:
u = u - b // d
v = v + a // d
print(u)
-19
Negative! So, our previous solution was the one we needed. We can get it back:
# switch addition/subtraction
u = u + b // d
v = v - a // d
u, v
(39, -197)
Of course, it is the same solution we found before. (We just started at a different initial solution.)
Homework
In your homework, you will automate this process: you will write a function that given a particular solution \((u_0, v_0)\), it finds the one with the least positive \(u\).
You can do it as above: have two cases, depending on whether \(u_0\) is positive or negative. But there is a clever way to find this smallest positive \(u\) directly using long division! Can you find it?
Hint
What happens if you divide \(u_0\) by \(b/d\)?
3.11. The Fundamental Theorem of Arithmetic#
The following result is an important property of primes:
Lemma 3.2 (Euclid’s Lemma)
If \(p\) is a prime and \(p \mid a \cdot b\), then either \(p \mid a\) or \(p \mid b\) (or both, as in mathematics/logic, “or” is not exclusive).
More generally, if \(p \mid a_1 \cdot a_2 \cdots a_k\), then \(p\) divides (at least) one of the \(a_i\)’s in the product.
This a defining property of primes and shows how they are indivisible parts of integers: if \(p\) is a part of \(a \cdot b\), then it is a part of either \(a\) or \(b\). We cannot “create” \(p\) from other parts from \(a\) and \(b\). For composite numbers, this can certainly happen. For example \(6 \mid 9 \cdot 10 = 90\) (as \(90 = 6 \cdot 15\)), but \(6 \mid 9\) and \(6 \nmid 15\). In this case we are “making” \(6\) from a factor of \(2\) from \(10\) and a factor of \(3\) from \(9\).
We can prove it using Bezout's Lemma:
Proof. Suppose that \(p \mid a \cdot b\), but \(p \nmid a\). We then need to show that \(p \mid b\).
Since \(p \nmid a\), we have that \(\gcd(a, p) = 1\), and by Bezout’s Lemma, there are integers \(u\) and \(v\) such that
Multiplying this by \(b\), we get
Now \(p\) divides \(ab\), so it divides \(uab\). Also, clearly it divides \(vpb\). Then it divides the right-side of the equation above, so it divides the left-side as well, i.e., \(p \mid b\).
The case for more than two factors can be done by breaking down the product in twos. For example, if \(p \mid abc = (ab) \cdot c\), then, by the case above with two factors, we have that either \(p \mid ab\) or \(p \mid c\). The former gives us that \(p\) divides either \(a\) or \(b\), so we have that \(p\) divides either \(a\), \(b\), or \(c\).
(The proper way to prove this more general case is to use Mathematical Induction, but the idea is hopefully clear.)
With Euclid's Lemma one can prove an important fact about primes you might remember from school:
Theorem 3.4 (Fundamental Theorem of Arithmetic (Prime Factorization))
Let \(n\) be an integers greater than or equal to \(2\). Then, there are primes \(p_1 < p_2 < \cdots < p_k\), for some \(k \geq 1\) and positive integers \(e_1, e_2, \ldots, e_k\) such that
Moreover, this representation is unique!
For example,
Hence, prime numbers are the (indivisible) building blocks for integers.
It is not too hard to prove this theorem, but we will skip for the sake of brevity. Instead, we will illustrate it with direct computations, e.g., writing a function to compute the factorization of a given integer.
Homework
You will write a function to compute the prime decomposition of an integer in your homework.
Although you will write your own function for prime factorization in your homework, of course, Sage has its own:
for year in range(2024, 2027):
print(f"{year} = {factor(year)}")
2024 = 2^3 * 11 * 23
2025 = 3^4 * 5^2
2026 = 2 * 1013
Note
You might remember from school how to compute GCDs using prime factorization. Although it works well for small numbers, the Euclidean Algorithm is a lot more efficient in practice!
3.12. Integers vs Real Numbers#
At first glance, it seems that since the integers are simpler than real numbers, problems involving integers should be easier than problems involving real numbers, but in fact, the opposite is true. In fact, the real numbers were created to “make things easier”, e.g., to add solutions to equations that do not have solutions with integers (or rational numbers), like \(x^2 = 2\).
As another example, consider finding all real solutions of
That is relatively easy: if \(z=0\), we only have the solution \((0,0,0)\). If not, can divide by \(z^2\) and get
a point in the unit circle. So, there is some angle \(\theta \in [0, 2 \pi)\) such that
Therefore, all solutions are then \((z \cos(\theta), z \sin(\theta), z)\) for all pairs of \(\theta \in [0, 2\pi)\) and real number \(z\).
But how about if we want \(x\), \(y\), and \(z\) be integers? We can still solve this problem, but it is a lot harder! We can spot some easy ones, like \((3k, 4k, 5k)\) for all integers \(k\), but it requires some clear ideas to find them all, like we easily did for real numbers.
Below, we illustrate other problems involving integers that are easy to state and understand but that are still quite difficult.
3.13. Some Open Problems#
3.13.1. The Goldbach Conjecture#
A conjecture is a statement believed to be true, but not proven to be true. In general, this means we have a lot of numerical evidence, but no mathematical proof. Here is a famous conjecture in Number Theory:
Conjecture (Goldbach Conjecture)
Every even integer greater than or equal to four is a sum of two prime numbers.
For instance, \(18 = 3 + 5\), or \(60 = 7 + 53\).
Let’s write a function that given an integer \(n \geq 4\), outputs two primes that add to it, or if there is no such pair, prints a message saying that such pair does not exist.
def goldbach(n):
"""
Given an even integer n >= 4, finds a pair of primes that add to n.
INPUT:
n: an even integer greater than or equal to 4.
OUTPUT:
Two primes that add to n.
"""
p = 2
while not is_prime(n - p) and p <= n / 2: # NOTE THE n / 2!
p = next_prime(p)
if p > n / 2:
print(f"{n} is a counter example!") # not going to happen!!!
return None
# print(f"{n} = {p} + {n - p}")
return p, n - p
Let’s test it for integers between \(4\) and \(100\):
for n in range(4, 101, 2):
p, q = goldbach(n)
print(f"{n:>3} = {p:>3} + {q:>3}")
4 = 2 + 2
6 = 3 + 3
8 = 3 + 5
10 = 3 + 7
12 = 5 + 7
14 = 3 + 11
16 = 3 + 13
18 = 5 + 13
20 = 3 + 17
22 = 3 + 19
24 = 5 + 19
26 = 3 + 23
28 = 5 + 23
30 = 7 + 23
32 = 3 + 29
34 = 3 + 31
36 = 5 + 31
38 = 7 + 31
40 = 3 + 37
42 = 5 + 37
44 = 3 + 41
46 = 3 + 43
48 = 5 + 43
50 = 3 + 47
52 = 5 + 47
54 = 7 + 47
56 = 3 + 53
58 = 5 + 53
60 = 7 + 53
62 = 3 + 59
64 = 3 + 61
66 = 5 + 61
68 = 7 + 61
70 = 3 + 67
72 = 5 + 67
74 = 3 + 71
76 = 3 + 73
78 = 5 + 73
80 = 7 + 73
82 = 3 + 79
84 = 5 + 79
86 = 3 + 83
88 = 5 + 83
90 = 7 + 83
92 = 3 + 89
94 = 5 + 89
96 = 7 + 89
98 = 19 + 79
100 = 3 + 97
This has been tested up to really huge numbers, so it is believed to be true, but no one has proved it yet.
3.13.2. The Twin Primes Conjecture#
Definition (Twin Primes)
Two primes are called twin primes if they differ by \(2\). This means that, apart from \(2\) and \(3\), they are as close as possible, since they will both be odd.
For example, \(11\) and \(13\) are twin primes, and \(101\) and \(103\) are as well.
Conjecture 3.1 (Twin Primes Conjecture)
There are infinitely many pairs of twin primes. In other words, given any integer \(N\) (no matter how large), there is a prime \(p \geq N\) such that \(p + 2\) is also prime.
So, let’s write a function that, given \(N\), finds a pair of twin primes with the smaller at least \(N\). Here we use the Sage’s handy function next_prime.
def twin_test(N):
"""
Given N, finds a pair of twin primes larger than n.
"""
p = next_prime(N)
while not is_prime(p + 2):
p = next_prime(p)
return p, p + 2
Let’s test it:
N = 10^12 # very large
twin_test(N)
(1000000000061, 1000000000063)