AIOJ

For the definitions of primes and composite numbers, see Elementary number theory.

The prime-counting function counts the primes less than or equal to $x$ and is denoted by $\pi(x)$ . As $x$ grows, it satisfies the approximation $\pi(x) \sim \dfrac{x}{\ln(x)}$ .

Primality tests

A primality test can be used to determine whether a given natural number is prime.

There are two types of primality tests:

Deterministic tests: determine with certainty whether a number is prime. Common examples include trial division, the Lucas–Lehmer test, and elliptic-curve primality proving.
Probabilistic tests: usually much faster than deterministic tests, but they may (although with very small probability) incorrectly identify a composite number as prime (the converse does not happen). A number that passes a probabilistic primality test is therefore called a probable prime until its primality is proved deterministically. A number that passes the test but is actually composite is called a pseudoprime. There are many specific types of pseudoprimes; the most common are Fermat pseudoprimes, which are composite numbers satisfying Fermat's little theorem. A common example of a probabilistic test is the Miller–Rabin test.

Trial division

The brute-force approach can simply enumerate numbers in increasing order and check whether each divides the given number.

Reference implementation

C++

bool isPrime(int a) {
  if (a < 2) return false;
  for (int i = 2; i < a; ++i)
    if (a % i == 0) return false;
  return true;
}

Python

def isPrime(a):
    if a < 2:
        return False
    for i in range(2, a):
        if a % i == 0:
            return False
    return True

This approach is quite reliable, but is it really necessary to test every number?

It is easy to observe that if $x$ is a divisor of $a$ , then $\frac{a}{x}$ is also a divisor of $a$ .

This tells us that for each pair $(x, \frac{a}{x} )$ , it is enough to test only one member. For convenience, consider the smaller member of each pair. It is easy to see that all these smaller members lie in the interval $[1, \sqrt{a}]$ .

Since $1$ is certainly a divisor, we do not test it.

Reference implementation

C++

bool isPrime(int a) {
  if (a < 2) return 0;
  for (int i = 2; (long long)i * i <= a; ++i)  // prevent overflow
    if (a % i == 0) return 0;
  return 1;
}

Python

def isPrime(a):
    if a < 2:
        return False
    for i in range(2, int(sqrt(a)) + 1):
        if a % i == 0:
            return False
    return True

Fermat primality test

The Fermat primality test is the simplest probabilistic primality test.

From Fermat's little theorem, we obtain the following idea for testing primality:

The basic idea is to repeatedly choose a base in $[2, n-1]$ , denoted $a$ , and check whether $a^{n-1} \equiv 1 \pmod n$ holds every time.

Reference implementation

C++

bool fermat(int n) {
  if (n < 3) return n == 2;
  // test_time is the number of tests; it is recommended to set it to at least 8
  // to ensure accuracy, but not too large, otherwise efficiency will suffer
  for (int i = 1; i <= test_time; ++i) {
    int a = rand() % (n - 2) + 2;
    if (quickPow(a, n - 1, n) != 1) return false;
  }
  return true;
}

Python

def fermat(n):
    if n < 3:
        return n == 2
    # test_time is the number of tests; it is recommended to set it to at least 8
    # to ensure accuracy, but not too large, otherwise efficiency will suffer
    for i in range(1, test_time + 1):
        a = random.randint(0, 32767) % (n - 2) + 2
        if quickPow(a, n - 1, n) != 1:
            return False
    return True

If $a^{n−1} \equiv 1 \pmod n$ but $n$ is not prime, then $n$ is called a Fermat pseudoprime to the base $a$ . In practice, when $a^{n−1} \equiv 1 \pmod n$ , $n$ is usually prime. However, there are counterexamples: for $n = 341$ and $a = 2$ , we have $2^{340}\equiv 1 {\pmod {341}}$ , but $341 = 11 \cdot 31$ is composite. In fact, for every fixed base $a$ , there are infinitely many such counterexamples¹.

