介绍

线性筛不仅能用 $O(n)$ 求出每一个数是否为质数，它还能求出每一个数的因数数量，因数和，甚至更加general的除数函数，莫比乌斯函数等等。

总结来说，线性筛可以在 $O(n)$ 中求出一些积性函数的值。

原理

线性筛之所以是 $O(n)$，是因为对于每一个数 $x$，我们只考虑它 最小 的那个质因数。

在当前使用 $i$ 来筛的时候，如果 i % primes[j] == 0，说明 $i$ 中包含了 primes[j]，那么再往后的话筛掉 cur = i * primes[j+1] 时用到的最小质因数就是 primes[j] 了，所以可以直接break。

const int M = 1e7;
bool isPrime[M+5];
int small[M+5];
vector<int> primes;
void preprocess() {   // 线性筛
    memset(isPrime, 1, sizeof(isPrime));
    small[1] = 1;
    for (int i = 2; i <= M; i++) {
        if (isPrime[i]) primes.push_back(i), small[i] = i;
        for (int j = 0; j < primes.size() && i * primes[j] <= M; j++) {
            int cur = i * primes[j];
            isPrime[cur] = 0;
            small[cur] = primes[j];   // 最小的质因子
            if (i % primes[j] == 0) break;
        }
    }
}

那么观察上述过程，其实我们发现了：

  1. 当 i % primes[j] == 0 时，说明 cur = i * primes[j] 里面至少有 $2$ 个 $p_j$（后文中 $p_j$ 代指 primes[j]），并且 cur 与 i 的最小质因子均为 $p_j$。
  2. 当 i % primes[j] != 0 时，说明 cur = i * primes[j] 里面只有 $1$ 个 $p_j$。

这个重要性质在后面我们求欧拉函数，因数个数时都有重大用处。

• 简单证明一下第一条中的 cur 与 $i$ 的最小质因子均为 $p_j$：首先 cur 肯定满足，因为线性筛的本质就是只用最小质因子筛掉一个数。对于 $i$ 来说，$i = \frac{cur}{p_j}$，说明 $i$ 只是比 cur 少了一个 $p_j$ 而已，所以 $i$ 的最小质因子也是 $p_j$。

应用

欧拉函数

欧拉函数的两个性质：

  1. 积性函数：$\varphi(mn) = \varphi(m)\varphi(n) \iff \gcd(m,n) = 1$
  2. 质因数：$\forall n = p_1^{k_1}p_2^{k_2}…p_r^{k_r}$, 如果 $~\exists ~i, ~s.t. ~k_i > 1$, 则 $\varphi(n) = \varphi(\frac{n}{p_i})*p_i$

运用上面的性质就可以用线性筛了。

代码 
const int M = 1e7;
int phi[M+5];
bool isPrime[M+5];
vector<int> primes;
void euler_function() {
    phi[1] = 1;
    memset(isPrime, 1, sizeof(isPrime));
    for (int i = 2; i <= M; i++) {
        if (isPrime[i]) {
            phi[i] = i-1;
            primes.push_back(i);
        }
        for (int j = 0; j < primes.size() && i * primes[j] <= M; j++) {
            int cur = i * primes[j];
            isPrime[cur] = 0;
            if (i % primes[j] == 0) {
                phi[cur] = phi[i] * primes[j];
                break;
            } else {
                phi[cur] = phi[i] * phi[primes[j]];
            }
        }
    }
}

因数个数

用 $d_i$ 表示 $i$ 的因数个数（有多少个数整除它），用 cnt[i] 表示 $i$ 的最小质因数出现次数。

$d_i$ 是积性函数，所以可以用线性筛。

注意到因为 $n=\prod\limits_{i=1}^m p_i^{c_i}$，则

$$d_n = \prod\limits_{i=1}^m (c_i+1) = (1+c_1)(1+c_2)…(1+c_m)$$

所以需要统计 cnt[i] 表示 $i$ 的最小质因数个数。

