min_25筛
min_25筛
前言：听说是一个叫 min_25 的日本人发明的，今天来口胡一下。
此算法是用来求一个特殊的函数的前缀和的，之所以说特殊，是因为要满足一些特殊条件。设函数为 \(f(x)\)，那我们要求出 \(\sum\limits_{i=1}^n f(x)\) 。
- 要满足以下条件
- \(f(p^k)\) （\(p\)为质数，k为常数）要求能在 O( 常数 ) 的复杂度内算出。
- 网上大多说 \(f(x)\) 要满足是积性函数，但我认为不太准确，首先若 \(f(x)\) 满足是积性函数，那一定是可以求的，若 \(f(x)\) 为常数个积性函数的和，其实也是可以分别求出再求和。
大致的思路是分别求先求质数的，再求合数的。
怎么求质数的呢，我认为这个思路应该是从埃氏筛法筛质数的思路上得到启发。
假设 \(f(p^k)=\sum\limits_{i=0}^m ai*p^i\) ，（\(p\)为质数，m很小）我们考虑分别求出质数们的某次方的和，到时候需要的时候就乘上相对应的系数\(ai\)再加起来即可，
那么问题就转化为求质数的m次方的前缀和，设 \(g(x)=x^m\) ，我们要求出 \(\sum\limits_{p为质数}g(p)\)，这时发现 \(g(x)\) 是一个完全积性函数，有很好的性质。
设 \(G(n,i)\) 表示，从所有数的和中，筛去了前 \(i\) 个质数的倍数后的结果（是不是很有埃筛的感觉），设 \(p[i]\) 表示第 \(i\) 小的质数。
初值是 \(G(n,0)=\sum\limits_{i=1}^ng(i)\) 你发现这可以 O( 1 ) 计算，而我们要的是 \(G(n,tot)\)。考虑怎么转移，其实很巧妙：
第一行显然，在埃筛里只需用<=\(\sqrt{n}\) 的质数去筛数。
第二行是在考虑筛掉 \(p[i]\) 的倍数，因为是完全积性函数，所以可以这么算，之所以减掉后面那玩意是因为之前已经用更小的质数筛过了。
我们再观察一下式子，发现 \(p[i]<=\sqrt{n}\) 即可，而 \(G(n,i)\) 的 \(n\) 只需保留 \(\lfloor \frac{n}{x} \rfloor\) 就行了，时间和状态都不多。
至此，质数的答案就求完了。
我们再考虑如何求合数。
有过之前求质数的经验，我们自然想到沿用思路，但要注意 \(f(x)\) 不是完全积性函数了。
设 \(F(n,i)\) 表示从所有数的和中，筛去了前 \(i-1\) 个质数的倍数后的结果，（注意这里与上面有所不同，主要是为了方便实现），可以知道，\(F(n,1)\) 即为所求。
先给出转移：
让我们慢慢解释一下式子，这个过程是将最小质因子为 \(p[j]\) 的数加进来。
- 首先前面大括号中的是之前已经求过的质数的答案。
- 其后是先枚举最小质因子\(p[j]^c\)，发现会漏算最后那个小玩意 \(f(p[j]^{c+1})\)，所以要加上，之所以这么麻烦，是因为 \(f(x)\) 不是完全积性函数。
- 这里你可以体会 \(f(x)\) 的那个条件，\(f(p[j]^c)*F(\frac{n}{p[j]^c},j+1)\) 是对的就行。
这样合数就求完了。
然而有一点小细节，我们发现这么求，\(f(1)\) 并没有被计算到，因为没有数的最小质因子是 1 ，
所以计算时要减去 \(f(1)\) ，你可以在算 \(G(n,i)\) 时提前减掉，
最后记得往答案里加上 \(f(1)\) 。
小结一下，感受一下整个过程，求质数的过程是从全部都有，慢慢地筛到只剩质数；而求合数的过程，则是从啥都没有，一点一点加到全部都有，
这中感觉很奇妙；设状态是沿用埃筛思想，在转移时巧用积性函数的性质进行DP转移。
来一道板题：P5325 【模板】Min_25筛
题目描述
定义积性函数\(f(x)\)，满足\(f(p^k)=p^k*(p^k-1)\)（\(p\)为质数），要求出\(\sum\limits_{i=1}^{n}f(i)\)，对\(10^9+7\)取模。
小解析
\(f(p)=p^2-p^1\) ，因此求质数时只需求出，\(\sum\limits_{p为质数}g(p^1)\) 和 \(\sum\limits_{p为质数}g(p^2)\)
//f(p^k)=p^k(p^k-1)
#include<cstdio>
#include<cmath>
#define ll long long
#define pl(x) ((x>=mo)? x-mo:x)
using namespace std;
const int N=2e5+10,mo=1e9+7;
int bz[N],su[N],id1[N],id2[N],a[N][3],g[N][3]; ll w[N];
int k,p,sqr,P,m,x,inv2,inv6,an; ll n,i,j;
inline void pre(){
int i,j;
for(i=2;i<=sqr;++i){
if (!bz[i]){
su[++P]=i;
a[P][0]=a[P-1][0]+1;
a[P][1]=pl(a[P-1][1]+i);
a[P][2]=pl(a[P-1][2]+1ll*i*i%mo);
}
for(j=1;j<=P&&1ll*su[j]*i<=sqr;++j){
bz[i*su[j]]=1;
if (i%su[j]==0) break;
}
}
}
inline int ksm(int x,int y){
int s=1;
while(y){
if (y&1) s=1ll*s*x%mo;
x=1ll*x*x%mo;
y>>=1;
}
return s;
}
int S(ll x,int t){
if (x<su[t]||x<=1) return 0;
int ss=0,k,i,s,p; ll pp;
if (x<=sqr) k=id1[x];
else k=id2[n/x];
ss=(1ll*(g[k][2]-a[t-1][2])-1ll*(g[k][1]-a[t-1][1]))%mo; ss=pl(ss+mo);
for(i=t;i<=P&&1ll*su[i]*su[i]<=x;++i)
for(pp=su[i],p=1ll*su[i]*su[i]%mo;pp*su[i]<=x;pp=pp*su[i],p=1ll*p*su[i]%mo){
s=pp%mo*(pp%mo-1)%mo; s=pl(s+mo);
ss=(ss+1ll*s*S(x/pp,i+1))%mo;
ss=(ss+1ll*p*(p-1))%mo; ss=pl(ss+mo);
}
return ss;
}
int main(){
scanf("%lld",&n); sqr=sqrt(n); pre();
inv2=ksm(2,mo-2); inv6=ksm(6,mo-2);
for(i=1;i<=n;i=j+1){
j=n/(n/i); w[++m]=n/i;
if (w[m]<=sqr) id1[w[m]]=m;
else id2[n/w[m]]=m;
x=w[m]%mo;
g[m][0]=x-1;
g[m][1]=1ll*(x+2)*(x-1)%mo*inv2%mo;
g[m][2]=1ll*x*(x+1)%mo*(x*2+1)%mo*inv6%mo;
g[m][2]=pl(g[m][2]-1+mo);
}
for(i=1;i<=P;++i){
p=su[i];
for(j=1;j<=m&&1ll*p*p<=w[j];++j){
if (w[j]/p<=sqr) k=id1[w[j]/p];
else k=id2[n/(w[j]/p)];
g[j][0]=(g[j][0]-1ll*(g[k][0]-a[i-1][0]))%mo; g[j][0]=pl(g[j][0]+mo);
g[j][1]=(g[j][1]-1ll*p*(g[k][1]-a[i-1][1]))%mo; g[j][1]=pl(g[j][1]+mo);
g[j][2]=(g[j][2]-1ll*p*p%mo*(g[k][2]-a[i-1][2]))%mo; g[j][2]=pl(g[j][2]+mo);
}
}
an=pl(S(n,1)+1);
printf("%d",an);
return 0;
}
浙公网安备 33010602011771号