If you've never encountered any number theory before, AoPS is a good place to start.

Resources

Prime Factorization

A positive integer $a$ is called a divisor or a factor of a non-negative integer $b$ if $b$ is divisible by $a$, which means that there exists some integer $k$ such that $b = ka$. An integer $n > 1$ is prime if its only divisors are $1$ and $n$. Integers greater than $1$ that are not prime are composite.

Every positive integer has a unique prime factorization: a way of decomposing it into a product of primes, as follows:

where the $p_i$ are distinct primes and the $a_i$ are positive integers.

Now, we will discuss how to find the prime factorization of any positive integer.

Copy
vector<int> factor(int n) {
	vector<int> ret;
	for (int i = 2; i * i <= n; i++) {
		while (n % i == 0) {
			ret.push_back(i);
			n /= i;
		}
	}
	if (n > 1) { ret.push_back(n); }
	return ret;
}

Copy
List<Integer> factor(int n) {
	List<Integer> factors = new ArrayList<>();
	for (int i = 2; i * i <= n; i++) {
		while (n % i == 0) {
			factors.add(i);
			n /= i;
		}
	}
	if (n > 1) { factors.add(n); }
	return factors;
}

Copy
from typing import List


def factor(n: int) -> List[int]:
	ret = []
	i = 2
	while i * i <= n:
		while n % i == 0:
			ret.append(i)
			n //= i
		i += 1
	if n > 1:
		ret.append(n)
	return ret

This algorithm runs in $\mathcal{O}(\sqrt{n})$ time, because the for loop checks divisibility for at most $\sqrt{n}$ values. Even though there is a while loop inside the for loop, dividing $n$ by $i$ quickly reduces the value of $n$, which means that the outer for loop runs less iterations, which actually speeds up the code.

Let's look at an example of how this algorithm works, for $n = 252$.

At this point, the for loop terminates, because $i$ is already 3 which is greater than $\lfloor \sqrt{7} \rfloor$. In the last step, we add $7$ to the list of factors $v$, because it otherwise won't be added, for a final prime factorization of $\{2, 2, 3, 3, 7\}$.

Solution - Counting Divisors

The most straightforward solution is just to do what the problem asks us to do - for each $x$, find the number of divisors of $x$ in $\mathcal{O}(\sqrt x)$ time.

Copy
#include <iostream>

using namespace std;

