质数的线性筛法和积性函数值的计算 | Eric Way's Personal Site

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

数学背景

对任意大于$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 \ge 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 \le i$，当且仅当$i$是质数取等号。）

重新整理一下思路：枚举所有的$i \ge 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} \text{Case 1} & \text{$n$ is prime,} \\ \text{Case 2} & p_1 \nmid n /p_1, \\ \text{Case 3} & p_1 \mid n/p_1. \end{cases} $$ 实际上第一种情况往往显然；第二种情况可以根据积性函数的性质得到（由于互质关系）$f(n)=f(p_1)f(n/p_1)$，只有第三种情况需要额外考虑。

欧拉函数的计算

已经知道 $$ \phi(n) = \prod_i p_i^{e_i}(1-\frac 1 {p_i}) = n\prod_i(1-\frac {1}{p_i}) $$ 当$n$是质数时， $$ \phi(n)=n-1 $$ 当$ p_1 \nmid n /p_1$时， $$ \phi(n) = \phi(n/p_1) (p_1-1) $$ 当$ p_1 \mid n /p_1$时， $$ \phi(n)=\phi(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];
            }
        }
    }
}

莫比乌斯函数的计算

已经知道 $$ \mu(n) = \begin{cases} 1 & n=1,\\ (-1)^{k} & n \text{ is a product of } k \text{ distinct primes}, \\ 0 & \text{otherwise}.\end{cases} $$ 当$n$是质数时， $$ \mu(n)=-1 $$ 当$ p_1 \nmid n /p_1$时， $$ \mu(n) = -\mu(n/p_1) $$ 当$ p_1 \mid n /p_1$时， $$ \mu(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];
            }
        }
    }
}

除数函数的计算

如果是计算除数的和：

已经知道 $$ \sigma(n) = \prod_i \left( 1+p_i+p_i^2+\dots+p_i^{e_i} \right) $$ 记 $$ g(n)=1+p_1+p_1^2+\dots+p_1^{e_{1}} $$ 当$n$是质数时， $$ g(n)=1+n \\ \sigma(n)=1+n $$ 当$ p_1 \nmid n /p_1$时，即$e_1=1$， $$ g(n)=1+p_1 \\ \sigma(n)=\sigma(p_1)\sigma(n/p_1) $$ 第二式由$\sigma$的积性可得。 当$ p_1 \mid n /p_1$时，即$e_1 \ge 2$， $$ g(n)=g(n/p_1)p_1+1 \\ \sigma(n)=\frac{\sigma(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]];
            }
        }
    }
}