筛法

· · 算法·理论

本文原载于作者的博客园。

前言

本文介绍了下述筛法：

一个自然的想法是对于每个 x，都进行一次素性检验。
显然，这一暴力的做法无法通过较大的数据范围。

于是，素数筛法就出现了。

埃拉托斯特尼筛（埃筛）

我们发现，对于一个正整数 x>1，它除本身之外的倍数必定是合数。
利用这一结论，我们可以避免很多不必要的检测。

若我们从小到大考虑每个数 x，同时将它除本身之外的倍数 y 都标记为合数，那么最后没被标记的就是质数了。

显然，我们只要考虑所有的 2\le x\le\lfloor\sqrt{n}\rfloor 即可。

实现

int P[N],cnt;
bool np[N];

void Eratosthenes(int n){
    cnt=0;
    np[0]=np[1]=true;
    for(int x=2;x*x<=n;++x)if(!np[x])for(int y=x*x;y<=n;y+=x)np[y]=true;
    for(int x=2;x<=n;++x)if(!np[x])P[++cnt]=x;
}

可以证明，埃筛的时间复杂度为 O(n\log\log n)。

一些优化

只筛奇数

由于除了 2 的偶数都是合数，因此我们只需考虑奇数即可。

这不仅能让我们的内存需求减半，还能让所需的操作减半。

压位

我们使用的 np 数组是 bool 类型的。一般一个 bool 类型变量占 1 字节（即 8 比特），但储存一个 bool 值只需 1 比特。
利用位运算相关知识，我们可以将每个 bool 值压到一个比特位里，这样我们仅需 n 比特（即 \displaystyle\frac{n}{8} 字节）而并非 n 字节，可以显著减少内存占用。这种方式被称为压位。

在 STL 中，vector<bool> 和 bitset<> 都会进行压位，并且有常数优化（后者稍微好一点）。
若觉得常数还是过大，可以手写 bitset。不过这就不是本篇文章要讨论的东西了。

线性筛（欧拉筛）

注意到在埃筛中会存在一个合数被多次标记的情况。
若能使每个合数只被标记一遍，就可以把时间复杂度改进为 O(n)。

int P[N],cnt;
bool np[N];

void Euler(int n){
    cnt=0;
    np[0]=np[1]=true;
    for(int x=2;x<=n;++x){
        if(!np[x])P[++cnt]=x;
        for(int i=1;i<=cnt&&P[i]*x<=n;++i){
            np[P[i]*x]=true;
            if(x%P[i]==0)break;
        }
    }
}

解释一下为什么每个合数只会被标记一次：

当 x\bmod P_i=0 时，换言之，x 之前被 P_i 筛过了。
由于 P 保存的质数是从小到大的，因此 x 乘上其他的质数的结果在后面一定会被 P_i 的倍数筛掉。
因此，这里不需要先筛一次，直接 break 即可。

线性筛求积性函数

在一类问题中，我们要预处理对于 \forall1\le x\le n（n 为给定的正整数），一些积性函数 f(x) 的值。
如：欧拉函数 \varphi(x)，莫比乌斯函数 \mu(x)，除数函数 \sigma_k(x)。

接下来将一一介绍如何用线性筛求出它们。

在以下的讲述中，皆有 \displaystyle x=\prod_{i=1}^qp_i^{\alpha_i}，p_1 为 x 最小的素因子，x'=\dfrac{x}{p_1}。
当提到 x\not\in\mathbb P 时，默认 x\neq1。

欧拉函数 \varphi(x)

首先，有 \varphi(1)=1。

注意到在线性筛中，每一个合数 x 都是被 p_1 筛掉的。
显然，在线性筛的过程中 x 是通过 p_1\times x' 筛掉的。

下面对 x'\bmod p_1 进行分类讨论：

\varphi(x)&=x\times\prod_{i=1}^q\frac{p_i-1}{p_i}\\ &=p_1\times x'\times\prod_{i=1}^q\frac{p_i-1}{p_i}\\ &=p_1\times\varphi(x') \end{aligned}

