Min25 筛学习笔记
用途
(Min25) 筛 可以在 (Oleft(frac{n^{frac{3}{4}}}{log n} ight)) 或者 (Thetaleft(n^{1 - epsilon} ight)) 的时间复杂度下解决一类积性函数 的前缀和问题。
要求 (f(p)) 是一个关于 (p) 的项数较少的多项式或可以快速求值；(f(p^c)) 可以快速求值。
实现
如果答案是一个多项式，那么我们要把这个多项式的每一项都拿出来求一次答案，最后把答案累加。
设 (s(n,k)) 表示 (1sim 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] imes 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
",addmod(1,solve(n,0)));
return 0;
}