Min25 筛学习笔记

Orz min25 Orz gtr

Min25 筛借助一个神奇的类似 DP 数组的手段降低了复杂度

BTW ，这里说一句题外话，每次看到 DP 的时候，总是不由得想起在 ATABC226F 题题解中看到的一句话，感到这句话是引人深思，并且蒟蒻现在也没有参透：（相信大家英文很好）

Before getting down to business, let us state the general rule, not limited to this problem: when trying to reduce the time complexity with DP, it is important to focus on “what aspect can be identified.” It is a very common technique to equate multiple states into one to reduce the computational complexity, enabling us to solve more problems.

让我们进入正题

扩展埃氏筛

学习自与扩展埃氏筛（min_25筛？）玩耍

其实感觉扩展埃氏筛就是 Min25 的简约版

我们知道埃氏筛的思想是说把每个素数的倍数全部筛掉，其累赘之处在于一个合数会被多个素数反复筛到

一个优化思路成就了线性筛，也就是保证每一个合数只会被筛一次

另外一个优化思路是在于有些题我们好像没必要知道具体有哪些素数

例1 loj6235 区间素数个数

给出一个 $n$ ，求出 $1\to n$ 的素数个数， $n\le 10^{11}$

设 $f(x)$ 表示 $1\dots x$ 之间的素数个数

考虑最开始我们只确定了 $1$ 不是素数，那么有 $f(x)=x-1$

假设现在我们处理到了素数 $p$ ，即需要把所有包含 $p$ 的合数筛掉

那么，对于 $f(x)$ 而言， $1\dots x$ 中，有 $\frac xp$ 个数是 $p$ 的倍数，考虑到直接 $f(x)-\frac xp$ 的话可能会算重，所以我们每次应该干掉的是以 $p$ 为最小质因子的合数（线性筛的思想）。

此时的 $f(\frac xp)$ 存的是什么呢？存的是所有的素数和最小质因子 $\ge p$ 的合数

那我们直接 $f(x)-=f(\frac xp)$ 就可以了？不行，我们损失了所有 $< p$ 的素数，那么最后得到了 $f(x)-=f(\frac xp)-f(p-1)$

埃氏筛只需要处理 $\sqrt n$ 内的素数

考虑最后答案是 $f(n)$ ， $f(n)$ 需要用到的是 $f(P)$ 和 $f(\frac nP)$ ，其中 $P$ 表示所有素数。

显然， $\frac np$ 的取值只有 $\sqrt n$ 种，所以我们记 $f2(x)$ 表示 $1\dots\frac nx$ 之间的素数，这样 $f,f2$ 只需要开 $\sqrt n$ 空间，也就是使用 $f(x)$ 记录 $1\dots\sqrt n$ ， $f2(x)$ 记录 $\sqrt n\dots n$

二者互相推求即可，答案 $f2(1)$ ，据说复杂度 $O(\dfrac{n^{\frac 34}}{ln\ n})$

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N = 1e6+5;

int read(){
    int s=0,w=1;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')w=-1;ch=getchar();}
    while(ch>='0'&&ch<='9')s=s*10+ch-'0',ch=getchar();
    return s*w;
}

int x,f1[N],f2[N];

signed main(){
	x=read(),lim=sqrt(x);
    for(int i=1;i<=lim;++i)f1[i]=sum(i),f2[i]=sum(x/i);
    for(int i=2;i<=lim;++i)if(f1[i]^f1[i-1]){// i 是素数 
        for(int j=1;j<=lim/i;++j)f2[j]-=f2[j*i]-f1[i-1];
        for(int j=lim/i+1;j<=x/i/i&&j<=lim;++j)f2[j]-=f1[x/j/i]-f1[i-1];
        //这个 x/j/i 怎么来的可以把 j=lim/i+1 带进去，容易发现和 f2[j] 意义相同 
        for(int j=lim;j>=i*i;--j)f1[j]-=f1[j/i]-f1[i-1];//从大到小更新 f1，防止重复
    }
	printf("%lld\n",f2[1]); 
    return 0;
}

类似题：loj6202 叶氏筛法

Min25 筛

此版块下的所有 $p$ 默认为质数， $p_i$ 表示第 $i$ 个质数，比如 $p_3=5$ 。定义 $lpf(i)$ 为 $i$ 的最小质因数（ $\text{Lowest Prime Factor}$ ）。

假设我们想求 $S(n)=∑_{i≤n}f(i)$