当 x\in\mathbb P 时，则有 \varphi(x)=x-1。

实现

int P[N],cnt;
bool np[N];
int phi[N];

void CalcPhi(int n){
    cnt=0;
    np[0]=np[1]=true;
    phi[1]=1;
    for(int x_=2;x_<=n;++x_){
        if(!np[x_]){
            P[++cnt]=x_;
            phi[x_]=x_-1;
        }
        for(int i=1;i<=cnt&&P[i]*x_<=n;++i){
            int x=P[i]*x_;
            np[x]=true;
            if(x_%P[i]==0){
                phi[x]=P[i]*phi[x_];
                break;
            }
            phi[x]=phi[P[i]]*phi[x_];
        }
    }
}

莫比乌斯函数 \mu(x)

当 x=1：\mu(x)=1。
当 x\in\mathbb P 时：\mu(x)=-1。
当 x\not\in\mathbb P 时：

void CalcMu(int n){ cnt=0; np[0]=np[1]=true; mu[1]=1; for(int x=2;x<=n;++x){ if(!np[x]){ P[++cnt]=x; mu[x]=-1; } for(int i=1;i<=cnt&&P[i]x_<=n;++i){ int x=P[i]x; np[x]=true; if(x%P[i]==0)break; mu[x]=-mu[x_]; } } }

## 除数函数 $\sigma_k(x)$
令 $f(x)=p_1^{\alpha_1}$。

