vii. Sublinear Sieves and Prime Counting

Prefix sums ${\sum }_{i\le n}f(i)$ of multiplicative or prime-supported functions for $n$ up to $10^{11}$–$10^{13}$, all built on the floor-value set $\{\lfloor n/i\rfloor \}$ of size $O(\sqrt{n})$ (Floor Sums).

Method	Computes	Time	Space
Lucy (Legendre / Meissel style)	$\pi (n)$, ${\sum }_{p\le n}{p}^k$	$O(n^{3/4}/\log n)$	$O(\sqrt{n})$
Min_25	$\sum f(i)$, $f$ multiplicative with $f(p)$ a low-degree polynomial	$O(n^{3/4}/\log n)$ (practical to $10^{11}$)	$O(\sqrt{n})$
Du’s sieve	$\sum f$ when $f\ast g$ and $g$ have easy prefix sums	$O(n^{2/3})$	$O(n^{2/3})$
powerful-number trick	$\sum f$ when $f$ agrees with an easy $g$ on primes	$O(\sqrt{n})$ plus $G$ evaluations	$O(\sqrt{n})$
squarefree count	$\#\{i\le n:i\text{ squarefree}\}$	$O(\sqrt{n})$	$O(\sqrt{n})$
Meissel–Lehmer / LMO	$\pi (n)$ to $10^{13}$+	$O(n^{2/3})$	$O(n^{1/3})$

Lucy’s prime counting

$S(v,p)$ = numbers in $[2,v]$ that are prime or have smallest prime factor $>p$. Start with $S(v)=v-1$ and, for each prime $p\le \sqrt{n}$ in increasing order, remove the numbers whose smallest prime factor is exactly $p$: $S(v)-=S(\lfloor v/p\rfloor )-S(p-1)$ for every $v\ge p^2$. At the end $S(v)=\pi (v)$ for all floor values $v$. Weighting by $i^k$ gives prime power sums.

struct Lucy {                                              // pi(v) for every v in {floor(n / i)}
    ll n, r; vector<ll> V, S; vector<int> lo, hi;
    int id(ll v) const { return v <= r ? lo[v] : hi[n / v]; }
    Lucy(ll n) : n(n) {
        r = sqrtl(n); while (r * r > n) r--; while ((r + 1) * (r + 1) <= n) r++;
        lo.assign(r + 2, 0); hi.assign(r + 2, 0);
        for (ll l = 1; l <= n; l = n / (n / l) + 1) { ll v = n / l; (v <= r ? lo[v] : hi[n / v]) = sz(V); V.push_back(v); S.push_back(v - 1); }
        for (ll p = 2; p <= r; p++) {
            if (S[id(p)] == S[id(p - 1)]) continue;         // p is not prime
            ll sp = S[id(p - 1)], p2 = p * p;
            for (int i = 0; i < sz(V) && V[i] >= p2; i++) S[i] -= S[id(V[i] / p)] - sp;   // V is decreasing
        }
    }
    ll pi(ll v) const { return S[id(v)]; }
};

Proof: Lucy's recurrence

Let $S_j(v)$ count $x\in [2,v]$ that are prime or have all prime factors $>p_j$. Going from $p_{j-1}$ to $p=p_j$ removes exactly the composites whose smallest prime factor is $p$: $x=py$ with $y\le v/p$ and $y$ having no prime factor $<p$, $y\ge p$. Those $y$ are counted by $S_{j-1}(\lfloor v/p\rfloor )$ minus the primes below $p$, i.e. $S_{j-1}(p-1)$. For $v<p^2$ nothing changes. After all primes $\le \sqrt{n}$, every composite $\le v\le n$ has been removed (its smallest prime factor is $\le \sqrt{v}$), so $S(v)=\pi (v)$. Updating $V$ in decreasing order uses $S(\lfloor v/p\rfloor )$ before it is changed in this round (smaller values come later). Cost: ${\sum }_{p\le \sqrt{n}}\#\{v\ge p^2\}=O(n^{3/4}/\log n)$.

Min_25 sieve

For multiplicative $f$ whose value on primes is a polynomial ${\sum }_k{c}_k{p}^k$ and with $f(p^e)$ computable: phase 1 runs Lucy for each power $p^k$ to get $G(v)={\sum }_{p\le v}f(p)$ at all floor values; phase 2 adds composites by smallest prime factor:
$S(x,j)=G(x)-G(p_{j-1})+{\sum }_{k\ge j,p_k^2\le x}{\sum }_{e\ge 1,p_k^{e+1}\le x}(f(p_k^e)S(\lfloor x/p_k^e\rfloor ,k+1)+f(p_k^{e+1}))$, and ${\sum }_{i\le n}f(i)=S(n,1)+f(1)$.