Min25 筛适用于一些积性函数 $f$ ，其中 $f(p)$ 可以用多项式表示，且 $f(p^k)$ 很好求值。

对于埃氏筛，我们是每次将最小质因子为 $p$ 的合数筛去，剩下的就是质数。

我们知道合数最小质因子至多为 $\sqrt n$ ，所以合数的贡献可以通过枚举最小质因子来计算，质数我们则使用另外的方法。

也就是说，我们把 $S(n)$ 拆成这样：

S (n) = \sum_{i \in P} f (i) + \sum_{i \in P, i \leq \sqrt{n}, i^{e} \leq n} f (i^{e}) \sum_{1 \leq x \leq n / i^{e}, l p f (x) > i} f (x)

$S(n)=\sum_{i\in P}f(i)+\sum_{i\in P,i\leq\sqrt n,i^e\leq n}f(i^e)\sum_{1\leq x\leq n/i^e,lpf(x)> i} f(x)$

即是素数的贡献+以素数 $i$ 为最小质因子的合数的贡献

素数贡献

假装现在我们正在解决这道例题： P5325 【模板】Min_25筛 - 洛谷

为了计算素数贡献，Min25 筛先把 $f$ 看成完全积性函数 $g$ ，其中 $f(x),g(x)$ 在 $x\in P$ 处相同。这也就是为什么要求 $f(p)$ 可以用多项式表示的原因：多项式可以被拆成单项式，而对单项式而言计算 $g$ 是容易的。

——GTR

因为 $f(p^k)=p^k(p^k-1)=p^{2k}-p^k$ ，所以我们不妨平方项和一次向单独看成两个多项式处理

令 $g(p)=p^k,g(x)=\Pi f(p_i^a)$ ,显然，当 $x\notin P$ 时， $g(x)\neq f(x)$

所谓“单项式而言，计算 $g$ 是容易的，也就是一点小学数学知识： $\sum_{i=1}^ni=\frac{n(n+1)}2,\sum_{i=1}^ni^2=\frac{n(n+1)(2n+1)}6,\sum_{i=1}^ni^3=(\frac{n(n+1)}2)^2$

开始计算：

我们可以从前面扩展埃氏筛那里毛一个“DP”数组过来：

设 $G(n,j)$ 表示前 $n$ 个数，经过前 $j$ 个质数的筛选，剩下的数的贡献

标准化的， $G(n,j)=\sum_{i\leq n}g(i)[i\in P \or lpf(i)>p_j]$

所谓剩下的数，就是指的是 $lpf>j$ 的合数和素数，大佬们可以自行理解一下这个 $G$ 和扩展埃氏筛里面那个 $f$ 的等价关系

显然 $G(n,0)=\sum_{i\leq n}g(i)$

显然素数的贡献就是 $G(n,|P|)$ ，我们考虑递推这个 $g$ ，发现 $G(n,j-1)$ 比 $G(n,j)$ 多的是

G (n, j - 1) - G (n, j) = \sum_{i \leq n} g (i) [i \notin P \land l p f (i) = p_{j}]

$G(n,j-1)-G(n,j)=\sum_{i\leq n}g(i)[i\notin P\and lpf(i)=p_j]$

考虑从小到大递推，辣么有：

