Number theory: gcd, divisors, sieve, and modular powers¶

~25 min · prerequisite: Binary search and the bisect module

What you will be able to do¶

Compute the greatest common divisor with Euclid's method, and the least common multiple from it.
List the divisors of a number in about sqrt(n) steps, and find all primes up to n with a sieve.
Compute a ** e % m without building the huge number a ** e, and use the modular inverse for a prime modulus.

The idea¶

Many Codeforces problems ask for an answer modulo a number, such as 10**9 + 7. The remainder of a divided by m is written a % m. Python's % always returns a value from 0 up to m - 1 when m is positive, even when a is negative. This lesson covers four tools that appear often.

The greatest common divisor, written gcd(a, b), is the largest number that divides both a and b. Euclid's method uses the rule gcd(a, b) = gcd(b, a % b), which holds because any common divisor of a and b also divides a % b. When b becomes 0, the answer is the other number. The least common multiple follows from it: lcm(a, b) = a // gcd(a, b) * b.

The number of Euclid steps is small. Let r be the remainder at one step. If the divisor is at most half of the number, the next remainder is below the divisor, which is at most half. If the divisor is more than half, the next remainder is a - b, which is less than half of a. Either way, every two steps at least halve the remainder, so there are about 2 * log2(min(a, b)) steps.

Divisors. If d divides n, then n // d also divides n, and one of the two is at most sqrt(n). So it is enough to test d from 1 to isqrt(n), and each hit gives two divisors. That is about sqrt(n) steps. For n = 10**9 this is about 31 623 steps.

The sieve of Eratosthenes finds all primes up to n. Start with every number marked as prime, then for each prime p up to sqrt(n), mark p * p, p * p + p, and so on as not prime. The marking for p touches about n / p numbers. Adding over the primes gives n times the sum of 1 / p. A known result, Mertens' second theorem, says that this sum grows like ln(ln(n)). So the sieve costs about O(n log log n) steps.

Modular powers. The number a ** e can be enormous. Instead, use fast exponentiation: square the base at each step and multiply it into the result whenever the current bit of e is 1. A number e has about log2(e) bits, so the loop runs that many times, with at most two multiplications each. Every intermediate value is reduced modulo m, so the numbers stay below m * m.

Worked example¶

Here is Euclid's method with a trace for gcd(48, 18).

def gcd(a, b):
    while b != 0:
        a, b = b, a % b
    return a


assert gcd(48, 18) == 6
assert gcd(7, 0) == 7

Line by line, this is what Python does with a = 48 and b = 18:

The loop runs while b is not zero.
The tuple assignment a, b = b, a % b computes the right side first. So the new a is 18, and the new b is 48 % 18 = 12.
The next pass gives a = 12 and b = 18 % 12 = 6.
The next pass gives a = 6 and b = 12 % 6 = 0.
The loop stops because b is zero. The function returns a, which is 6.

The tuple assignment is what makes this short. Two separate assignments would overwrite a before the new b was computed.

Try it¶

The first block lists divisors and primes. The second checks fast exponentiation against Python's built-in three-argument pow.

import math


def divisors(n):
    found = set()
    for d in range(1, math.isqrt(n) + 1):
        if n % d == 0:
            found.add(d)
            found.add(n // d)
    return sorted(found)


def primes_up_to(n):
    if n < 2:
        return []
    flags = bytearray([1]) * (n + 1)
    flags[0] = 0
    flags[1] = 0
    for p in range(2, math.isqrt(n) + 1):
        if flags[p]:
            flags[p * p::p] = bytes(len(range(p * p, n + 1, p)))
    return [i for i in range(n + 1) if flags[i]]


assert divisors(36) == [1, 2, 3, 4, 6, 9, 12, 18, 36]
assert divisors(1) == [1]
assert len(primes_up_to(100)) == 25
assert primes_up_to(20) == [2, 3, 5, 7, 11, 13, 17, 19]
assert primes_up_to(1) == []

import random


def fast_pow(base, exp, mod):
    result = 1 % mod
    base %= mod
    while exp > 0:
        if exp & 1:
            result = result * base % mod
        base = base * base % mod
        exp >>= 1
    return result


for _ in range(2000):
    b = random.randint(0, 10**6)
    e = random.randint(0, 500)
    m = random.randint(1, 10**9)
    assert fast_pow(b, e, m) == pow(b, e, m)

p = 10**9 + 7
inv3 = pow(3, p - 2, p)
assert 3 * inv3 % p == 1
assert pow(3, -1, 7) == 5

Common mistakes¶

Building the full power before taking the remainder. The expression 2 ** (10**9) % m first creates 2 ** (10**9), a number with about a billion bits, which is roughly 125 megabytes. Use pow(2, 10**9, m) instead, which reduces at every step.
Expecting a negative remainder. Python returns a non-negative remainder when the divisor is positive. This check fails:
```
assert -7 % 3 == -1
```
```
AssertionError
```
The value is 2. If you need the C-style result, compute it yourself.
Dividing with / where an integer is needed. The operator / always returns a float, so n / d can not be used as an index or a count:
```
values = [10, 20, 30]
print(values[4 / 2])
```
```
TypeError: list indices must be integers or slices, not float
```
Use // for integer division.
Taking the remainder with a zero divisor. Euclid's method is safe because the loop stops at zero, but a direct a % 0 fails:
```
print(7 % 0)
```
```
ZeroDivisionError: division by zero
```
Check that the divisor is not zero before you use it.

Practice¶

Easy: compute gcd(84, 36) with math.gcd and check the result is 12.
Easy: count the primes up to 100 by trial division, using all(p % d for d in range(2, p)). The answer is 25.
Medium: write lcm_list(values) that returns the least common multiple of a list, using functools.reduce with math.gcd as in the lcm formula. Check lcm_list([4, 6, 10]) is 60.
Medium: for a prime p, write mod_inverse(a, p) that returns pow(a, p - 2, p). Check that mod_inverse(2, 7) is 4, because 2 * 4 = 8 and 8 % 7 = 1.
Stretch: write count_divisors(n) using the square-root loop from divisors, and check it for every n from 1 to 200 against a brute-force count.

Solution

```python import math

assert math.gcd(84, 36) == 12 ```

python primes = [p for p in range(2, 101) if all(p % d for d in range(2, p))] assert len(primes) == 25
```python from functools import reduce import math

def lcm_list(values): return reduce(lambda a, b: a // math.gcd(a, b) * b, values)

assert lcm_list([4, 6, 10]) == 60 ```

```python def mod_inverse(a, p): return pow(a, p - 2, p)

assert mod_inverse(2, 7) == 4 ```

Fermat's little theorem says that for a prime p and a not divisible by p, a ** (p - 2) is an inverse of a modulo p.

```python import math

def count_divisors(n): count = 0 for d in range(1, math.isqrt(n) + 1): if n % d == 0: count += 2 if d * d != n else 1 return count

for n in range(1, 201): brute = sum(1 for d in range(1, n + 1) if n % d == 0) assert count_divisors(n) == brute ```

A divisor below sqrt(n) pairs with one above it. The exact square root is counted once.

Recap¶

Euclid's method finds gcd in about 2 * log2(min(a, b)) steps, and lcm follows from it.
Divisors and primes need only loops up to sqrt(n), and a sieve marks composites in about n log log n steps.
Use pow(base, exp, mod), which keeps every step small, and use // rather than / for integers.