质数的线性筛法和积性函数值的计算

质数的线性筛法比埃氏筛法复杂度更优，且可以被用来求积性函数的函数值。本文给出了相应的数学背景和代码实现，给出了三个常用积性函数（欧拉函数、莫比乌斯函数和除数函数）的线性筛法求法。

数学背景

对任意大于 $1$ 的自然数 $n$ 进行质因数分解： $n = \prod_{i} p_{i}^{e_{i}}$ 均存在最小的 $p_{i}$ ，设为 $p_{1}$ 。则：

$p_{1}$ 是 $n$ 最小的质因数，也是 $n$ 最小的除 $1$ 以外的因数（最小非平凡因数）；
$n / p_{1}$ 是 $n$ 最大的除 $n$ 以外的因数（最大非平凡因数）；
当且仅当 $p_{1} = n$ ，或者 $n / p_{1} = 1$ ， $n$ 是质数。
$n / p_{1}$ 的最小的质因数不小于 $p_{1}$ ，当且仅当 $p_{1}$ 在 $n$ 的质因数分解中对应的次数 $e_{1}$ 满足 $e_{1} \geq 2$ 时 $n / p_{1}$ 的最小的质因数等于 $p_{1}$ 。这个结论的证明可以考虑将 $n$ 的所有质因数放在一个可重集合 $S$ 中，该可重集合的最小元素是 $p_{1}$ ；而 $n / p_{1}$ 的所有质因数构成的可重集合 $S^{'}$ 是 $S$ 去掉一个 $p_{1}$ ，则显然 $S^{'}$ 的最小元素不小于 $p_{1}$ 。当 $n / p_{1}$ 的最小质因数小于 $p_{1}$ 时，显然 $gcd (n / p_{1}, p_{1}) = 1$ 。

质数的线性筛

现在我们转换一下视角，把重点放在上面的 $n / p_{1}$ ，即 $n$ 的最大非平凡因数上。对于每个大于 $1$ 的自然数 $i$ ，枚举所有的自然数 $n$ ，使得 $i$ 是 $n$ 的最大非平凡因数。这需要 $i$ 乘上一个质数 $p_{j}$ ，其中 $p_{j}$ 扮演的角色是 $n$ 的最小质因数，来实现。但这个质数 $p_{j}$ 是有范围限制的，因为根据上面的讨论， $p_{j}$ 不应当超过 $i$ 的最小质因数。由于 $i$ 大于 $1$ ，被枚举到的 $n$ 显然是合数。而且由于这种枚举方式的唯一性（最大非平凡因数乘以最小质因数），每个合数 $n$ 只会被枚举到 $1$ 次。这样的筛法是最优的线性时间复杂度，我们称之为线性筛。在上面的算法中，在遍历 $i$ 的过程中，我们需要维护一个不超过 $i$ 的质数列表，这样才能枚举 $p_{j}$ 。（因为 $p_{j}$ 不超过 $i$ 的最小质因数，显然有 $p_{j} \leq i$ ，当且仅当 $i$ 是质数取等号。）

重新整理一下思路：枚举所有的 $i \geq 2$ 。如果 $i$ 此前没有被标记过，则它是一个质数，将其加入到我们维护的质数列表中。然后从小到大遍历质数列表中的所有质数 $p_{j}$ ，将 $n = i \times p_{j}$ 标记为合数。当 $p_{j} | i$ 时， $p_{j}$ 是 $i$ 的最小质因数，此时应标记最后一个 $n$ 并停止循环。

代码实现如下：

const int N = 1e7+5;
short flg[N]; int p[N], tot;
void seive(int n) {
    rep(i,2,n) {
        if (!flg[i]) p[++tot]=i;
        for(int j=1; j<=tot && i*p[j]<=n; ++j) {
            flg[i * p[j]] = 1;
            if (i % p[j] == 0) break;
        }
    }
}

积性函数值的计算