int main() {
	int n;
	cin >> n;
	for (int q = 0; q < n; q++) {
		int x;
		int div_num = 0;

Copy
import java.io.BufferedReader;
import java.io.IOException;
import java.io.InputStreamReader;

public class Divisors {
	public static void main(String[] args) throws IOException {
		BufferedReader read = new BufferedReader(new InputStreamReader(System.in));
		int queryNum = Integer.parseInt(read.readLine());
		StringBuilder ans = new StringBuilder();
		for (int q = 0; q < queryNum; q++) {

Copy
ans = []
for _ in range(int(input())):
	div_num = 0
	x = int(input())
	i = 1
	while i * i <= x:
		if x % i == 0:
			div_num += 1 if i**2 == x else 2
		i += 1
	ans.append(div_num)

print("\n".join(str(i) for i in ans))

This solution runs in $\mathcal{O}(n \sqrt x)$ time, which is just fast enough to get AC. However, we can actually speed this up to get an $\mathcal{O}((x + n) \log x)$ solution!

First, let's discuss an important property of the prime factorization. Consider:

Then the number of divisors of $x$ is simply $(a_1 + 1) \cdot (a_2 + 1) \cdots (a_k + 1)$.

Why is this true? The exponent of $p_i$ in any divisor of $x$ must be in the range $[0, a_i]$ and each different exponent results in a different set of divisors, so each $p_i$ contributes $a_i + 1$ to the product.

$x$ can have $\mathcal{O}(\log x)$ distinct prime factors, so if we can find the prime factorization of $x$ efficiently, we can use it with the above property to answer queries in $\mathcal{O}(\log x)$ time instead of the previous $\mathcal{O}(\sqrt x)$ time.

Here's how we find the prime factorization of $x$ in $\mathcal{O}(\log x)$ time with $\mathcal{O}(x \log x)$ preprocessing:

1. For each $k \leq 10^6$, find any prime number that divides $k$. To find this, we can use the Sieve of Eratosthenes which runs in $\mathcal{O}(n \log n)$, where $n$ is the larger numbers we consider. There's also a version of the sieve that runs in linear time, but we won't be needing it.
2. We can find the prime factorization of $x$ by repeatedly dividing it by the prime numbers we calculated earlier until $x = 1$.

Using this method gives us the following code:

Copy
#include <iostream>

using namespace std;

const int MAX_N = 1e6;
// max_div[i] contains the largest prime that goes into i
int max_div[MAX_N + 1];

int main() {
	for (int i = 2; i <= MAX_N; i++) {

Copy
import java.io.BufferedReader;
import java.io.IOException;
import java.io.InputStreamReader;

public class Divisors {
	private static final int MAX_N = (int)Math.pow(10, 6);
	public static void main(String[] args) throws IOException {
		// maxDiv[i] contains the largest prime that can divide i
		int[] maxDiv = new int[MAX_N + 1];
		for (int i = 2; i <= MAX_N; i++) {

Copy
MAX_N = 10**6

# max_div[i] contains the largest prime that can go into i
max_div = [0 for _ in range(MAX_N + 1)]
for i in range(2, MAX_N + 1):
	if max_div[i] == 0:
		for j in range(i, MAX_N + 1, i):
			max_div[j] = i

ans = []

GCD & LCM

GCD

The greatest common divisor (GCD) of two integers $a$ and $b$ is the largest integer that is a factor of both $a$ and $b$. In order to find the GCD of two non-negative integers, we use the Euclidean Algorithm, which is as follows:

This algorithm can be implemented with a recursive function as follows:

Copy
public int gcd(int a, int b) { return b == 0 ? a : gcd(b, a % b); }

Copy
int gcd(int a, int b) { return b == 0 ? a : gcd(b, a % b); }

Copy
def gcd(a: int, b: int) -> int:
	return a if b == 0 else gcd(b, a % b)

This function runs in $\mathcal{O}(\log ab)$ time because $a\le b \implies b\pmod a <\frac{b}{2}$.

The worst-case scenario for the Euclidean algorithm is when $a$ and $b$ are consecutive Fibonacci numbers $F_n$ and $F_{n + 1}$. In this case, the algorithm will calculate that $\gcd(F_n, F_{n + 1}) = \gcd(F_{n - 1}, F_n) = \dots = \gcd(0, F_1)$. This takes a total of $n+1$ calls, which is proportional to $\log \left(F_n F_{n+1}\right)$.

LCM

The least common multiple (LCM) of two integers $a$ and $b$ is the smallest integer divisible by both $a$ and $b$. The LCM can be calculated with the GCD using this property:

Also, these two functions are associative, meaning that if we want to take the GCD or LCM of more than two elements, we can do so two at a time, in any order. For example,

Euler's Totient Function

Properties

Euler's totient function - written using phi $\phi(n)$ - counts the number of positive integers in the interval $[1,n]$ which are coprime to $n$. Two numbers $a$ and $b$ are coprime if their greatest common divisor is equal to 1 i.e. $gcd(a,b)=1$.

Here are the values of $\phi(n)$ for the first 20 numbers:

The totient function is multiplicative meaning that $\phi(nm)=\phi(n) \cdot \phi(m)$, where $n$ and $m$ are coprime - $gcd(n, m)=1$. For example $\phi(15)=\phi(3 \cdot 5)=\phi(3) \cdot \phi(5) = 2 \cdot 4 = 8$.

Let's take a look at some edge cases for $\phi(n)$:

• If n is a prime number then $\phi(n)=n-1$ because $gcd(n, x)=1$ for all $1 \leq x < n$
• If n is a power of a prime number, $n=p^q$ where p is a prime number and $1 \leq q$ then
  there are exactly $p^{q-1}$ numbers divisible by $p$, so $\phi(p^q)=p^{q} - p^{q-1} = p^{q-1}(p - 1)$

Using the multiplicative property and the last edge case we can compute the value of $\phi(n)$ from the factorization of number $n$. Let the factorization be $n=p_1^{q_1} \cdot p_2^{q_2} \cdot \ldots \cdot p_k^{q_k}$ where $p_i$ is a prime factor of $n$, then:

Below is an implementation for factorization in $\mathcal{O}(\sqrt{n})$. It can be a little bit tricky to understand why we subtract $ans/p$ from $ans$. For example $ans=p^q \cdot x$ ,where $p$ is a prime factor and $x$ is the rest of the prime factorization. By subtracting $\frac{ans}{p}=p^{q-1} \cdot x$ we end up with: $p^q \cdot x - p^{q-1} \cdot x = p^{q-1} \cdot x \cdot (p - 1)$ which is exactly the form of $\phi(n)$ described a few lines above.

Copy
int phi(int n) {
	int ans = n;
	for (int p = 2; p * p <= n; p++) {
		if (n % p == 0) {
			while (n % p == 0) { n /= p; }
			ans -= ans / p;
		}
	}
	if (n > 1) { ans -= ans / n; }
	return ans;
}

Copy
public static int phi(int n) {
	int ans = n;
	for (int p = 2; p * p <= n; p++) {
		if (n % p == 0) {
			while (n % p == 0) { n /= p; }
			ans -= ans / p;
		}
	}

	if (n > 1) { ans -= ans / n; }
	return ans;
}

Copy
def phi(n: int) -> int:
	ans = n
	p = 2
	while p * p <= n:
		if n % p == 0:
			while n % p == 0:
				n //= p
			ans -= ans // p
		p += 1

	if n > 1:
		ans -= ans // n

	return ans

Usually in problems we need to precompute the totient of all numbers between $1$ and $n$, then factorizing is not efficient. The idea is the same as the Sieve of Eratosthenes. Since it's almost the same with the Sieve of Eratosthenes the time complexity will be: $\mathcal{O}(N\log\log{N})$.

Copy
void precompute() {
	for (int i = 1; i < MAX_N; i++) { phi[i] = i; }
	for (int i = 2; i < MAX_N; i++) {
		// If i is prime
		if (phi[i] == i) {
			for (int j = i; j < MAX_N; j += i) { phi[j] -= phi[j] / i; }
		}
	}
}

Copy
public static void precompute() {
	for (int i = 1; i < MAX_N; i++) { phi[i] = i; }

	for (int i = 2; i < MAX_N; i++) {
		// If i is prime
		if (phi[i] == i) {
			for (int j = i; j < MAX_N; j += i) { phi[j] -= phi[j] / i; }
		}
	}
}

Copy
def precompute():
	for i in range(1, MAX_N):
		phi[i] = i

	for i in range(2, MAX_N):
		# If i is prime
		if phi[i] == i:
			for j in range(i, MAX_N, i):
				phi[j] -= phi[j] // i

Solution

We are asked to compute the sum of GCDs for all pairs $(i, j)$ such that $1 \le i < j \le n$:

Let's define a helper function $f(n)$ which computes the sum of GCDs for a fixed second element $n$:

This function sums $\gcd(i, n)$ for all $i \le n$. The terms where $\gcd(i, n) = d$ are exactly those where $d$ divides $n$ and $\gcd(\frac{i}{d}, \frac{n}{d}) = 1$. The number of such integers $i$ is given by Euler's totient function $\phi(\frac{n}{d})$. Thus, we can rewrite $f(n)$ as a sum over the divisors of $n$:

We can compute $f(n)$ for all $n$ up to $10^6$ efficiently. Instead of factoring every number, we iterate through every possible divisor $i$ and update all its multiples $j$. For a fixed $i$, we add the contribution $i \cdot \phi(\frac{j}{i})$ to $f(j)$ for all $j$ that are multiples of $i$.

Finally, the problem asks for pairs where $i < j$. The value $f(j)$ includes the case $i=j$ (where $\gcd(j, j) = j$), which we must exclude. The answer for a given $n$ is the prefix sum of these adjusted values:

The total complexity of this precomputation is bounded by the harmonic series sum: $\sum_{i=1}^{N} \frac{N}{i} \approx N \ln N$

Copy
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
const int MAX_N = 1e6;

ll phi[MAX_N + 1];
ll f[MAX_N + 1];
ll sum[MAX_N + 1];

int main() {

Problems