Problem 347: Largest Integer Divisible by Two Primes

Mathematical Foundation

Theorem 1 (Structure of $p$-$q$-smooth numbers). The set of positive integers whose only prime factors are $p$ and $q$ is exactly $\{p^a q^b : a \ge 0, b \ge 0\} \setminus \{1\}$. The subset with both primes present is $\{p^a q^b : a \ge 1, b \ge 1\}$.

Proof. By the fundamental theorem of arithmetic, the prime factorization of any integer in this set involves only $p$ and $q$. Requiring both to appear forces $a \ge 1$ and $b \ge 1$. $\square$

Lemma 1 (Finite Enumeration). For fixed $p < q$ with $pq \le N$, the number of integers $p^a q^b \le N$ with $a \ge 1, b \ge 1$ is at most $\log_p(N) \cdot \log_q(N/p)$. For $N = 10^7$, this is $O(\log^2 N)$.

Proof. The exponent $a$ ranges from 1 to $\lfloor \log_p(N/q) \rfloor$ and for each $a$, $b$ ranges from 1 to $\lfloor \log_q(N/p^a) \rfloor$. Both bounds are at most $\log_2(N) \approx 23$. $\square$

Theorem 2 (Correctness of Exhaustive Search). $M(p,q,N) = \max\{p^a q^b : a \ge 1, b \ge 1, p^a q^b \le N\}$, and this maximum is attained since the set is finite and non-empty (as $pq \le N$).

Proof. The set $\{p^a q^b \le N : a,b \ge 1\}$ is non-empty because $pq \le N$ by hypothesis. It is finite because $p^a q^b \ge 2^{a+b}$ grows exponentially. A finite non-empty subset of $\mathbb{Z}$ has a maximum. $\square$

Lemma 2 (Prime Pair Bound). The number of prime pairs $(p,q)$ with $p < q$ and $pq \le N$ is $O(N / \log^2 N)$.

Proof. For each prime $p$, the number of primes $q \in (p, N/p]$ is $O(N/(p \log N))$ by the prime number theorem. Summing over primes $p \le \sqrt{N}$ gives $O(N/\log^2 N)$. $\square$

Editorial

We but p * next_prime > N means no valid q. Finally, find max p^a * q^b <= N with a >= 1, b >= 1.

Pseudocode

but p * next_prime > N means no valid q
Find max p^a * q^b <= N with a >= 1, b >= 1
For this power_p = p^a, find max b

Complexity Analysis

• Time: $O\!\left(\frac{N}{\log^2 N} \cdot \log^2 N\right) = O(N)$. There are $O(N/\log^2 N)$ prime pairs, each requiring $O(\log^2 N)$ work to enumerate powers.
• Space: $O(N)$ for the prime sieve.

Answer

#include <bits/stdc++.h>
using namespace std;

int main() {
    const long long N = 10000000;

    // Sieve of Eratosthenes
    vector<bool> is_prime(N + 1, true);
    is_prime[0] = is_prime[1] = false;
    for (long long i = 2; i * i <= N; i++) {
        if (is_prime[i]) {
            for (long long j = i * i; j <= N; j += i)
                is_prime[j] = false;
        }
    }

    vector<long long> primes;
    for (long long i = 2; i <= N; i++) {
        if (is_prime[i]) primes.push_back(i);
    }

    long long total = 0;

    for (int i = 0; i < (int)primes.size(); i++) {
        long long p = primes[i];
        if (p * p > N) break; // p < q, so p*q > p*p

        for (int j = i + 1; j < (int)primes.size(); j++) {
            long long q = primes[j];
            if (p * q > N) break;

            long long best = 0;
            // Enumerate p^a * q^b with a>=1, b>=1
            long long pa = p;
            while (pa * q <= N) {
                long long val = pa * q;
                // Multiply by q as much as possible
                while (val * q <= N) val *= q;
                best = max(best, val);
                pa *= p;
            }

            total += best;
        }
    }

    cout << total << endl;
    return 0;
}

def solve():
    N = 10**7

    # Sieve of Eratosthenes
    is_prime = [False, False] + [True] * (N - 1)
    for i in range(2, int(N**0.5) + 1):
        if is_prime[i]:
            for j in range(i * i, N + 1, i):
                is_prime[j] = False

    primes = [i for i in range(2, N + 1) if is_prime[i]]

    total = 0
    for i, p in enumerate(primes):
        if p * p > N:
            break
        for j in range(i + 1, len(primes)):
            q = primes[j]
            if p * q > N:
                break

            best = 0
            pa = p
            while pa * q <= N:
                val = pa * q
                while val * q <= N:
                    val *= q
                best = max(best, val)
                pa *= p

            total += best

    print(total)

if __name__ == "__main__":
    solve()