G (n, j - 1) - G (n, j) = {\begin{cases} 0, & p_{j}^{2} > n \\ g (p_{j}) (G (\frac{n}{p_{j}}, j - 1) - G (p_{j}, j - 1)), & p_{j}^{2} \leq n \end{cases}

$G(n,j-1)-G(n,j)=\begin{cases} 0,&p_j^2>n\\ g(p_j)(G(\frac n{p_j},j-1)-G(p_j,j-1)),&p_j^2\leq n \end{cases}$

怎么理解呢，埃氏筛的式子叫做： $f(x)-=f(\frac xp)-f(p-1)$

而我们这里仅仅是多加了一维（反正埃氏筛定义的 $f$ 其实也是滚动的），而且多算上一个 $g(p_j)$ 的贡献

即是，提出一个最小质因子 $p_j$ 后，那么在剩下的数中就不能有小于它的质因子了，也就是 $G(\frac n{p_j},j-1)$ ，但是这样多减掉了一些素数，哪些？小于 $p_j$ 的所有素数的贡献都没了

所以可以发现 $G(p_j,j-1)=\sum_{i\leq j}g(p_i)$

显然，我们没有必要求出每一个 $G(i,x)$ ，有一个结论叫做： $\lfloor\dfrac{\lfloor\frac na\rfloor}{b}\rfloor=\lfloor\dfrac n{ab}\rfloor$

所以我们只需要求出 $\lfloor\frac nx\rfloor$ 形式的数就可以了，这样的数有 $\sqrt n$ 个，仿照埃氏筛 $f2$ 式的处理，我们如此离散化:

可以用 $ind1[x]$ 表示 $x$ 这个数对应的数组下标， $ind2[x]$ 表示 $n/x$ 这个数对应的下标。这样两个 $ind$ 数组最大都只会到 $\sqrt n$ 。

复杂度 ~~O(能过)~~ $O(\dfrac{n^{3/4}}{\log n})$ ，证明可以看这里

总贡献

注意这里我们把素数和合数贡献揉在一起算了

rua，不同地，我们记 $F(n,j)$ 表示 $\sum_{i\leq n}f(i)[lpf(i)>j]$ ，答案 $F(n,0)+f(1)$

F (n, j) = G (n, | P |) - G (p_{j}, j - 1) + \sum_{i > j, p_{i}^{e} \leq n} f (p_{i}^{e}) (F (\frac{n}{p_{i}^{e}}, i) + [e \neq 1])

$F(n,j)=G(n,|P|)-G(p_j,j-1)+\sum_{i>j,p_i^e\leq n}f(p_i^e)(F(\frac n{p_i^e},i)+[e\neq 1])$

即是大于 $p_j$ 的质数，然后枚举最小质因子 $p_i>p_j$ ， $[e=1]$ 时，就是素数 $p_i$ 贡献，算过了。

递归求解即可，总复杂度 $O(n^{0,75}/ln\ n+n/poly(ln\ n))$

例题代码

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N = 1e6+5;
const int mod = 1e9+7;

char buf[1<<23],*p1=buf,*p2=buf;
#define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++)
int read(){
    int s=0,w=1; char ch=getchar();
    while(!isdigit(ch)){if(ch=='-')w=-1;ch=getchar();}
    while(isdigit(ch))s=s*10+(ch^48),ch=getchar();
    return s*w;
}

int prime[N],sp1[N],sp2[N];
int w[N],g1[N],g2[N],ind1[N],ind2[N];
int n,sqr,cnt,tot;
bool vis[N];

void getprime(int n){
    for(int i=2;i<=n;++i){
        if(!vis[i])
            prime[++cnt]=i,
            sp1[cnt]=(sp1[cnt-1]+i)%mod,
            sp2[cnt]=(sp2[cnt-1]+i*i%mod)%mod;
        for(int j=1;j<=cnt&&prime[j]*i<=n;++j){
            vis[prime[j]*i]=true;
            if(i%prime[j]==0)break;
        }
    }
}
int calc(int x,int y){
    if(prime[y]>=x)return 0;
    int k=x<=sqr?ind1[x]:ind2[n/x];
    int res=(g2[k]-g1[k]+mod-sp2[y]+sp1[y]+mod)%mod;
    for(int i=y+1;i<=cnt&&prime[i]*prime[i]<=x;++i){
        int mp=prime[i],t=mp%mod;
        for(int j=1;mp<=x;++j,mp=mp*prime[i],t=mp%mod)
            (res+=t*(t-1)%mod*(calc(x/mp,i)+(j!=1))%mod)%=mod;
    }
    return res;
}

signed main(){
    n=read(),sqr=sqrt(n),getprime(sqr);
    for(int i=1;i<=n;){
        int j=n/(n/i);
        w[++tot]=n/i,g1[tot]=w[tot]%mod;
        g2[tot]=(g1[tot]*(g1[tot]+1)%mod*(2*g1[tot]+1)%mod*(166666668)%mod-1+mod)%mod;
        g1[tot]=(g1[tot]*(g1[tot]+1)/2%mod-1+mod)%mod;
        w[tot]<=sqr?ind1[w[tot]]=tot:ind2[n/w[tot]]=tot;
        i=j+1;
    }
    for(int i=1;i<=cnt;++i)
        for(int j=1;j<=tot&&prime[i]*prime[i]<=w[j];++j){
            int k=w[j]/prime[i];
            k=k<=sqr?ind1[k]:ind2[n/k];
            g1[j]-=prime[i]*(g1[k]-sp1[i-1]+mod)%mod;
            g2[j]-=prime[i]*prime[i]%mod*(g2[k]-sp2[i-1]+mod)%mod;
            (g1[j]+=mod)%=mod,(g2[j]+=mod)%=mod;
        }
    printf("%lld\n",(calc(n,0)+1)%mod);
    return 0;
}