Since the Fermat test cannot guarantee correctness for a single base, a natural idea is to check several bases. However, even checking every possible $n$ -coprime base $a$ still cannot guarantee that $n$ is prime. In other words, the converse of Fermat's little theorem is false: even if for every $a\perp n$ we have $a^{n-1}\equiv 1\pmod n$ , $n$ need not be prime. Such numbers are called Carmichael numbers. There are infinitely many of them as well. This forces us to seek a stricter primality test.

Miller–Rabin primality test

The Miller–Rabin primality test is a better method for testing primality. It was developed by Miller and Rabin as an improvement of the Fermat test. Like other probabilistic primality tests, it can only detect pseudoprimes. A much slower deterministic algorithm is needed to guarantee primality. However, no known number has passed the Miller–Rabin test or another advanced probabilistic test while actually being composite, so it is safe to use in practice.

Ignoring the cost of multiplication, testing a number $n$ for $k$ rounds takes $O(k \log n)$ time. Miller–Rabin is often used for high-precision integers, where the time complexity is $O(k \log^3n)$ ; techniques such as FFT can optimize it to $O(k \log^2n \log \log n \log \log \log n)$ .

To address the challenge posed by Carmichael numbers, the Miller–Rabin test further uses the following property of primes:

Theorem on square roots of 1

If $p$ is an odd prime, then the solutions to $x^2 \equiv 1 \pmod p$ are $x \equiv 1 \pmod p$ and $x \equiv p - 1 \pmod p$ .

Proof

It is easy to verify that when $p$ is an odd prime, both $x\equiv 1\pmod p$ and $x\equiv p-1\pmod p$ satisfy the equation above. By Lagrange's theorem, these are all its solutions.

Combining Fermat's little theorem with the theorem on square roots of one gives the Miller–Rabin primality test:

Factor the exponent of the equation $a^{n-1} \equiv 1 \pmod n$ , namely $n−1$ , as $n−1=u \times 2^t$ ;
In each round, first choose a random $a$ , compute $v = a^{u} \bmod n$ , and then square this value at most $t$ times;
During the process, if a nontrivial square root of $1$ is found (that is, a root other than $\pm 1$ ), the number is not prime;
Otherwise, use the Fermat primality test.

There are also a few implementation details:

In one round, if $a^{u \times 2^s} \equiv n-1 \pmod n$ at some point, all subsequent squarings will yield $1$ , so the round can be accepted immediately.
If a nontrivial square root $a^{u \times 2^s} \not\equiv n-1 \pmod n$ is found, all subsequent squarings will yield $1$ . One may return false immediately or wait until after $t$ squarings to return false.

This gives a more correct Miller–Rabin implementation (from fjzzq2002):

Reference implementation

C++

bool millerRabin(int n) {
  if (n < 3 || n % 2 == 0) return n == 2;
  if (n % 3 == 0) return n == 3;
  int u = n - 1, t = 0;
  while (u % 2 == 0) u /= 2, ++t;
  // test_time is the number of tests; it is recommended to set it to at least 8
  // to ensure accuracy, but not too large, otherwise efficiency will suffer
  for (int i = 0; i < test_time; ++i) {
    // 0, 1, and n-1 pass directly; a is chosen from [2, n-2]
    int a = rand() % (n - 3) + 2, v = quickPow(a, u, n);
    if (v == 1) continue;
    int s;
    for (s = 0; s < t; ++s) {
      if (v == n - 1) break;  // found the trivial square root n-1; pass this round
      v = (long long)v * v % n;
    }
    // If a nontrivial square root is found, the loop reaches s == t because it cannot break early
    // If the Fermat test fails, v never equals -1 before the loop reaches s == t
    if (s == t) return 0;
  }
  return 1;
}

Python

