筛法

TKXZ133's Blog / 2023-06-05 / 原文

埃氏筛

如何求出 $1\sim n$ 中有几个质数？

考虑这样一件事情：对于任意一个大于 $1$ 的正整数 $n$，那么它的 $x$ 倍就是合数 $(x > 1)$。利用这个结论，我们可以避免很多次不必要的检测。

如果我们从小到大考虑每个数，然后同时把当前这个数的所有（比自己大的）倍数记为合数，那么运行结束的时候没有被标记的数就是素数了。

void Eratosthenes(){
    for(int i=2;i<=n;i++)
        if(!v[i]){
            prime[++cnt]=i;
            if(i*i>n) continue;
            for(int j=i*i;j<=n;j+=i)
                v[j]=1;
        }
}

其时间复杂度为 $O(n\log\log n)$。

如何优化？

筛至平方根

显然，要找到直到 $n$ 为止的所有素数，仅对不超过 $\sqrt n$ 的素数进行筛选就足够了。

只筛奇数

因为除 $2$ 以外的偶数都是合数，所以我们可以直接跳过它们，只用关心奇数就好。

首先，这样做能让我们内存需求减半；其次，所需的操作大约也减半。

减少内存的占用

我们注意到筛法只需要 $n$ 比特的内存。因此我们可以通过将变量声明为布尔类型，只申请 $n$ 比特而不是 $n$ 字节的内存，来显著的减少内存占用。即仅占用 $\frac{n}{8}$ 字节的内存。

但是，这种称为位级压缩的方法会使这些位的操作复杂化。任何位上的读写操作都需要多次算术运算，最终会使算法变慢。

因此，这种方法只有在 $n$ 特别大，以至于我们不能再分配内存时才合理。在这种情况下，我们将牺牲效率，通过显著降低算法速度以节省内存（减小 $8$ 倍）。

当然，可以用 vector<bool> 和 bitset<>。

分块筛选

由优化「筛至平方根」可知，不需要一直保留整个 is_prime[1...n] 数组。为了进行筛选，只保留到 $\sqrt n$ 的素数就足够了，即 prime[1...sqrt(n)]。并将整个范围分成块，每个块分别进行筛选。这样，我们就不必同时在内存中保留多个块，而且 CPU 可以更好地处理缓存。

设 $s$ 是一个常数，它决定了块的大小，那么我们就有了 $\lceil {\frac n s} \rceil$ 个块，而块 $k(k = 0 \dots \lfloor {\frac n s} \rfloor)$ 包含了区间 $[ks, ks + s - 1]$ 中的数字。我们可以依次处理块，也就是说，对于每个块 $k$，我们将遍历所有质数（从 $1$ 到 $\sqrt n$）并使用它们进行筛选。

值得注意的是，我们在处理第一个数字时需要稍微修改一下策略：首先，应保留 $[1, \sqrt n] $ 中的所有的质数；第二，数字 $0$ 和 $1$ 应该标记为非素数。在处理最后一个块时，不应该忘记最后一个数字 $n$ 并不一定位于块的末尾。

int Eratosthenes(int n){
    const int S=10000;
    int nsqrt=sqrt(n);
    vector <int> prime;
    vector <bool> v(nsqrt+1,1);
    for(int i=2;i<=nsqrt;i++)
        if(v[i]){
            prime.push_back(i);
            for(int j=i*i;j<=nsqrt;j++) v[j]=0;
        }
    int res=0;
    vector <bool> block(S);
    for(int k=0;k*S<=n;k++){
        fill(block.begin(),block.end(),1);
        int start=k*S;
        for(int p:prime){
            int started=(start+p-1)/p;
            int j=max(started,p)*p-start;
            for(;j<S;j+=p) block[j]=0;
        }
        if(k==0) block[0]=block[1]=0;
        for(int i=0;i<S&&start+i<=n;i++)
            if(block[i]) res++;
    }
    return res;
}