template <class FPE> mint min25(ll n, const vector<mint>& coef, FPE fpe) {   // f(p) = sum coef[k] p^k; fpe(p, e) = f(p^e)
    ll r = sqrtl(n); while (r * r > n) r--; while ((r + 1) * (r + 1) <= n) r++;
    vector<ll> V; vector<int> lo(r + 2), hi(r + 2);
    for (ll l = 1; l <= n; l = n / (n / l) + 1) { ll v = n / l; (v <= r ? lo[v] : hi[n / v]) = sz(V); V.push_back(v); }
    auto id = [&](ll v) { return v <= r ? lo[v] : hi[n / v]; };
    vector<int> primes; { vector<char> c(r + 1); for (ll i = 2; i <= r; i++) if (!c[i]) { primes.push_back(i); for (ll j = i * i; j <= r; j += i) c[j] = 1; } }
    int K = sz(coef); vector<mint> G(sz(V), 0);
    for (int k = 0; k < K; k++) {                          // Lucy on p^k: g(v) = sum_{2 <= i <= v} i^k, then sieve
        vector<mint> g(sz(V));
        for (int i = 0; i < sz(V); i++) {
            mint v = V[i];
            if (k == 0) g[i] = v - 1; else if (k == 1) g[i] = v * (v + 1) / 2 - 1; else g[i] = v * (v + 1) * (2 * v + 1) / 6 - 1;
        }
        vector<mint> pp(sz(primes) + 1, 0);
        for (int j = 0; j < sz(primes); j++) pp[j + 1] = pp[j] + mint(primes[j]).pow(k);
        for (int j = 0; j < sz(primes); j++) {
            ll p = primes[j]; mint pk = mint(p).pow(k);
            for (int i = 0; i < sz(V) && V[i] >= p * p; i++) g[i] -= pk * (g[id(V[i] / p)] - pp[j]);
        }
        for (int i = 0; i < sz(V); i++) G[i] += coef[k] * g[i];
    }
    vector<mint> Gp(sz(primes) + 1, 0);                    // sum of f(p) over the first j primes
    for (int j = 0; j < sz(primes); j++) { mint s = 0, pw = 1; for (int k = 0; k < K; k++) s += coef[k] * pw, pw *= primes[j]; Gp[j + 1] = Gp[j] + s; }
    auto S = [&](auto& self, ll x, int j) -> mint {        // sum of f(i) over 2 <= i <= x with spf(i) >= primes[j]
        if (j < sz(primes) && primes[j] > x) return 0;
        mint res = G[id(x)] - Gp[j];
        for (int k = j; k < sz(primes) && (ll)primes[k] * primes[k] <= x; k++) {
            ll p = primes[k], pe = p;
            for (int e = 1; pe * p <= x; e++, pe *= p) res += fpe(p, e) * self(self, x / pe, k + 1) + fpe(p, e + 1);
        }
        return res;
    };
    return S(S, n, 0) + 1;                                 // + f(1)
}

Supported directly: $\phi$ ($f(p)=p-1$: coef = {-1, 1}), $\mu$ ({-1}), ${\sigma }_1$ ({1, 1}), $d$ ({2}), sums of primes ({0, 1} with phase 1 only). Powers up to $k=2$ are built in; add more power-sum formulas for higher degrees.

Proof: Min_25's two phases

Phase 1 is Lucy applied to the completely multiplicative weight $i^k$: removing composites with smallest prime factor $p$ from ${\sum }_{i\le v}{i}^k$ subtracts $p^k(g(\lfloor v/p\rfloor )-{\sum }_{q<p}{q}^k)$ because $(py{)}^k=p^k{y}^k$; at the end $g(v)={\sum }_{p\le v}{p}^k$, and linearity in the coefficients gives $G(v)={\sum }_{p\le v}f(p)$. Phase 2: every $i\ge 2$ with smallest prime factor $\ge p_j$ is either a prime $\ge p_j$ (the $G(x)-G(p_{j-1})$ term) or $p_k^e\cdot y$ with $p_k\ge p_j$ the smallest prime factor, $p_k^e\Vert i$, and either $y=1$ with $e\ge 2$ (the $f(p_k^{e+1})$ terms, shifted index) or $y\ge 2$ with smallest prime factor $>p_k$ (the recursive term, multiplicative because $\gcd (p_k^e,y)=1$). Composite $i\le x$ need $p_k^2\le x$. The recursion visits only floor values, and its cost is $O(n^{3/4}/\log n)$ in the standard analysis (practically $O(n^{1-\epsilon })$ fast).

Du’s sieve (generic)