def millerRabin(n):
    if n < 3 or n % 2 == 0:
        return n == 2
    if n % 3 == 0:
        return n == 3
    u, t = n - 1, 0
    while u % 2 == 0:
        u = u // 2
        t = t + 1
    # test_time is the number of tests; it is recommended to set it to at least 8
    # to ensure accuracy, but not too large, otherwise efficiency will suffer
    for i in range(test_time):
        # 0, 1, and n-1 pass directly; a is chosen from [2, n-2]
        a = random.randint(2, n - 2)
        v = pow(a, u, n)
        if v == 1:
            continue
        s = 0
        while s < t:
            if v == n - 1:
                break
            v = v * v % n
            s = s + 1
        # If a nontrivial square root is found, the loop reaches s == t because it cannot break early
        # If the Fermat test fails, v never equals -1 before the loop reaches s == t
        if s == t:
            return False
    return True

It can be proved² that for an odd composite number $n > 9$ , the probability of passing the Miller–Rabin test for one randomly chosen base $a$ is at most one fourth. Therefore, after choosing $k$ bases at random, the probability of misclassifying a composite number as prime is at most $1/4^k$ .

Proof

Let $n-1=u2^t$ , where $u$ is odd and $t$ is a positive integer. Then the integer $n$ passes the Miller–Rabin test with base $a$ exactly when

a^u\equiv 1{\textstyle\pmod n},\text{ or }a^{u2^i}\equiv -1{\textstyle\pmod n}\text{ for some }0\le i < t.

Let $a$ range over the set $S$ (more precisely, its residue classes). We will show that

|S| \le \dfrac14\varphi(n).

Here, $\varphi(n)$ is the Euler totient function. The proof has three steps.

Step 1: Let $\ell$ be the largest positive integer such that $2^\ell \mid p-1$ for every prime factor of $n$ , denoted $p$ . We can prove that

S\subseteq S' = \{a\bmod n:a^{u2^{\ell-1}}\equiv\pm 1{\textstyle\pmod n}\}.

For $S$ , an element $a$ has only two possibilities. If $a^u\equiv 1\pmod n$ , then clearly $a^{u2^{\ell-1}}\equiv 1\pmod n$ , so $a\in S'$ . If, for some $0\le i < t$ , $a^{u2^i}\equiv -1\pmod n$ , then for every prime factor $p\mid n$ , $a^{u2^i}\equiv-1\pmod p$ . Let $\delta_p(a)$ denote the order of $a$ modulo $p$ ; see order. Then $\delta_p(a)\mid u2^{i+1}$ but $\delta_p(a)\nmid u2^{i}$ , so in the prime factorization of $\delta_p(a)$ , the exponent of $2$ is exactly $i+1$ , and hence $2^{i+1}\mid\delta_p(a)$ . Fermat's little theorem gives $\delta_p(a)\mid p-1$ , so $2^{i+1}\mid p-1$ . This holds for every prime factor of $n$ , denoted $p$ , so $i+1\le\ell$ . Consequently, $a^{u2^{\ell-1}} = (a^{u2^i})^{2^{\ell-1-i}} \equiv \pm 1 \pmod n$ , and again $a\in S'$ . Combining the two cases gives $S\subseteq S'$ .

Step 2: Compute $|S'|$ .

Suppose that the prime factorization of $n$ is $n = p_1^{e_1}p_2^{e_2}\cdots p_k^{e_k}$ . By the Chinese remainder theorem, $a^{u2^{\ell - 1}}\equiv 1\pmod n$ is equivalent to $a^{u2^{\ell - 1}}\equiv 1\pmod{p_i^{e_i}}$ for every $p_i^{e_i}$ . Since every odd prime power $p_i^{e_i}$ has a primitive root, the number of solutions to $a^{u2^{\ell - 1}}\equiv 1\pmod{p_i^{e_i}}$ , given by the number of solutions, is

\gcd(u2^{\ell-1},p_i^{e_i-1}(p_i-1)) = \gcd(u2^{\ell-1},p_i-1) = 2^{\ell-1}\gcd(u,p_i-1).

The first equality holds because $u$ is a divisor of $n-1$ and therefore cannot be divisible by $p_i$ ; the second holds by the choice of $\ell$ . Thus, by the Chinese remainder theorem, the number of solutions to $a^{u2^{\ell-1}}\equiv 1\pmod n$ is

