iii. Primes and Sieves

Sieving flips “for each $n$, look at its divisors” into “for each $d$, touch its multiples”: ${\sum }_{d\le n}n/d=O(n\log n)$ over all $d$, ${\sum }_{p\le n}n/p=O(n\log \log n)$ over primes, $O(n)$ for the linear sieve (each composite touched once, by its smallest prime factor). Testing one big number: Miller–Rabin; counting primes up to $10^{11}$: Sublinear Sieves.

$n$	$10^3$	$10^5$	$10^6$	$10^7$	$10^8$	$10^9$
$\pi (n)$	168	9592	78498	664579	5761455	50847534

$\pi (n)\sim n/\ln n$; the largest prime gap below $10^9$ is 282 and below $10^{18}$ is under 1500, so “scan upward until a prime” is cheap.

Eratosthenes, plain and fast

vector<int> eratosthenes(int n) {
    vector<char> comp(n + 1); vector<int> pr;
    for (int i = 2; i <= n; i++) if (!comp[i]) { pr.push_back(i); for (ll j = (ll)i * i; j <= n; j += i) comp[j] = 1; }
    return pr;
}
vector<int> fastSieve(int lim) {                           // primes <= lim, odd-only, cache-sized blocks (1e9 in ~1.5 s)
    const int S = (int)round(sqrt(lim)), R = lim / 2;
    vector<int> pr{2}; if (lim < 2) return {}; pr.reserve(int(lim / log(lim) * 1.1));
    vector<char> sv(S + 1); vector<pair<int, int>> cp;     // (prime, next odd index to cross out)
    for (int i = 3; i <= S; i += 2) if (!sv[i]) { cp.push_back({i, i * i / 2}); for (int j = i * i; j <= S; j += 2 * i) sv[j] = 1; }
    vector<char> block(S);
    for (int L = 1; L <= R; L += S) {
        fill(all(block), 0);
        for (auto& [p, idx] : cp) for (int i = idx; i < S + L; idx = (i += p)) block[i - L] = 1;
        for (int i = 0; i < min(S, R - L + 1); i++) if (!block[i]) pr.push_back((L + i) * 2 + 1);
    }
    while (sz(pr) && pr.back() > lim) pr.pop_back();
    return pr;
}
vector<char> segmentedSieve(ll L, ll R) {                 // isPrime for [L, R], R up to ~1e12, R - L up to ~1e7
    vector<char> ok(R - L + 1, 1);
    for (ll p : eratosthenes((int)sqrtl(R) + 1)) for (ll j = max(p * p, (L + p - 1) / p * p); j <= R; j += p) ok[j - L] = 0;
    for (ll x = L; x <= min(R, 1LL); x++) ok[x - L] = 0;
    return ok;
}

Proof: Eratosthenes and its cost

A composite $n$ has a prime factor $p\le \sqrt{n}$, and $n\ge p^2$ is a multiple of $p$ reached by $p$‘s loop, so it is crossed out; primes are never crossed out because every crossed number is a proper multiple. Starting at $p^2$ is safe because smaller multiples $kp$ with $k<p$ have a prime factor $<p$ and were crossed earlier. Cost: ${\sum }_{p\le n}n/p=O(n\log \log n)$ by Mertens’ theorem (${\sum }_{p\le n}1/p=\ln \ln n+O(1)$). The fast version stores only odd numbers (index $i$ represents $2i+1$), keeps for each base prime the next index to cross (stepping by $p$ in index space is stepping by $2p$ in value space), and processes blocks of $\sqrt{n}$ indices that fit in cache; correctness is the same argument restricted to odd values. The segmented sieve crosses, for every prime $p\le \sqrt{R}$, its multiples $\ge \max (p^2,L)$ in the window.

Linear sieve and multiplicative functions

struct LinearSieve {                                      // spf, primes, and the smallest-prime-power part of each n
    vector<int> spf, primes, pw, ex;                      // n = pw[n] * rest with pw[n] = spf(n)^ex[n], gcd(pw, rest) = 1
    LinearSieve(int n) : spf(n + 1), pw(n + 1, 1), ex(n + 1) {
        for (int i = 2; i <= n; i++) {
            if (!spf[i]) spf[i] = i, primes.push_back(i), pw[i] = i, ex[i] = 1;
            for (int p : primes) {
                if (p > spf[i] || (ll)i * p > n) break;
                int x = i * p; spf[x] = p;
                if (p == spf[i]) pw[x] = pw[i] * p, ex[x] = ex[i] + 1; else pw[x] = p, ex[x] = 1;
            }
        }
    }
    template <class T, class F> vector<T> multiplicative(F fpe) const {   // fpe(p, e) = f(p^e); f(1) = 1
        int n = sz(spf) - 1; vector<T> f(n + 1); if (n >= 1) f[1] = 1;
        for (int i = 2; i <= n; i++) f[i] = pw[i] == i ? fpe(spf[i], ex[i]) : f[i / pw[i]] * f[pw[i]];
        return f;
    }
    vector<pair<int, int>> factor(int x) const { vector<pair<int, int>> r; while (x > 1) { int p = spf[x], e = 0; while (x % p == 0) x /= p, e++; r.push_back({p, e}); } return r; }
};
// e.g. phi: S.multiplicative<ll>([](ll p, int e) { ll q = 1; for (int i = 1; i < e; i++) q *= p; return q * (p - 1); });
//      mu:  S.multiplicative<int>([](int, int e) { return e == 1 ? -1 : 0; });
//      d:   S.multiplicative<int>([](int, int e) { return e + 1; });

