Unit 14: Numbers and Primes

Testing one number: stop at the square root

To check whether n is prime, try dividing by 2, 3, 4… but only up to the square root. If n = a * b then one of a and b is at most the square root, so a factor above it would already have been found below.

So one number costs about sqrt(n), which for a million is a thousand.

Python

def is_prime(n):
    if n < 2:
        return False
    i = 2
    while i * i <= n:
        if n % i == 0:
            return False
        i += 1
    return True

print([n for n in range(2, 30) if is_prime(n)])

It prints

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29]

Testing every number: cross out instead

Checking a million numbers one at a time costs a thousand million. The sieve of Eratosthenes does the lot in about three million: take each prime and cross out its multiples.

Start the crossing at p * p — anything smaller has a smaller factor and is already gone.

Python

limit = 30
prime = [True] * (limit + 1)
prime[0] = False
prime[1] = False

for p in range(2, limit + 1):
    if prime[p]:
        for multiple in range(p * p, limit + 1, p):
            prime[multiple] = False

print([n for n in range(limit + 1) if prime[n]])
print(sum(prime))

It prints

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29]
10

How many primes are there?

Roughly n / log n of the numbers below n are prime — 25 below 100, 168 below 1000, about 78 thousand below a million.

That matters when a problem says "for each prime": it is far fewer than n, and often the difference between fitting and not.

Every number is a product of primes, one way only

Divide out each factor as often as it goes. Whatever survives past the square root is itself a prime and goes on the end.

360 is 2 times 2 times 2 times 3 times 3 times 5. Nothing else multiplies to 360 in primes.

Python

def factorise(n):
    factors = []
    d = 2
    while d * d <= n:
        while n % d == 0:
            factors.append(d)
            n //= d
        d += 1
    if n > 1:
        factors.append(n)
    return factors

print(factorise(360))
print(factorise(97))
print(factorise(1000000007))

It prints

[2, 2, 2, 3, 3, 5]
[97]
[1000000007]

The factors tell you how many divisors there are

360 is 2 to the 3, times 3 squared, times 5. A divisor picks 0 to 3 twos, 0 to 2 threes, and 0 or 1 five — so 4 * 3 * 2 = 24 of them.

Counting them by trying every number up to 360 gives the same answer and takes 360 times longer.

Python

print(4 * 3 * 2)
print(len([d for d in range(1, 361) if 360 % d == 0]))

It prints

24
24

Euclid, still the fastest thing in the unit

The greatest common divisor of a and b is the gcd of b and a % b, until the second one is 0. It takes about log steps and it is two lines.

The lowest common multiple comes free: divide the product by the gcd — and divide first if the numbers are big, or the product may be enormous.

Python

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

print(gcd(360, 84))
print(360 // gcd(360, 84) * 84)

It prints

12
2520

Why every counting problem says "modulo 1000000007"

The answers are astronomically large. Rather than making you print a thousand digits, the setter asks for the remainder after dividing by a big prime — usually 1000000007.

Adding and multiplying survive that perfectly: you may take the remainder at every step and the answer is the same.

Python

MOD = 1000
a = 123456
b = 789

print((a * b) % MOD)
print(((a % MOD) * (b % MOD)) % MOD)

It prints

784
784

Powers by doubling

3 ** 100 under a modulus, one multiply at a time, is a hundred steps. Squaring instead does it in seven: read the exponent in binary — unit 3 again — and square as you go.

This is called fast exponentiation and it is everywhere.

Python

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

MOD = 1000000007
print(power(2, 10, MOD))
print(power(3, 100, MOD))
print(pow(3, 100, MOD))

It prints

1024
886041711
886041711

Python has it built in, and you should still know the loop

pow(base, exponent, mod) is the same thing and is what you would type in a contest. The loop is here because it is the shape of unit 10's jumping and unit 16's matrix powers — the same doubling, three times over.

And in C++ there is no built-in, so people write this loop from memory.

Subtracting can go negative

In Python -3 % 7 is 4, which is what you want. In most other languages it is -3, which is not — so contest code from elsewhere is full of (a - b + mod) % mod.

Worth knowing when you read someone else's solution and wonder what the extra + mod is for.

Python

print(-3 % 7)
print((4 - 6) % 7)

It prints

4
5

The one operation that does not survive

Adding, subtracting and multiplying all work fine with remainders. Division does not: 10 divided by 5 is 2, but the remainders of 10 and 5 tell you nothing useful.

And counting problems are full of division — every "choose k from n" is a fraction.

So multiply by the opposite instead

The inverse of a is the number that multiplies with it to give 1. Under 11, the inverse of 7 is 8, because 56 leaves 1.

Dividing by 7 is then multiplying by 8, and everything stays a whole number.

Python

MOD = 11
for a in range(1, 11):
    for b in range(1, 11):
        if a * b % MOD == 1:
            print(a, b)

It prints

1 1
2 6
3 4
4 3
5 9
6 2
7 8
8 7
9 5
10 10

Fermat's little theorem finds it for you

When the modulus is prime, a to the power mod - 1 leaves 1. So a times a to the power mod - 2 leaves 1 — which makes a ** (mod - 2) the inverse.

That is one call to fast exponentiation, and it is why contests always pick a prime modulus.

Python

MOD = 1000000007
inverse = pow(3, MOD - 2, MOD)
print(inverse)
print(3 * inverse % MOD)
print(pow(7, 9, 11), 7 * pow(7, 9, 11) % 11)

It prints

333333336
1
8 1

Now binomial coefficients work

"How many ways to choose 50 from 100" is a fraction of factorials with about 160 digits in it. Under the modulus it is three numbers multiplied together.

Work out all the factorials once, then every question is two inverses and a multiply.

Python

MOD = 1000000007
limit = 200

fact = [1] * (limit + 1)
for i in range(1, limit + 1):
    fact[i] = fact[i - 1] * i % MOD

def choose(n, k):
    return fact[n] * pow(fact[k], MOD - 2, MOD) % MOD * pow(fact[n - k], MOD - 2, MOD) % MOD

print(choose(10, 3))
print(choose(20, 10))
print(choose(100, 50))

It prints

120
184756
538992043

Why the answer looks like nonsense

choose(100, 50) really is a 30-digit number; 538992043 is its remainder. That is the answer the judge wants, and there is no way to look at it and tell it is right.

So check your code on choose(10, 3), where you can count 120 by hand. Always test a modular program on a case small enough to verify.