\prod_{p\mid n}2^{\ell-1}\gcd(u,p-1).

Similarly, $a^{u2^{\ell - 1}}\equiv -1\pmod n$ is equivalent to $a^{u2^{\ell - 1}}\equiv -1\pmod{p_i^{e_i}}$ for every $p_i^{e_i}$ . For any factor $p_i^{e_i}$ , the condition $a^{u2^{\ell - 1}}\equiv -1\pmod{p_i^{e_i}}$ is equivalent to $a^{u2^{\ell - 1}}\not\equiv 1\pmod{p_i^{e_i}}$ and $a^{u2^{\ell}}\equiv 1\pmod{p_i^{e_i}}$ . As above, the number of solutions to $a^{u2^{\ell}}\equiv 1\pmod{p_i^{e_i}}$ is $2^{\ell}\gcd(u,p_i-1)$ , so the number of solutions to $a^{u2^{\ell - 1}}\equiv -1\pmod{p_i^{e_i}}$ is also

2^{\ell}\gcd(u,p_i-1) - 2^{\ell-1}\gcd(u,p_i-1) = 2^{\ell-1}\gcd(u,p_i-1).

Applying the Chinese remainder theorem again, the number of solutions to the congruence $a^{u2^{\ell - 1}}\equiv -1\pmod n$ is

\prod_{p\mid n}2^{\ell-1}\gcd(u,p-1).

Combining the two cases gives

|S'| = 2\prod_{p\mid n}2^{\ell-1}\gcd(u,p-1).

Step 3: Prove that $|S'|\le\varphi(n)/4$ .

Using the expression $\varphi(n)=\prod_ip_i^{e_i-1}(p_i-1)$ for the Euler totient function, we have