分块筛法的渐进时间复杂度与埃氏筛法是一样的（除非块非常小），但是所需的内存将缩小为 $O(\sqrt{n} + S)$，并且有更好的缓存结果。另一方面，对于每一对块和区间 $[1, \sqrt{n}]$ 中的素数都要进行除法，而对于较小的块来说，这种情况要糟糕得多。因此，在选择常数 $S$ 时要保持平衡。

块大小 $S$ 取 $10^4 \sim 10^5$ 之间，可以获得最佳的速度。

线性筛

埃氏筛法仍有优化空间，它会将一个合数重复多次标记。有没有什么办法省掉无意义的步骤呢？答案是肯定的。

如果能让每个合数都只被标记一次，那么时间复杂度就可以降到 $O(n)$ 了。

void sieve(){
    for(int i=2;i<=n;i++){
        if(!v[i]){prime[++cnt]=i;}
        for(int j=1;j<=cnt&&i*prime[j]<=n;j++){
            v[i*prime[j]]=1;
            if(i%prime[j]==0) break;
        }
    }
}

注意到筛法求素数的同时也得到了每个数的最小质因子。

筛法求欧拉函数

注意到在线性筛中，每一个合数都是被最小的质因子筛掉。比如设 $p_1$ 是 $n$ 的最小质因子，$n' = \frac{n}{p_1}$，那么线性筛的过程中 $n$ 通过 $n' \times p_1$ 筛掉。

观察线性筛的过程，我们还需要处理两个部分，下面对 $n' \bmod p_1$ 分情况讨论。

如果 $n' \bmod p_1 = 0$，那么 $n'$ 包含了 $n$ 的所有质因子。

\[\begin{aligned} \varphi(n) & = n \times \prod_{i = 1}^s{\frac{p_i - 1}{p_i}} \\\\ & = p_1 \times n' \times \prod_{i = 1}^s{\frac{p_i - 1}{p_i}} \\\\ & = p_1 \times \varphi(n') \end{aligned} \]

那如果 $n' \bmod p_1 \neq 0$ 呢，这时 $n'$ 和 $p_1$ 是互质的，根据欧拉函数性质，我们有：

\[\begin{aligned} \varphi(n) & = \varphi(p_1) \times \varphi(n') \\\\ & = (p_1 - 1) \times \varphi(n') \end{aligned} \]

void sieve(){
    phi[1]=1;
    for(int i=2;i<=n;i++){
        if(!v[i]){prime[++cnt]=i;phi[i]=i-1;}
        for(int j=1;j<=cnt&&i*prime[j]<=n;j++){
            v[i*prime[j]]=1;
            if(i%prime[j]==0){
                phi[i*prime[j]]=phi[i]*prime[j];
                break;
            }
            phi[i*prime[j]]=phi[i]*phi[prime[j]];
        }
    }
}

筛法求莫比乌斯函数

根据莫比乌斯函数的定义，设 $n$ 是一个合数，$p_1 $ 是 $ n $ 的最小质因子，$n'=\frac{n}{p_1}$，有：

\[\mu(n)= \begin{cases} 0 & n' \bmod p_1 = 0\\\\ -\mu(n') & \text{otherwise} \end{cases} \]

若 $n$ 是质数，有 $\mu(n)=-1$。

void sieve(){
    mu[1]=1;
    for(int i=2;i<=n;i++){
        if(!v[i]){prime[++cnt]=i;mu[i]=-1;}
        for(int j=1;j<=cnt&&i*prime[j]<=n;j++){
            v[i*prime[j]]=1;
            if(i%prime[j]==0) break;
            mu[i*prime[j]]=-mu[i];
        }
    }
}

筛法求约数个数

用 $d_i$ 表示 $i$ 的约数个数，$num_i$ 表示 $i$ 的最小质因子出现次数。

约数个数定理：

定理：若 $n=\prod_{i=1}^m p_i^{c_i}$ 则 $d_i=\prod_{i=1}^m (c_i+1)$。

证明：我们知道 $p_i^{c_i}$ 的约数有 $p_i^0,p_i1,\dots ,p_i^{c_i} $ 共 $c_i+1$ 个，根据乘法原理，$n$ 的约数个数就是 $\prod_{i=1}^m (c_i+1)$。
实现：

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