当 $x=1$：$\sigma_k(x)=1$。  
当 $x\in\mathbb P$：$f(x)=x,\sigma_k(x)=x^k+1$。  
当 $x\not\in\mathbb P$：
- $x'\bmod p_1=0$：
$f(x)=p_1\times f(x'),\sigma_k(x)=\begin{cases}\sigma_k(f(x))\times\sigma_k(\dfrac{x}{f(x)}),&f(x)\neq x\\ p_1^k\times\sigma_k(x') +1,&f(x)=x\end{cases}$。
- $x'\bmod p_1\neq0$：
$f(x)=p_1,\sigma_k(x)=\sigma_k(x')\times\sigma_k(p_1)$。
### 实现
```cpp
int QuickPow(int x,int y){
    int ans=1;
    while(y){
        if(y&1)ans=ans*x;
        x=x*x;
        y>>=1;
    }
    return ans;
}

int P[N],cnt;
bool np[N];
int f[N],Pow[N],sigma[N];

void CalcSigma(int n,int k){
    cnt=0;
    np[0]=np[1]=true;
    sigma[1]=1;
    for(int x_=2;x_<=n;++x_){
        if(!np[x_]){
            ++cnt;
            P[cnt]=x_;
            Pow[cnt]=QuickPow(x_,k);
            f[x_]=x_;
            sigma[x_]=Pow[cnt]+1;
        }
        for(int i=1;i<=cnt&&P[i]*x_<=n;++i){
            int x=P[i]*x_;
            np[x]=true;
            if(x_%P[i]==0){
                f[x]=P[i]*f[x_];
                if(f[x]!=x)sigma[x]=sigma[f[x]]*sigma[x/f[x]];
                else sigma[x]=Pow[i]*sigma[x_]+1;
                break;
            }
            f[x]=P[i];
            sigma[x]=sigma[P[i]]*sigma[x_];
        }
    }
}

一般的积性函数 f(x)

和求 \sigma_k(x) 类似，现在将过程概括如下。

令 g(x)=p_1^{\alpha_1}。

当 x=1：显然有 f(x)=1。
当 x\in\mathbb P：g(x)=x，f(x) 根据函数的性质 O(1) 计算。
当 x\not\in\mathbb P：

对于较小的 n，我们可以使用线性筛预处理出来；
若 n 的最大值为 N，则可以筛出 N^{\frac{2}{3}} 内的值，此时时间复杂度为 O(n^{\frac{2}{3}})。

为了进一步加快计算，可以使用记忆化：将计算过的值用 map/unordered_map 存下来。

实现

int P[N],cnt;
bool np[N];
int phi[N],SumOfPhi[N];

void EulerSieveSumOfPhi(int n){
    cnt=0;
    np[0]=np[1]=true;
    phi[1]=1;
    for(int x_=2;x_<=n;++x_){
        if(!np[x_]){
            P[++cnt]=x_;
            phi[x_]=x_-1;
        }
        for(int i=1;i<=cnt&&P[i]*x_<=n;++i){
            int x=P[i]*x_;
            np[x]=true;
            if(x_%P[i]==0){
                phi[x]=P[i]*phi[x_];
                break;
            }
            phi[x]=phi[P[i]]*phi[x_];
        }
    }
    for(int x=1;x<=n;++x)SumOfPhi[x]=SumOfPhi[x-1]+phi[x];
}

unordered_map<int,int> MapSumOfPhi;

int CalcSumOfPhi(int n){
    if(n<=MAXN)return SumOfPhi[n];
    if(MapSumOfPhi[n])return MapSumOfPhi[n];
    int ans=n*(n+1)/2;
    for(int l=2,r;l<=n;l=r+1){
        r=n/(n/l);
        ans-=(r-l+1)*CalcSumOfPhi(n/l);
    }
    return MapSumOfPhi[n]=ans;
}

莫比乌斯函数 \mu(n)

原理与欧拉函数 \varphi(n) 类似，只不过利用的是 \epsilon=\mu*1：

1 &=\sum_{i=1}^n\epsilon(i)\\ &=\sum_{i=1}^n\sum_{d\mid i}\mu(\dfrac{n}{d})\\ &=\sum_{d=1}^n\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\mu(i)\\ &=\sum_{d=1}^nS(\lfloor\dfrac{n}{d}\rfloor)\\ &=S(n)+\sum_{d=2}^nS(\lfloor\dfrac{n}{d}\rfloor) \end{aligned}

因此 S(n)=1-\sum_{d=2}^nS(\lfloor\dfrac{n}{d}\rfloor)。

实现

int P[N],cnt;
bool np[N];
int mu[N],SumOfMu[N];

void EulerSieveSumOfMu(int n){
    cnt=0;
    np[0]=np[1]=true;
    mu[1]=1;
    for(int x_=2;x_<=n;++x_){
        if(!np[x_]){
            P[++cnt]=x_;
            mu[x_]=-1;
        }
        for(int i=1;i<=cnt&&P[i]*x_<=n;++i){
            int x=P[i]*x_;
            np[x]=true;
            if(x_%P[i]==0)break;
            mu[x]=-mu[x_];
        }
    }
    for(int x=1;x<=n;++x)SumOfMu[x]=SumOfMu[x-1]+mu[x];
}

unordered_map<int,int> MapSumOfMu;

int CalcSumOfMu(int n){
    if(n<=MAXN)return SumOfMu[n];
    if(MapSumOfMu[n])return MapSumOfMu[n];
    int ans=1;
    for(int l=2,r;l<=n;l=r+1){
        r=n/(n/l);
        ans-=(r-l+1)*CalcSumOfMu(n/l);
    }
    return MapSumOfMu[n]=ans;
}

一般的数论函数 f(n)

找到一个合适的数论函数 g(n)，令 h=f*g，则有：

\sum_{i=1}^nh(i) &=\sum_{i=1}^n\sum_{d\mid i}f(\dfrac{n}{d})g(d)\\ &=\sum_{d=1}^n\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}g(d)f(i)\\ &=\sum_{d=1}^ng(d)S(\lfloor\dfrac{n}{d}\rfloor)\\ &=g(1)S(n)+\sum_{d=2}^ng(d)S(\lfloor\dfrac{n}{d}\rfloor) \end{aligned}

因此 \displaystyle S(n)=\dfrac{\sum_{i=1}^nh(i)-\sum_{d=2}^ng(d)S(\lfloor\dfrac{n}{d}\rfloor)}{g(1)}。

这就要求：