与线性筛的过程契合，在线性筛中求积性函数值一般需要考虑函数在以下三种情况的值，其中后两种考虑 $n$ 和 $n / p_{1}$ 的递推关系。 $f (n) = {\begin{cases} Case 1 & n is prime, \\ Case 2 & p_{1} ∤ n / p_{1}, \\ Case 3 & p_{1} ∣ n / p_{1} . \end{cases}$ 实际上第一种情况往往显然；第二种情况可以根据积性函数的性质得到（由于互质关系） $f (n) = f (p_{1}) f (n / p_{1})$ ，只有第三种情况需要额外考虑。

欧拉函数的计算

已经知道 $ϕ (n) = \prod_{i} p_{i}^{e_{i}} (1 - \frac{1}{p_{i}}) = n \prod_{i} (1 - \frac{1}{p_{i}})$ 当 $n$ 是质数时， $ϕ (n) = n - 1$ 当 $p_{1} ∤ n / p_{1}$ 时， $ϕ (n) = ϕ (n / p_{1}) (p_{1} - 1)$ 当 $p_{1} ∣ n / p_{1}$ 时， $ϕ (n) = ϕ (n / p_{1}) p_{1}$ 代码实现：

const int N = 1e7+5;
int phi[N]; short flg[N]; int p[N], tot;
void getPhi(int n) {
    phi[1] = 1;
    rep(i,2,n) {
        if (!flg[i]) phi[i] = i-1, p[++tot]=i;
        for(int j=1; j<=tot && i*p[j]<=n; ++j) {
            flg[i * p[j]] = 1;
            if (i % p[j] == 0) {
                phi[i * p[j]] = p[j] * phi[i];
                break;
            } else {
                phi[i * p[j]] = phi[p[j]] * phi[i];
            }
        }
    }
}

莫比乌斯函数的计算

已经知道 $μ (n) = {\begin{cases} 1 & n = 1, \\ (- 1)^{k} & n is a product of k distinct primes, \\ 0 & otherwise . \end{cases}$ 当 $n$ 是质数时， $μ (n) = - 1$ 当 $p_{1} ∤ n / p_{1}$ 时， $μ (n) = - μ (n / p_{1})$ 当 $p_{1} ∣ n / p_{1}$ 时， $μ (n) = 0$ 代码实现如下：

const int N = 1e7+5;
short mu[N]; short flg[N]; int p[N], tot;
void getMu(int n) {
    mu[1] = 1;
    rep(i,2,n) {
        if (!flg[i]) mu[i] = -1, p[++tot]=i;
        for(int j=1; j<=tot && i*p[j]<=n; ++j) {
            flg[i * p[j]] = 1;
            if (i % p[j] == 0) {
                mu[i * p[j]] = 0;
                break;
            } else {
                mu[i * p[j]] = -mu[i];
            }
        }
    }
}

除数函数的计算

如果是计算除数的和：

已经知道 $σ (n) = \prod_{i} (1 + p_{i} + p_{i}^{2} + \dots + p_{i}^{e_{i}})$ 记 $g (n) = 1 + p_{1} + p_{1}^{2} + \dots + p_{1}^{e_{1}}$ 当 $n$ 是质数时， $g (n) = 1 + n σ (n) = 1 + n$ 当 $p_{1} ∤ n / p_{1}$ 时，即 $e_{1} = 1$ ， $g (n) = 1 + p_{1} σ (n) = σ (p_{1}) σ (n / p_{1})$ 第二式由 $σ$ 的积性可得。当 $p_{1} ∣ n / p_{1}$ 时，即 $e_{1} \geq 2$ ， $g (n) = g (n / p_{1}) p_{1} + 1 σ (n) = \frac{σ (n / p_{1}) g (n)}{g (n / p_{1})}$

const int N = 1e7+5;
int sigma[N], g[N]; short flg[N]; int p[N], tot;
void getMu(int n) {
    sigma[1] = 1;
    rep(i,2,n) {
        if (!flg[i]) sigma[i]=g[i]=i+1, p[++tot]=i;
        for(int j=1; j<=tot && i*p[j]<=n; ++j) {
            flg[i * p[j]] = 1;
            if (i % p[j] == 0) {
                g[i * p[j]] = g[i] * p[j] + 1;
                sigma[i * p[j]] = sigma[i] / g[i] * g[i * p[j]];
                break;
            } else {
                g[i * p[j]] = p[j]+1;
                sigma[i * p[j]] = sigma[i] * sigma[p[j]];
            }
        }
    }
}

如果计算除数的个数，思路类似：

const int N = 1e7 + 5;
short flg[N]; int p[N], tot;
int d[N]; int g[N];
void init(int n) {
    d[1] = 1;
    rep(i,2,n) {
        if (!flg[i]) d[i] = g[i] = 2, p[++tot] = i;
        for(int j=1; j<=tot && i*p[j]<=n; j++) {
            flg[i * p[j]] = 1;
            if (i % p[j] == 0) {
                g[i * p[j]] = g[i] + 1;
                d[i * p[j]] = d[i] / (g[i]) * (g[i * p[j]]);
                break;
            } else {
                g[i * p[j]] = 2;
                d[i * p[j]] = d[i] * d[p[j]];
            }
        }
    }
}