Proof: the linear sieve marks every composite exactly once

Every composite $x$ factors uniquely as $x=i\cdot p$ with $p=\text{spf}(x)$ and $i=x/p$, where $\text{spf}(i)\ge p$. The loop for this $i$ reaches $p$ (it runs over primes up to $\text{spf}(i)$) and sets $\text{spf}(x)=p$. Any other pair $(i^{\prime },p^{\prime })$ with $i^{\prime }{p}^{\prime }=x$ and $p^{\prime }>\text{spf}(x)$ has $\text{spf}(i^{\prime })=\text{spf}(x)<p^{\prime }$, so the loop for $i^{\prime }$ breaks before $p^{\prime }$. Hence each composite is written once and the total work is $O(n)$. Prime-power part: if $p=\text{spf}(i)$ then $x=ip$ has the same smallest prime with one more factor, so pw and ex extend; otherwise $p$ appears exactly once in $x$. Multiplicative functions: $f(n)=f(p^e)f(n/p^e)$ for $p^e\Vert n$ by multiplicativity, and $n/p^e<n$ was computed earlier.

	Eratosthenes	Linear sieve
Time	$O(n\log \log n)$	$O(n)$
Smallest prime factor	an extra pass	built in
Multiplicative functions ($\phi ,\mu ,d,\sigma ,\dots$)	a separate pass per function	computed during the sieve, for any $f$ evaluable on prime powers

Sieving over multiples

vector<ll> exactGcdPairs(const vector<int>& a, int V) {    // ans[d] = #pairs i < j with gcd(a_i, a_j) = d, values <= V
    vector<ll> freq(V + 1), cnt(V + 1), ans(V + 1);
    for (int x : a) freq[x]++;
    for (int d = 1; d <= V; d++) for (int m = d; m <= V; m += d) cnt[d] += freq[m];   // elements divisible by d
    for (int d = V; d >= 1; d--) { ans[d] = cnt[d] * (cnt[d] - 1) / 2; for (int m = 2 * d; m <= V; m += d) ans[d] -= ans[m]; }
    return ans;
}
pair<vector<int>, vector<ll>> divisorCountAndSum(int n) {  // d(k), sigma(k) for all k <= n, O(n log n)
    vector<int> d(n + 1); vector<ll> s(n + 1);
    for (int i = 1; i <= n; i++) for (int j = i; j <= n; j += i) d[j]++, s[j] += i;
    return {d, s};
}

Idiom	Cost
for each $d$, loop over multiples (divisor lists, ${\sum }_{d\mid n}$ transforms)	$O(n\log n)$
“exactly $d$” from “divisible by $d$” by subtracting multiples (descending)	$O(n\log n)$
store divisor lists for all $n\le N$	$O(N\log N)$ memory ($\approx 14N$ for $N=10^6$)
factor every $n\le N$ by spf	$O(\log n)$ each
Möbius transform (subset-sum over divisors), body f[j] += g[i]	$O(N\log N)$, a Dirichlet-style transform
count multiples of each $i$ present in an array, body cnt[i] += present[j]	$O(N\log N)$

Proof: the exact-gcd recursion

cnt[d] counts elements divisible by $d$, so $\left(\genfrac{}{}{0ex}{}{\text{cnt}[d]}{2}\right)$ counts pairs whose gcd is a multiple of $d$. These split by the exact gcd $g\in \{d,2d,\dots \}$, so the count for exactly $d$ is that binomial minus the exact counts for proper multiples, which are final when processed in decreasing order. Both loops cost ${\sum }_{d}V/d=O(V\log V)$. This is Möbius inversion over the divisibility poset (Arithmetic Functions).

Cache blocking: for $N\ge 10^8$, process the sieve array in blocks of about 32 KB so each block fits in L1 cache (fastSieve above already does this) — typically a 2-3× speedup over a flat sieve, which thrashes memory at this size.

The idea generalises well beyond primes: sieving over subsets is exactly SOS DP ($O(2^{n}n)$ over the subset lattice instead of the divisibility lattice); sieving over a range of exponents handles problems framed around $p^k$; “for each state, push to all reachable states” in a DP is the same harmonic-sum pattern whenever transitions are multiples.