void sieve(){
    d[1]=1;
    for(int i=2;i<=n;i++){
        if(!v[i]){prime[++cnt]=i;d[i]=2;num[i]=1;}
        for(int j=1;j<=cnt&&i*prime[j]<=n;j++){
            v[i*prime[j]]=1;
            if(i%prime[j]==0){
                num[i*prime[j]]=num[i]+1;
                d[i*prime[j]]=d[i]/num[i*prime[j]]*(num[i*prime[j]]+1);
                break;
            }
            num[i*prime[j]]=1;
            d[i*prime[j]]=d[i]*2;
        }
    }
}

筛法求约数之和

$f_i$ 表示 $i$ 的约数和，$g_i$ 表示 $i$ 的最小质因子的 $p^0+p^1+p^2+\dots p^k$。

void sieve(){
    f[1]=g[1]=1;
    for(int i=2;i<=n;i++){
        if(!v[i]){prime[++cnt]=i;g[i]=i+1;f[i]=i+1;}
        for(int j=1;j<=cnt&&i*prime[j]<=n;j++){
            v[i*prime[j]]=1;
            if(i%prime[j]==0){
                g[i*prime[j]]=g[i]*prime[j]+1;
                f[i*prime[j]]=f[i]/g[i]*g[i*prime[j]];
                break;
            }
            f[i*prime[j]]=f[i]*f[prime[j]];
            g[i*prime[j]]=prime[j]+1;
        }
    }
}

线性筛条件

假如一个积性函数 $f$ 满足：对于任意质数 $p$ 和正整数 $k$，可以在 $O(1)$ 时间内计算 $f(p^k)$，那么可以在 $O(n)$ 时间内筛出 $f(1),f(2),\dots,f(n)$ 的值。

设合数 $n$ 的质因子分解是 $\prod_{i=1}^k p_i^{\alpha_i}$，其中 $p_1<p_2<\dots<p_k$ 为质数，我们在线性筛中记录 $g_n=p_1^{\alpha_1}$，假如 $n$ 被 $x\cdot p$ 筛掉（$p$ 是质数），那么 $g$ 满足如下递推式：

\[g_n= \begin{cases} g_x\cdot p & x\bmod p=0\\\\ p & \text{otherwise} \end{cases} \]

假如 $n=g_n$，说明 $n$ 就是某个质数的次幂，可以 $O(1)$ 计算 $f(n)$；否则，$f(n)=f(\frac{n}{g_n})\cdot f(g_n)$。

杜教筛

杜教筛被用来处理数论函数的前缀和问题。对于求解一个前缀和，杜教筛可以在低于线性时间的复杂度内求解。

对于数论函数 $f$，要求我们计算 $S(n)=\sum_{i=1}^{n}f(i)$。

我们想办法构造一个 $S(n)$ 关于 $S\left(\left\lfloor\frac{n}{i}\right\rfloor\right)$ 的递推式

对于任意一个数论函数 $g$，必满足

\[\begin{aligned}\sum_{i=1}^{n}\sum_{d \mid i}g(d)f\left(\frac{i}{d}\right)&=\sum_{i=1}^{n}g(i)S\left(\left\lfloor\frac{n}{i}\right\rfloor\right)\\\iff\sum_{i=1}^{n}(f\ast g)(i)&=\sum_{i=1}^{n}g(i)S\left(\left\lfloor\frac{n}{i}\right\rfloor\right)\end{aligned} \]

证明：

$g(d)f\left(\frac{i}{d}\right)$ 就是对所有 $i\leq n$ 的做贡献，因此变换枚举顺序，枚举 $d,\frac{i}{d}$（分别对应新的 $i,j$）

\[\begin{aligned} &\sum_{i=1}^n\sum_{d \mid i}g(d)f\left(\frac{i}{d}\right)\\ =&\sum_{i=1}^n\sum_{j=1}^{\left\lfloor\frac{n}{i}\right\rfloor}g(i)f(j)\\ =&\sum_{i=1}^ng(i)\sum_{j=1}^{\left\lfloor\frac{n}{i}\right\rfloor}f(j)\\ =&\sum_{i=1}^ng(i)S\left(\left\lfloor\frac{n}{i}\right\rfloor\right) \end{aligned} \]

