min_25筛学习笔记
min25筛是一个神奇的筛法
是较为常用的筛法
假设要求\(f(x)\)
要求
- \(f\)是积性函数
- \(f(p^e)\)是关于\(p\)的低阶多项式
- 存在一个在质数处取值与f一样的完全积性函数
以下以洛谷上的板子题为例
求的积性函数是\(f(p^e)=p^e(p^e-1)=(p^e)^2-p^e\)
主要思路是分情况讨论,按照实数和合数进行分类
问题变成了求
对后面那个再次进行分类讨论,枚举最小质因子及其指数
因为最小质因子是根号级别的所以可以直接枚举
分完情况了我们来进行计算
因为\(f(p)\)是低阶多项式,所以可以拆成单项式,计算\(f(p)=p^k\)的前缀和
这里就是拆成一次和二次
质数部分
因为
首先肯定是不能进行线性筛的
我们来整一个DP
设\(g(n,i)=\sum\limits_{j=1}^n[j\leq n \&\& mindiv(j)>prime_j]j^k\)
设\(sp_i=\sum\limits_{i=1}^nf(p_i)\)
用人话说就是所有小于\(n\)的数中没有小于第\(i\)个质数的质因子的数
后面的\(i^k\)是一个和\(f\)在质数处取值一样的完全积性函数
我们发现\(g(n,x)\)就是最后的答案
其中\(p_x\)是最后一个小于等于\(\sqrt{n}\)的质数
考虑怎样递推
考虑这个递推式的含义
从\(i-1\)到\(i\)少了一些数
因为是完全积性函数,所以可以直接把最小质因子提出来
后面那个是为了防止把质数给弄没了
因为我们需要的\(p\)是根号级别的,所以可以直接线性筛暴力算
但是\(n\)的范围还是太大了,无法都算出来
但我们有结论
所以我们的\(n\)在转移的过程中一直除来除去,也只会用到能表示成\(\lfloor\frac{n}{x}\rfloor\)的值
这显然只有\(\sqrt{n}\)个
我们需要在存储dp数组的时候进行一些奥妙重重的实现
发现第二维每次只会增加1,所以可以直接滚掉
计算答案
同样使用dp来计算答案
设\(S(n,i)\)表示最小质因子大于\(p_i\)的数的函数值之和
那么最终要求得就是\(S(n,0)\)
这里的\(g(n)\)指\(g(n,x)\),\(x\)是小于等于\(\sqrt{n}\)的最大值
后面那个中括号的意思差不多是那个数不是质数
然后就做完了
直接爆搜就可以,复杂度是对的
不能处理多组询问
洛谷板子
#include<cstdio>
#include<cmath>
using namespace std;
namespace Patchouli{
const int N=10000005;
const int MOD=1000000007;
const int INV2=500000004;
const int INV3=333333336;
long long prime[N],mindiv[N],sp1[N],sp2[N],cnt;
long long n;
long long sqr,tot,g1[N],g2[N],w[N],pos1[N],pos2[N];
inline long long read(){
long long a=1,b=0;char t;
do{t=getchar();if(t=='-')a=-1;}while(t>'9'||t<'0');
do{b=b*10-'0'+t;t=getchar();}while(t>='0'&&t<='9');
return a*b;
}
inline void init(int n){
mindiv[1]=1;
for(int i=2;i<=n;++i){
if(!mindiv[i]){
mindiv[i]=i;
prime[++cnt]=i;
sp1[cnt]=(sp1[cnt-1]+i)%MOD;
sp2[cnt]=(sp2[cnt-1]+1ll*i*i)%MOD;
}
for(int j=1;j<=cnt&&i*prime[j]<=n&&prime[j]<=mindiv[i];++j)
mindiv[i*prime[j]]=prime[j];
}
}
long long S(long long x,int y){
if(prime[y]>=x)return 0;
long long k=x<=sqr?pos1[x]:pos2[n/x];
long long ans=(g2[k]-g1[k]+MOD-(sp2[y]-sp1[y])+MOD)%MOD;
for(int i=y+1;i<=cnt&&prime[i]*prime[i]<=x;++i){
long long tmp=prime[i];
for(int e=1;tmp<=x;++e,tmp=tmp*prime[i]){
long long res=tmp%MOD;
ans=(ans+res*(res-1)%MOD*(S(x/tmp,i)+(e!=1)))%MOD;
}
}
return ans%MOD;
}
int QAQ(){
n=read();
sqr=sqrt(n);
init(sqr);
for(long long l=1,r;l<=n;l=r+1){
r=n/(n/l);
w[++tot]=n/l;
long long tmp=(n/l)%MOD;
//求出从2到tmp的自然数之和
g1[tot]=(tmp*(tmp+1)/2)%MOD-1;
//求出从2到tmp的平方和,式子写这么奇怪只是为了少取模
g2[tot]=tmp*(tmp+1)/2%MOD*(2*tmp+1)%MOD*INV3%MOD-1;
//记录一下各个值出现的位置
if(n/l<=sqr)pos1[n/l]=tot;
else pos2[n/(n/l)]=tot;
}
for(int i=1;i<=cnt;++i){
for(int j=1;j<=tot&&prime[i]*prime[i]<=w[j];++j){
long long k=w[j]/prime[i];
k=k<=sqr?pos1[k]:pos2[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;
g2[j]%=MOD;
}
}
printf("%lld\n",(S(n,0)+1+MOD)%MOD);
return false;
}
}
int main(){
return Patchouli::QAQ();
}