代码 
const int M = 1e7;
int d[M+5], cnt[M+5];  // d[i]: i 的因数数量，cnt[i]: i 的最小质因数次数
bool isPrime[M+5];
vector<int> primes;
void divisor_count() {
    cnt[1] = d[1] = 1;
    memset(isPrime, 1, sizeof(isPrime));
    for (int i = 2; i <= M; i++) {
        if (isPrime[i]) {
            cnt[i] = 1;
            d[i] = 2;
            primes.push_back(i);
        }
        for (int j = 0; j < primes.size() && i * primes[j] <= M; j++) {
            int cur = i * primes[j];
            isPrime[cur] = 0;
            if (i % primes[j] == 0) {
                cnt[cur] = cnt[i] + 1;
                d[cur] = d[i] / (cnt[i] + 1) * (cnt[i] + 2);  // 已经有的质因数，所以把原来的次数+1除掉，再乘上次数+2
                break;
            } else {
                cnt[cur] = 1;
                d[cur] = d[i] * 2;  // 出现了一个新的质因数，出现次数为 1，所以直接 *2
            }
        }
    }
}

因数和

对于一个数 $n$ 的因数和，我们将它质因数分解后有 $n=\prod\limits_{i=1}^m p_i^{c_i}$，则

$$sum_n = \prod\limits_{i=1}^m (\sum\limits_{j=1}^{c_i} p_j) = (1+p_1+p_1^2+…+p_1^{c_1})(1+p_2+p_2^2+…+p_2^{c_2})…(1+p^m+p_m^2+…+p_m^{c_m})$$

所以与因数个数不同的地方在于，我们需要处理 small_sum[i] 代表：

若 $i$ 的最小质因数为 $p$，出现次数为 $c$，那么

$$smallsum[i] = p^0 + p^1 + … + p^c$$

代码 
const int M = 1e7;
ll sum[M+5], small_sum[M+5];  // sum[i]: i 的因数和
// small_sum[i]: i 的最小质因数为 p，出现次数为 c，那么 small_sum[i] = p^0 + p^1 + ... + p^c
bool isPrime[M+5];
vector<int> primes;
void divisor_sum() {
    sum[1] = small_sum[1] = 1;
    memset(isPrime, 1, sizeof(isPrime));
    for (int i = 2; i <= M; i++) {
        if (isPrime[i]) {
            small_sum[i] = i + 1;
            sum[i] = i + 1;
            primes.push_back(i);
        }
        for (int j = 0; j < primes.size() && i * primes[j] <= M; j++) {
            int cur = i * primes[j];
            isPrime[cur] = 0;
            if (i % primes[j] == 0) {
                sum[cur] = sum[i] / small_sum[i];
                small_sum[cur] = small_sum[i] * primes[j] + 1;
                sum[cur] *= small_sum[cur];
                break;
            } else {
                small_sum[cur] = primes[j] + 1;
                sum[cur] = sum[i] * small_sum[cur];
            }
        }
    }
}

例题

例1 loj124. 除数函数求和 1

题意

定义 $\sigma_k(n) = \sum_{d | n} d ^ k$，给定 $k$，求 $\sum\limits_{i=1}^n \sigma_k(i)$。

其中，$1 \leq n,k \leq 10^7$，答案对 $10^9+7$ 取模。

法一 因数贡献 - 题解

考虑每一个因数 $d$ 对于答案的贡献，很明显贡献为 $d^k * \lfloor \frac{n}{d} \rfloor$。

那么枚举因数即可。按理来说这个复杂度应该不好过，但是居然过了？

• 如果过不了，那么可以预处理每一个质因数的 $p^k$，然后通过分解质因数来计算 $d^k$。

代码 
#include <bits/stdc++.h>
using namespace std;
ll k;
int n;
int main() {
    cin >> n >> k;
    Z ans = 0;
    for (int i = 1; i <= n; i++) {
        ans += qpow((Z)(i), k) * (n / i);
    }
    cout << ans << endl;
}
法二 线性筛 - 题解

参考上面的因数和计算方式，由于 $n=\prod\limits_{i=1}^m p_i^{c_i}$，那么要求出一个数 $n$ 的 $\sigma_k(n) = \sum_{d | n} d ^ k$，我们实际上有：

$$\sigma_k(n) = \sum_{d | n} d ^ k = (1^k + p_1^k + p_1^{2k} + … + p_1^{c_1k})(1^k + p_2^k + p_2^{2k} + … + p_2^{c_2k}) …$$

那么可以用线性筛如法炮制上面的做法，不过这个还是需要预处理每个质因数的 $k$ 次方，不仅复杂度没改善，还难写。