If $f\ast g=h$ with prefix sums $S_g$, $S_h$ cheap: $g(1)S_f(n)=S_h(n)-{\sum }_{d=2}^{n}g(d)S_f(\lfloor n/d\rfloor )$. Sieve $S_f$ up to $n^{2/3}$, memoise above. The concrete $\Phi$ / Mertens version is in Arithmetic Functions.

template <class T> struct DuGeneric {                      // small[x] = S_f(x) for x <= L; Sg, Sh prefix sums of g, h; g(1) = 1
    vector<T> small; function<T(ll)> Sg, Sh; unordered_map<ll, T> memo;
    T operator()(ll n) {
        if (n < sz(small)) return small[n];
        if (auto it = memo.find(n); it != memo.end()) return it->second;
        T r = Sh(n);
        for (ll l = 2, rr; l <= n; l = rr + 1) { rr = n / (n / l); r -= (Sg(rr) - Sg(l - 1)) * (*this)(n / l); }
        return memo[n] = r;
    }
};

Powerful-number trick

If $f$ is multiplicative and $g$ is multiplicative with $f(p)=g(p)$ on primes and an easy prefix sum $G$, then $h=f\ast g^{-1}$ vanishes on non-powerful numbers (every prime exponent $\ge 2$), so ${\sum }_{i\le n}f(i)={\sum }_{d\text{ powerful}}h(d)G(\lfloor n/d\rfloor )$; there are only $O(\sqrt{n})$ powerful numbers up to $n$, enumerated by DFS over primes $\le \sqrt{n}$. $h(p^e)$ comes from $f(p^e)={\sum }_{i=0}^{e}g(p^i)h(p^{e-i})$.

template <class FPE, class GPE, class GS> mint powerfulSum(ll n, FPE fpe, GPE gpe, GS G) {
    ll r = sqrtl(n); while (r * r > n) r--; while ((r + 1) * (r + 1) <= n) r++;
    vector<int> pr; { vector<char> c(r + 1); for (ll i = 2; i <= r; i++) if (!c[i]) { pr.push_back(i); for (ll j = i * i; j <= r; j += i) c[j] = 1; } }
    mint res = 0;
    auto dfs = [&](auto& self, int k, ll d, mint hd) -> void {   // d powerful, h(d) = hd; primes from index k on
        res += hd * G(n / d);
        for (int j = k; j < sz(pr) && (i128)d * pr[j] * pr[j] <= n; j++) {
            ll p = pr[j]; vector<mint> h{1};               // h(p^0..)
            ll pe = p * p;
            for (int e = 2; (i128)d * pe <= n; e++, pe *= p) {
                while (sz(h) <= e) { int t = sz(h); mint v = fpe(p, t); for (int i = 1; i <= t; i++) v -= gpe(p, i) * h[t - i]; h.push_back(v); }
                if (h[e] != mint(0)) self(self, j + 1, d * pe, hd * h[e]);
                if ((i128)pe * p > n) break;
            }
        }
    };
    dfs(dfs, 0, 1, 1);
    return res;
}
ll squarefreeCount(ll n) {                                 // sum_{d <= sqrt n} mu(d) floor(n / d^2)
    ll r = sqrtl(n); while (r * r > n) r--; while ((r + 1) * (r + 1) <= n) r++;
    vector<int> mu(r + 1, 1), mark(r + 1); ll res = 0;
    for (ll i = 2; i <= r; i++) if (!mark[i]) for (ll j = i; j <= r; j += i) { mark[j] = 1; mu[j] = (j / i) % i == 0 ? 0 : -mu[j]; }
    for (ll d = 1; d <= r; d++) res += mu[d] * (n / (d * d));
    return res;
}

Proof: powerful numbers, and the squarefree formula

$h=f\ast g^{-1}$ is multiplicative; at a prime $f(p)=g(p)h(1)+g(1)h(p)=g(p)+h(p)$, so $h(p)=0$ and $h$ is zero on any number with an exponent 1. Then ${\sum }_{i\le n}f(i)={\sum }_{d}h(d){\sum }_{e\le n/d}g(e)={\sum }_{d}h(d)G(\lfloor n/d\rfloor )$. A powerful number is $a^2{b}^3$ with $b$ squarefree, so there are at most ${\sum }_b\sqrt{n/b^3}=O(\sqrt{n})$ of them. The DFS assigns each powerful $d$ its prime factorisation with exponents $\ge 2$ exactly once. Squarefree: $[i\text{ squarefree}]={\sum }_{d^2\mid i}\mu (d)$ (check on prime powers), and summing over $i\le n$ gives ${\sum }_d\mu (d)\lfloor n/d^2\rfloor$; the inner mu sieve marks $\mu (j)=0$ when $i^2\mid j$ and flips the sign per prime factor.