\dfrac{\varphi(n)}{|S'|} = \dfrac{1}{2}\prod_ip_i^{e_i-1}\dfrac{p_i-1}{2^{\ell-1}\gcd(u,p_i-1)}.

For every $i$ , the corresponding factor $p_i^{e_i-1}\dfrac{p_i-1}{2^{\ell-1}\gcd(u,p_i-1)}$ is even, so $\varphi(n)/|S'|$ is an integer. Suppose that $|S'|\le\varphi(n)/4$ is false. Then necessarily $\varphi(n)/|S'|=1,2,3$ , that is,

\prod_ip_i^{e_i-1}\dfrac{p_i-1}{2^{\ell-1}\gcd(u,p_i-1)} = 2,4,6.

Since every factor in the product is even, the product either has one factor equal to $2,4,6$ , or has two factors both equal to $2$ .

First consider the case of two factors. Neither factor has an odd prime factor, so $p_i^{e_i-1}=1$ ; that is, $n$ is square-free. Assume $n=p_1p_2$ with primes $p_1<p_2$ . Both factors equal $2$ , so $p_i-1=2^{\ell}\gcd(u,p_i-1)$ always holds. Hence $p_i=1+2^\ell m_i$ , where $m_i$ is odd and $m_i\mid u$ . Reducing $p_1p_2=n=1+u2^t$ modulo $m_1$ gives $p_1p_2\equiv 1\pmod{m_1}$ , and hence $p_2\equiv 1\pmod{m_1}$ , so $m_1\mid m_2$ . The reverse divisibility also holds. Thus $m_1=m_2$ , which means $p_1=p_2$ , contradicting $p_1<p_2$ . This case is impossible.

Finally, consider the case of only one factor, namely the composite number $n=p^e$ with $e>1$ . Then $p^{e-1}\mid 2,4,6$ must hold. The only possibility is $p=3,e=2$ , namely $n=9$ , contradicting the assumption of the proposition. This case is impossible as well.

All cases together show that $|S'|\le\varphi(n)/4$ .

Combining the three steps gives $|S|\le |S'|\le \varphi(n)/4$ for every odd composite $n>9$ .

In addition, assuming the Generalized Riemann hypothesis (GRH), for a number $n$ , it suffices to test every integer in $[2, \min\{n-2, \lfloor 2\ln^2 n \rfloor\}]$ to deterministically establish whether $n$ is prime.³

In the OI setting, primality is usually tested for numbers in $[1, 2^{64})$ . For numbers in $[1, 2^{32})$ , the three bases $\{2, 7, 61\}$ make Miller–Rabin deterministic; for numbers in $[1, 2^{64})$ , the seven bases $\{2, 325, 9375, 28178, 450775, 9780504, 1795265022\}$ suffice.⁴

One may also use $\{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37\}$ , the first $12$ primes, to test primes in $[1, 2^{64})$ .

When using a number $a$ from the sequences above as a base to test whether $n$ is prime, note that:

Every number must be used; do not choose only those less than $n$ ;
Replace $a$ by $a \bmod n$ ;
If $a \equiv 0 \pmod n$ or $a \equiv \pm 1 \pmod n$ , accept the round immediately.

Highly composite numbers

As the name suggests, a prime has exactly two divisors. A highly composite number is a number with the greatest number of divisors (and the smallest value among numbers with the same number of divisors), so the notion is relative to a set.

An intuitive definition is: in a set of positive integers, the number with the most divisors and the smallest value is a highly composite number.

Highly composite number

For a positive integer $n$ , if every positive integer smaller than $n$ has fewer divisors than $n$ , then it is called a highly composite number (anti-prime, also known as a highly composite number).

Note

Do not confuse this with an emirp, a prime whose digit reversal is a different prime (for example, both 149 and 941 are emirps, while 101 is not).

Procedure

How can we find highly composite numbers?

Since we need the number of divisors, the first step is prime factorization. Write $n$ as $n=p_{1}^{k_{1}}p_{2}^{k_{2}} \cdots p_{n}^{k_{n}}$ , where $p$ is prime and $k$ is its exponent. The total number of divisors is then $(k_1+1) \times (k_2+1) \times (k_3+1) \cdots \times (k_n+1)$ .

However, prime factorization is clearly expensive, and the result for one number cannot be reused for the next. We therefore need another method.

Let us examine the properties of highly composite numbers.

A highly composite number is necessarily a product of powers of consecutive primes starting from $2$ .
Smaller primes have exponents at least as large as larger primes: in $n=p_{1}^{k_{1}}p_{2}^{k_{2}} \cdots p_{n}^{k_{n}}$ , we have $k_1 \geq k_2 \geq k_3 \geq \cdots \geq k_n$ .

Explanation:

If the primes are not consecutive starting from $2$ , replacing a prime by a smaller one while keeping its exponent leaves the divisor count unchanged but decreases $n$ . When the primes are exchanged to become consecutive starting from $2$ , the value of $n$ is smallest.
If a smaller prime has a smaller exponent than a larger prime, swapping the two primes while keeping their exponents leaves the number of divisors of $n$ unchanged but decreases $n$ .

There are two further questions:

For a given $n$ , how far into the primes should we enumerate?

The most extreme case is $n=p_{1}p_{2} \cdots p_{n}$ , so it is enough to multiply consecutive primes until the product is just at most $n$ . If we enumerate a larger prime, some earlier prime must have exponent $0$ , so the result cannot be highly composite.
How many powers should we enumerate?

Consider an extreme case. If a power of the smallest prime already exceeds the given maximum $n$ , then in every other form the largest exponent is smaller than that power. In the extreme case where $n$ is a power of $2$ , it is enough to enumerate up to $\lfloor\log_2 n\rfloor$ .

With these details in place, how do we implement the search?

Each stage before considering a prime can be viewed as a node in a tree, and we search level by level. When should the search stop?

The current number is already greater than the target number;
The factor currently being enumerated can no longer be used;
The current divisor count is greater than the target divisor count;
The current divisor count is exactly the target count (at this point, determine whether the minimum $\mathit{ans}$ should be updated).

Then the DFS can continue iterating downward, enumerating the exponent one level at a time.

Examples

Codeforces 27E. A number with a given number of divisors

Find the smallest natural number with a given number of divisors. The answer is guaranteed not to exceed $10^{18}$ .

Solution idea

For this type of problem, use the divisor count as the DFS stopping condition and continually update the smallest value found.

Reference code

#include <iostream>
constexpr int p[15] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47};
// The product of the primes in p exceeds 1e18