代码 
#include <bits/stdc++.h>
using namespace std;
const int maxn = 1e7+5;

ll k;
const int M = 1e7;
Z sum[M+5], small_sum[M+5];  // sum[i]: i 的所有因数 d^k 之和
// small_sum[i]: i 的最小质因数为 p，出现次数为 c，那么 small_sum[i] = p^0k + p^1k + ... + p^ck
bool isPrime[M+5];
vector<int> primes;
vector<Z> primes_k;  // primes_k[j] = qpow(primes[j], k)
void divisor_sum() {
    sum[1] = small_sum[1] = 1;
    memset(isPrime, 1, sizeof(isPrime));
    for (int i = 2; i <= M; i++) {
        if (isPrime[i]) {
            small_sum[i] = qpow((Z)i, k) + 1;
            sum[i] = small_sum[i];
            primes.push_back(i);
            primes_k.push_back(qpow((Z)i, k));
        }
        for (int j = 0; j < primes.size() && i * primes[j] <= M; j++) {
            int cur = i * primes[j];
            isPrime[cur] = 0;
            if (i % primes[j] == 0) {
                sum[cur] = sum[i] / small_sum[i];
                small_sum[cur] = small_sum[i] * primes_k[j] + 1;
                sum[cur] *= small_sum[cur];
                break;
            } else {
                small_sum[cur] = primes_k[j] + 1;
                sum[cur] = sum[i] * small_sum[cur];
            }
        }
    }
}

int n;
int main() {
    cin >> n >> k;
    divisor_sum();
    Z ans = 0;
    for (int i = 1; i <= n; i++) ans += sum[i];
    cout << ans << endl;
}

例2 loj125. 除数函数求和 2

题意

给定 $n$，求 $\sum_{i = 1} ^ n 2\sigma_2(i)+3\sigma_1(i)+5\sigma_0(i)$，其中 $\sigma_k(n) = \sum_{d | n} d ^ k$。

答案对 998244353 取模。

其中，$1 \leq n \leq 10^9$。

题解

其实就是上一题的加强版，仍然是考虑贡献。

考虑每一个因数 $d$ 对于答案的贡献，很明显贡献为 $d^k * \lfloor \frac{n}{d} \rfloor$。

注意到 $\lfloor \frac{n}{d} \rfloor$ 很明显可以用数论分块。

我们只要能求出在一个块 $[L,R]$ 中，$d^2$ 的和即可。

这个其实就是平方和公式：

$$\sum\limits_{i=1}^n i^2 = 1^2+2^2+…+n^2 = \frac{n(n+1)(2n+1)}{6}$$

平方和公式推导

$$(n+1)^3 - n^3 = 3n^2+3n+1$$

从 $1$ 开始的话有

$$2^3 - 1^3 = 3\times1^2 + 3\times1 + 1$$

那么一直累加，可以得到

$$(n+1)^3 - 1^3 = 3 \sum\limits_{i=1}^n i^2 + 3 \sum\limits_{i=1}^n i + \sum\limits_{i=1}^n 1$$

$$= 3 \sum\limits_{i=1}^n i^2 + 3 \frac{n(n+1)}{2} + n$$

那么移项可以得到

$$3\sum\limits_{i=1}^n i^2 = (n+1)^3 + 3 \frac{n(n+1)}{2} + (n + 1)$$

化简可得到

$$\sum\limits_{i=1}^n i^2 = \frac{n(n+1)(2n+1)}{6}$$

参考链接

代码 
#include <bits/stdc++.h>
using namespace std;
ll n;
Z inv6 = (Z)(1) / 6;

Z sum2(ll x) {
    return (Z)(x) * (x+1) * (2*x+1) * inv6;
}
Z sum2(ll L, ll R) {
    return sum2(R) - sum2(L-1);
}
Z sum1(ll x) {
    return (x * (x+1)) / 2;
}
Z sum1(ll L, ll R) {
    return sum1(R) - sum1(L-1);
}
Z sum0(ll L, ll R) {
    return R - L + 1;
}
int main() {
    cin >> n;
    ll L = 1;
    Z ans = 0;
    while (L <= n) {
        ll R = n / (n / L);
        ans += (2 * sum2(L, R) + 3 * sum1(L, R) + 5 * sum0(L, R)) * (n / L);
        L = R + 1;
    }
    cout << ans << endl;
}