Min25 筛学习笔记
Min25 筛学习笔记
用途
\(Min25\) 筛 可以在 \(O\left(\frac{n^{\frac{3}{4}}}{\log n}\right)\) 或者 \(\Theta\left(n^{1 - \epsilon}\right)\) 的时间复杂度下解决一类积性函数 的前缀和问题。
要求 \(f(p)\) 是一个关于 \(p\) 的项数较少的多项式或可以快速求值；\(f(p^c)\) 可以快速求值。
实现
如果答案是一个多项式，那么我们要把这个多项式的每一项都拿出来求一次答案，最后把答案累加。
设 \(s(n,k)\) 表示 \(1\sim n\) 中所有的质数和最小质因子大于 \(pri[k]\) 的所有合数的答案。
最终我们要求的就是 \(s(n,0)\)。
枚举最小质因子，\(s(n,k)=f(i)(i\ is\ a\ prime \& i>pri[k+1])+\sum_{i=k+1}^{pri[i] \times pri[i] \leq n}\sum_{j=1}^{pri[i]^j \leq n} f(pri[i]^j) (s(\frac{n}{pri[i]^j},i)+j!=1)\)
之所以要加上一个 \(j!=1\)是因为类似于 \(pri[i]^j,j \geq 2\) 的贡献并不会被计算。
对于后面的这个部分，我们直接递归求解即可，不用记忆化复杂度也是正确的，而且 \(f(pri[i]^j)\) 根据我们的要求是可以快速计算的。
问题就在于前面的部分怎么求，可以用一下简单的容斥。
用 \(1 \sim n\) 中所有质数的答案减去小于等于 \(pri[k]\) 的答案。
对于小于等于 \(pri[k]\) 的答案，因为合法的质数大小是小于等于 \(\sqrt{n}\) 的，所以可以用线性筛预处理。
\(1 \sim n\) 中所有质数的答案可以用 \(dp\) 预处理。
设 \(g(n,k)\) 为 \(1 \sim n\) 中所有质数和最小质因子大于 \(pri[k]\) 的合数的答案。
因为只要求出所有质数的答案，所以我们可以把这个函数看成一个完全积性函数。
\(g(n,0)\) 的答案是一个自然数幂和的形式，指数较小的时候可以直接用一个前缀和的公式 \(O(1)\) 计算，指数较大的时候就要考虑用伯努利数。
转移方程
\(g(n,k)=g(n,k-1)-f(pri[k])(g(\frac{n}{pri[k]},k-1)-g(pri[k-1],k-1))\)。
含义就是去掉最小质因子恰好为 \(pri[k]\) 的合数的贡献。
因为小于等于 \(pri[k-1]\) 的质数和 \(pri[k]\) 相乘得到的合数最小质因子小于 \(pri[k]\) ,所以它的贡献不能被计算，所以我们要把它减去，减去的部分恰好就是我们用线性筛处理出的前缀和。
同时因为是完全积性函数，所以我们不用去考虑后面的地方有没有把这个东西除尽。
发现可以把第二维滚掉，同时有用的第一维只有 \(\frac{n}{k}\) 的形式，也就是说我们可以把第一维的下标映射到两个数组中。
设 \(maxcnt\) 为 \(\leq \sqrt{n}\) 的质数的个数，我们要的只是 \(g(n,maxcnt)\) 的值，所以第一维和第二维总共的状态不会很多。
直接递推即可，注意边界。
最后的时候还要把 \(f(1)\) 的值加上，因为我们并没有统计到它。
代码
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<iostream>
#define rg register
const int maxn=1e6+5,mod=1e9+7,inv2=500000004,inv6=166666668;
inline int delmod(rg int now1,rg int now2){
return now1-=now2,now1<0?now1+mod:now1;
}
inline int addmod(rg int now1,rg int now2){
return now1+=now2,now1>=mod?now1-mod:now1;
}
inline int mulmod(rg long long now1,rg int now2){
return now1*=now2,now1>=mod?now1%mod:now1;
}
typedef long long ll;
bool not_pri[maxn];
int pri[maxn],sp1[maxn],sp2[maxn],g1[maxn],g2[maxn],tot,v1[maxn],v2[maxn],sqr;
ll w[maxn];
void pre(rg int mmax){
not_pri[0]=not_pri[1]=1;
for(rg int i=2;i<=mmax;i++){
if(!not_pri[i]){
pri[++pri[0]]=i;
sp1[pri[0]]=addmod(sp1[pri[0]-1],i);
sp2[pri[0]]=addmod(sp2[pri[0]-1],mulmod(i,i));
}
for(rg int j=1;j<=pri[0] && 1LL*i*pri[j]<=mmax;j++){
not_pri[i*pri[j]]=1;
if(i%pri[j]==0) break;
}
}
}
ll n;
int solve(rg ll x,rg int y){
if(y && x<=pri[y]) return 0;
rg int now=x<=sqr?v1[x]:v2[n/x];
rg int ans=delmod(delmod(g2[now],g1[now]),delmod(sp2[y],sp1[y]));
for(rg int i=y+1;i<=pri[0] && 1LL*pri[i]*pri[i]<=x;i++){
rg ll np=pri[i];
for(rg int j=1;np<=x;j++,np*=pri[i]){
rg ll tmp=np%mod;
ans=addmod(ans,tmp*(tmp-1)%mod*addmod(solve(x/np,i),j!=1)%mod);
}
}
return ans;
}
int main(){
scanf("%lld",&n);
sqr=sqrt(n);
pre(sqr+100);
rg ll now;
for(rg ll l=1,r;l<=n;l=r+1){
r=(n/(n/l));
now=n/l;
w[++tot]=now;
if(now<=sqr) v1[now]=tot;
else v2[n/now]=tot;
now%=mod;
g1[tot]=delmod((now+1)*now%mod*inv2%mod,1);
g2[tot]=delmod((2*now+1)%mod*(now+1)%mod*now%mod*inv6%mod,1);
}
for(rg int i=1;1LL*pri[i]*pri[i]<=n;i++){
for(rg int j=1;j<=tot && 1LL*pri[i]*pri[i]<=w[j];j++){
rg ll now=w[j]/pri[i];
now=now<=sqr?v1[now]:v2[n/now];
g1[j]=delmod(g1[j],mulmod(pri[i],delmod(g1[now],sp1[i-1])));
g2[j]=delmod(g2[j],mulmod(1LL*pri[i]*pri[i]%mod,delmod(g2[now],sp2[i-1])));
}
}
printf("%d\n",addmod(1,solve(n,0)));
return 0;
}