int n;
long long ans = 2e18;

// u: index in p of the prime currently under consideration
// num: value constructed so far
// cnt: number of divisors of the current value
// pre: exponent of the previous factor, limiting the exponent chosen now
void dfs(int u, long long num, long long cnt, int pre) {
  if (cnt > n || u >= 15) return;
  if (cnt == n) return ans = std::min(ans, num), void();
  for (int i = 1; i <= pre; ++i) {
    if (num * p[u] > ans) break;  // pruning
    dfs(u + 1, num *= p[u], cnt * (i + 1), i);
  }
}

int main() {
  std::cin >> n;
  dfs(0, 1, 1, 59);  // floor(log2(1e18))=19
  std::cout << ans << std::endl;
  return 0;
}

ZOJ 2562 More Divisors

Among the numbers not exceeding $n$ , find the one with the most divisors.

Solution idea

The idea is the same as above; only the DFS stopping condition needs to be changed. Pay attention to the data range, since a 32-bit integer may overflow.

Reference code

#include <iostream>

int p[16] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53};
unsigned long long n;
unsigned long long ans,
    ans_num;  // ans is the largest highly composite number at most n (updated continuously); ans_sum is
              // its number of divisors.

// depth: which prime is currently being enumerated
// temp: value when the current divisor count is num
// num: current divisor count
// up: exponent of the previous prime, limiting the current exponent
void dfs(int depth, unsigned long long temp, unsigned long long num, int up) {
  if (depth >= 16 || temp > n) return;
  if (num > ans_num) {  // update the answer
    ans = temp;
    ans_num = num;
  }
  if (num == ans_num && ans > temp) ans = temp;  // update the answer
  for (int i = 1; i <= up; i++) {
    if (temp * p[depth] > n)
      break;  // pruning: if multiplying by this factor exceeds ans, it cannot be optimal
    dfs(depth + 1, temp *= p[depth], num * (i + 1),
        i);  // take one copy of this factor and search for the next factor
  }
  return;
}

using std::cin;
using std::cout;

int main() {
  cin.tie(nullptr)->sync_with_stdio(false);
  while (cin >> n) {
    ans_num = 0;
    dfs(0, 1, 1, 60);
    cout << ans << '\n';
  }
  return 0;
}

References and notes

Rui-Juan Jing, Marc Moreno-Maza, Delaram Talaashrafi, "Complexity Estimates for Fourier-Motzkin Elimination", Journal of Functional Programming 16:2 (2006), pp. 197–217.
Number theory, Part I: Primes and primality testing
Notes on Miller–Rabin and Pollard–Rho - Bill Yang's Blog
Primality test - Wikipedia
Fermat pseudoprime - Wikipedia
Taozi's algorithm notes — an introduction to highly composite numbers (ACM/OI)
The Rabin-Miller Primality Test
Highly composite number - Wikipedia

Footnotes

Theorem 1 of Pomerance, Carl, John L. Selfridge, and Samuel S. Wagstaff, "The pseudoprimes to 25⋅ 10⁹," Mathematics of Computation 35, no. 151 (1980): 1003–1026, shows that for any fixed base $a$ , infinitely many composite numbers can also pass the stronger Miller–Rabin primality test. ↩
This result and its proof are based on Section 3.5 of Crandall, Richard, and Carl Pomerance. Prime numbers: a computational perspective. New York, NY: Springer New York, 2005. ↩
Bach, Eric, "Explicit bounds for primality testing and related problems", Mathematics of Computation, 55:191 (1990), pp. 355–380. ↩
For more results of this kind, see Deterministic variant of the Miller–Rabin primality test. ↩