那么可以得到递推式

\[g(1)S(n)=\sum_{i=1}^n(f\ast g)(i)-\sum_{i=2}^ng(i)S\left(\left\lfloor\frac{n}{i}\right\rfloor\right) \]

那么假如我们可以快速对 $\sum_{i=1}^n(f \ast g)(i)$ 求和，并用数论分块求解 $\sum_{i=2}^ng(i)S\left(\left\lfloor\frac{n}{i}\right\rfloor\right)$ 就可以在较短时间内求得 $g(1)S(n)$。

求莫比乌斯函数前缀和

易知，$\epsilon=\mu*1$，即 $\epsilon(n)=\sum_{d|n}\mu(d)$，故构造 $g=1$，则：

\[S(n)=\sum_{i=1}^n\mu(i)=\sum_{i=1}^n\epsilon(i)-\sum_{i=2}^nS(\lfloor\frac{n}{i}\rfloor)=1-\sum_{i=2}^{n}S(\lfloor\frac{n}{i}\rfloor) \]

观察到 $\lfloor \frac n i \rfloor$ 最多只有 $O(\sqrt n)$ 种取值，我们就可以应用整除分块来计算每一项的值了。

直接计算的时间复杂度为 $O(n^{\frac 3 4})$。考虑先线性筛预处理出前 $n^{\frac{2}{3}}$ 项，剩余部分的时间复杂度为 $O(\int_{0}^{n^{\frac{1}{3}}} \sqrt{\frac{n}{x}} ~ dx)=O(n^{\frac 2 3})$。

对于较大的值，需要用 map 存下其对应的值，方便以后使用时直接使用之前计算的结果。

求欧拉函数前缀和

易知，$\text{id}=\varphi*1$，即 $\text{id}(n)=\sum_{d|n}\mu(d)$，故构造 $g=1$，则：

\[S(n)=\sum_{i=1}^n\varphi(i)=\sum_{i=1}^n\text{id}(i)-\sum_{i=2}S(\lfloor\frac{n}{i}\rfloor)=\frac{1}{2}n(n+1)-\sum_{i=2}^nS(\lfloor\frac{n}{i}\rfloor) \]

#include <bits/stdc++.h>
using namespace std;
const int N=2002000,L=2000000;
#define int long long

int prime[N],phi[N],mu[N],v[N];
int cnt,T,n;

map<int,int> mp_phi,mp_mu;

void sieve(){
    phi[1]=mu[1]=1;
    for(int i=2;i<=L;i++){
        if(!v[i]){prime[++cnt]=i;phi[i]=i-1;mu[i]=-1;}
        for(int j=1;j<=cnt&&i*prime[j]<=L;j++){
            v[i*prime[j]]=1;
            if(i%prime[j]==0){
                phi[i*prime[j]]=phi[i]*prime[j];
                break;
            }
            phi[i*prime[j]]=phi[i]*phi[prime[j]];
            mu[i*prime[j]]=-mu[i];
        }
    }
    for(int i=1;i<=L;i++) phi[i]+=phi[i-1];
    for(int i=1;i<=L;i++) mu[i]+=mu[i-1];
}

int get_mu(int n){
    if(n<=L) return mu[n];
    if(mp_mu[n]) return mp_mu[n];
    int res=1;
    for(int l=2,r;l<=n;l=r+1){
        r=n/(n/l);
        res-=get_mu(n/l)*(r-l+1);
    }
    return mp_mu[n]=res;
}

int get_phi(int n){
    if(n<=L) return phi[n];
    if(mp_phi[n]) return mp_phi[n];
    int res=n*(n+1)/2;
    for(int l=2,r;l<=n;l=r+1){
        r=n/(n/l);
        res-=get_phi(n/l)*(r-l+1);
    }
    return mp_phi[n]=res;
}

signed main(){
    sieve();
    scanf("%lld",&T);
    while(T--){
        scanf("%lld",&n);
        cout<<get_phi(n)<<' '<<get_mu(n)<<'\n';
    }
    return 0;
}