数论筛法-min_25筛-学习笔记
i207M
2019-01-12 16:54:27
整除分块是**2**$\sqrt{n}$个数！！！数组要开大
~~今天突然被学习了min_25筛...~~
## 用途
在$O(\frac{n^{\frac 34}}{\log n})$的时间内求**积性函数**$f(x)$的前缀和，要求：
1. 积性函数
2. $f(p)$是一个关于p的简单多项式
3. $f(p^c)$可以快速计算
PS：关于复杂度
min_25筛的复杂度其实是$O(n^{1-\epsilon})$的，只是$10^{12}$内能很快出解
升级版min_25筛可以做到$O(n^{2/3})$~~但是我不会~~
min_25筛重要的一点是可以通过最小质因子遍历每个数。有时可以筛一些非积性函数。
## 内容
[by zhoushuyu](https://www.cnblogs.com/zhoushuyu/p/9187319.html)
首先我们要对每个$x=\lfloor\frac ni\rfloor$求出$\sum_{i=1}^x[\text{i是质数}]f(i)$。
怎么求呢？
先线性筛出$\sqrt{N}$范围内的质数，设$p_j$表示从小到大第$j$个质数。
设$g(n,j)=\sum_{i=1}^{n}[i \in P \ or\ \min(p)>p_j]f(i)$
**说人话就是**：$i$是质数，或者$i$的最小质因子大于$p_j$，把$1-n$内满足条件的$f(i)$加起来就是$g(n,j)$。
**再说人话就是**：埃氏筛（对于质数i，从$i^2$开始筛），用前$j$个质数筛，剩下每个数的$f(i)$之和。
### 首先考虑求1-n的质数个数
*接下来如无特殊说明，都把1开除数字籍*
$g(n,j)=g(n,j-1)-(g(\lfloor\frac n{p_j}\rfloor,j-1)-(j-1))$
考虑埃氏筛的过程，首先前j-1轮的贡献我们已经算出来了，看看第j轮筛掉哪些数，因为我们是从$p^2$开始筛的，所以我们把每个p的倍数都$/p$，这样就递归到了一个子问题。但是我们算重了一些数，事实上，我们所有形如$p_ip_j,i<j$的数又被算了一遍，所以减去它们的数量。
诶？但是**$p_i^2p_j$不用去重吗**？不用。因为除掉$p_j$后，它的最小质因数不大于$p_{j-1}$，它也不是质数，所以没有被算。
这里有一个**复杂度的优化**：当$p_j^2>n$时，后面减去的那一项恒为0，所以式子变为：
$\large g(n,j)=\begin{cases} g(n,j-1)&p_j^2\gt n\\ g(n,j-1)-f(p_j)[g(\frac{n}{p_j},j-1)-\sum_{i=1}^{j-1}f(p_i)]&p_j^2\le n\end{cases}$
**滚动数组**。
同理，根据积性函数的性质，我们也可以筛$g(n,j)=\sum_{i=1}^{n}[i \in P \ or\ \min(p)>p_j]f(i)$。（刚才讨论的式子即为$f(i)=1$）
我们只要找到$p_j^2\le n,p_{j+1}^2>n$的$j$，此时的$g(n,j)$即为1-n的质数的答案。
**关于$g(n,j)$的初值问题**：$g(n,0)$表示所有数的和，也就是把所有数都当作是质数带入$f(p)$的那个多项式中算出的结果。
复杂度被证明是$O(\frac{n^{\frac 34}}{\log n})$
## 关键
现在我们已经对于$x=\lfloor\frac ni\rfloor$求出$\sum_{i=1}^x[\text{i是质数}]f(i)$。
我们设$S(n,j)=\sum_{i=1}^n[\min(p)\ge p_j]f(i)$，也就是所有满足最小质因子大于等于$p_j$的$f$值之和。
那么最终的答案就是$S(n,1)+f(1)$。
鉴于质数的答案我们已经算出来了，是$g(n,j)-\sum_{i=1}^{j-1}f(p_i)$。（因为要保证最小质因子大于等于$p_j$所以要把小于它的质数减掉）
考虑合数。**我们枚举这个合数的最小质因子及其出现次数，然后直接乘即可**。
$\large S(n,j)=g(n,|P|)-\sum_{i=1}^{j-1}f(p_i)+\sum_{k=j}^{p_k^2\le n}\sum_{e=1}^{p_k^{e+1}\le n}S(\frac{n}{p_k^e},k+1)\times f(p_k^e)+f(p_k^{e+1})$
第一重循环是枚举最小的质因数，第二重枚举次数。
注意，最后要加上$f(p_k^{e+1})$，原因是这样才能筛出$p^k$
然后这个的复杂度也被证明是$O(\frac{n^{\frac 34}}{\log n})$的。
## 实战
大家的入门题好像都是LOJ 6053简单的函数
筛$\sum_{i=1}^x[\text{i是质数}]i^k$时，对于每个指数必须单独筛。
~~注意区分n和x~~
```cpp
#define N 200005
int p[N],cntp;
bool notp[N];
int sp[N];
void prework(int n)
{
notp[1]=1;
for(ri i=2; i<=n; ++i)
{
if(!notp[i]) p[++cntp]=i,sp[cntp]=(sp[cntp-1]+i)%md;
for(ri j=1,t; j<=cntp&&(t=p[j]*i)<=n; ++j)
{
notp[t]=1;
if(i%p[j]==0) break;
}
}
}
LL w[N];
int m,g[N],h[N];
int block,id1[N],id2[N];
LL n;
int S(LL x,int y)
{
if(x<=1||p[y]>x) return 0;
int k=(x<=block?id1[x]:id2[n/x]),res=((g[k]-h[k]-(sp[y-1]-y+1))%md+md)%md;
if(y==1) res=(res+2)%md;
for(ri i=y; i<=cntp&&(LL)p[i]*p[i]<=x; ++i)
{
LL t1=p[i],t2=(LL)p[i]*p[i];
for(ri e=1; t2<=x; ++e,t1=t2,t2*=p[i])
res=(res+(LL)(p[i]^e)*S(x/t1,i+1)+(p[i]^(e+1)))%md;
}
return res;
}
signed main()
{
#ifdef M207
freopen("in.in","r",stdin);
// freopen("out.out","w",stdout);
#endif
in(n); block=sqrt(n);
prework(block);
for(LL i=1,j; i<=n; i=j+1)
{
j=n/(n/i);
w[++m]=n/i;
h[m]=(w[m]-1)%md,g[m]=(w[m]%md*((w[m]+1)%md)%md*inv2%md-1)%md;
if(w[m]<=block) id1[w[m]]=m;
else id2[j]=m;
}
for(ri j=1; j<=cntp; ++j)
for(ri i=1; i<=m&&(LL)p[j]*p[j]<=w[i]; ++i)
{
int k=(w[i]/p[j]<=block?id1[w[i]/p[j]]:id2[n/(w[i]/p[j])]);
g[i]=(g[i]-(LL)p[j]*(g[k]-sp[j-1])%md+md)%md;
h[i]=(h[i]-(h[k]-(j-1))%md+md)%md;
}
int ans=S(n,1)+1;
out(ans);
return 0;
}
```
### SPOJ DIVCNTK
```cpp
#define N 200005
int p[N],cntp;
bool notp[N];
void prework(int n)
{
notp[1]=1;
for(ri i=2; i<=n; ++i)
{
if(!notp[i]) p[++cntp]=i;
for(ri j=1,t; j<=cntp&&(t=p[j]*i)<=n; ++j)
{
notp[t]=1;
if(i%p[j]==0) break;
}
}
}
LL w[N];
int m,block=100000;
int id1[N],id2[N];
LL n;
ull K,f[N];
ull S(LL x,int y)
{
if(x<=1||x<p[y]) return 0;
int k=(x<=block)?id1[x]:id2[n/x]; ull res=(K+1)*(f[k]-y+1);
for(ri i=y; i<=cntp&&(LL)p[i]*p[i]<=x; ++i)
{
LL t1=p[i],t2=(LL)p[i]*p[i];
for(ri e=1; t2<=x; ++e,t1=t2,t2*=p[i])
res+=(e*K+1)*S(x/t1,i+1)+((e+1)*K+1);
}
return res;
}
void solve()
{
scanf("%lld%llu",&n,&K);
m=0;
for(LL i=1,j; i<=n; i=j+1)
{
j=n/(n/i); w[++m]=n/i;
f[m]=w[m]-1;
if(w[m]<=block) id1[w[m]]=m;
else id2[j]=m;
}
for(ri j=1; j<=cntp; ++j)
for(ri i=1; i<=m&&(LL)p[j]*p[j]<=w[i]; ++i)
{
int k=(w[i]/p[j]<=block)?id1[w[i]/p[j]]:id2[n/(w[i]/p[j])];
f[i]=f[i]-f[k]+j-1;
}
printf("%llu\n",S(n,1)+1);
}
```
杜教筛的模板题，跑得飞快：
```cpp
LL Phi(int x,int y)
{
if(x<=1||x<p[y]) return 0;
int k=(x<N)?id1[x]:id2[n/x]; LL res=g[k]-h[k]-(fp[y-1]-(y-1));
for(ri i=y; (LL)p[i]*p[i]<=x; ++i)
{
LL t1=p[i],t2=(LL)p[i]*p[i];
for(ri e=1; t2<=x; t1=t2,t2*=p[i],++e)
res+=t1/p[i]*(p[i]-1)*Phi(x/t1,i+1)+t1*(p[i]-1);
}
return res;
}
LL Mu(int x,int y)
{
if(x<=1||x<p[y]) return 0;
int k=(x<N)?id1[x]:id2[n/x]; LL res=-(h[k]-(y-1));
for(ri i=y; (LL)p[i]*p[i]<=x; ++i)
res-=Mu(x/p[i],i+1);
return res;
}
void solve()
{
m=0;
for(LL i=1,j; i<=n; i=j+1)
{
j=n/(n/i); w[++m]=n/i;
g[m]=w[m]*((LL)w[m]+1)/2-1,h[m]=w[m]-1;
if(w[m]<N) id1[w[m]]=m;
else id2[j]=m;
}
for(ri j=1; j<=cntp; ++j)
for(ri i=1; (LL)p[j]*p[j]<=w[i]; ++i)
{
int k=(w[i]/p[j]<N)?id1[w[i]/p[j]]:id2[n/(w[i]/p[j])];
g[i]-=p[j]*(g[k]-fp[j-1]), h[i]-=h[k]-(j-1);
}
out(Phi(n,1)+1,Mu(n,1)+1);
}
```