Min_25筛
前言
最近神仙Min_25筛的题越来越多了,555~只好跑过来学一下~(´ー`~)
听说时间复杂度是$O(\frac{n^\frac{3}{4}}{log n})$,然而并不会证。
Min_25筛可以用来计算积性函数$f(x)$的前缀和,要求可快速计算$f(p^c)$。
做法
首先我们要对每个$x=\left\lfloor\frac{n}{i}\right\rfloor$快速计算$\sum\limits_{i=1}^{x}[i∈P]f(i)$。
我们先线性筛出$[1,\sqrt{n}]$内的所有质数,设第j小的质数为$P_j$。
设$g(n,j)=\sum\limits_{i=1}^{n}[i∈P||mn(i)>P_j]f(i)$。$mn(i)$表示i最小的质因子。
$g(n,j)$的意义为,运行了j次线性筛后,没被筛掉的数的f的和。
容易发现,$\sum\limits_{i=1}^{x}[i∈P]f(i)=g(x,|P|)$。$|P|$为小于等于$\sqrt{x}$的质数个数。
我们考虑$g(n,j)$的转移过程:
当$P_j^2>n$时,不会筛掉任何数,因此$g(n,j)=g(n,j-1)$。
当$P_j^2<=n$时,被筛掉的数一定含有质因子$P_j$,且最小质因子$≥P_j$。考虑减去$f(P_j)g(\frac{n}{P_j},j-1)$,由$g(n,j)$的表达式可发现多减去了$\sum\limits_{i=1}^{j-1}f(P_i)$,即除以$P_j$后为小于$P_j$的质数的数。维护质数的f的前缀和,每次转移时加上多减的部分即可。
$g(n,j)=\begin{cases} g(n,j-1)&P_j^2>n\\ g(n,j-1)-f(P_j)[g(\frac{n}{P_j},j-1)-\sum\limits_{i=1}^{j-1}f(P_i)]&P_j^2<=n\end{cases}$
初始情况下,$g(n,0)$为$[1,n]$所有数的f值。
求质数的个数
以求质数个数为例,即$f(x)=1$的情况:
1 const int mod=998244353; 2 namespace min_25{ 3 int cnt=0,tot=0,sqr,pr[1000010],g[1000010],id1[1000010],id2[1000010]; 4 bool vis[1000010]={false};ll w[1000010]; 5 void get_prime(int n){ 6 for(int i=2;i<=n;i++){ 7 if(!vis[i]) 8 pr[++cnt]=i; 9 for(int j=1;j<=cnt&&i*pr[j]<=n;j++){ 10 vis[i*pr[j]]=true; 11 if(i%pr[j]==0) 12 break; 13 } 14 } 15 return; 16 } 17 int solve(ll n){ 18 get_prime(sqr=(int)sqrt(n)); 19 for(ll i=1,j;i<=n;i=j+1){ 20 j=n/(n/i);w[++tot]=n/i; 21 if(w[tot]<=sqr) 22 id1[w[tot]]=tot; 23 else 24 id2[n/w[tot]]=tot; 25 g[tot]=(w[tot]-1)%mod; 26 } 27 for(int j=1;j<=cnt;j++) 28 for(int i=1;i<=tot&&(ll)pr[j]*pr[j]<=w[i];i++){ 29 int id=(w[i]/pr[j]<=sqr)?id1[w[i]/pr[j]]:id2[n/(w[i]/pr[j])]; 30 g[i]=(g[i]-g[id]+j-1+mod)%mod; 31 } 32 return g[1]; 33 } 34 }//f(x)=1
代码注释:
$(1)$$w$存储所有$\left\lfloor\frac{n}{i}\right\rfloor$的值。由于$w$的值过大,分为$<=\sqrt{n}$(即$id1$)和$>\sqrt{n}$(即$id2$)两部分存储下标。
$(2)$我们最后需要的只有$g(n,|P|)$,且$g(x,j)$的值只与$g(y,j-1)$有关,因此开一维数组就够了。
$(3)$$f(1)$另算。
求$f(x)$的前缀和
上面我们求出了范围内所有质数的f的和,那么怎么求f的前缀和呢?
设$S(n,j)=\sum\limits_{i=1}^{n}[mn(i)>=P_j]f(i)$。显然所求为$S(n,1)+f(1)$。
我们将$S(n,j)$分为质数和合数两部分处理。
质数部分:即大于等于$P_j$的质数的f的和,为$g(n,|P|)-\sum\limits_{i=1}^{j-1}f(P_i)$。
合数部分:枚举最小质因子及其次数,根据积性函数的性质直接递归计算。
$\sum\limits_{i=j}^{P_i^2<=n}\sum\limits_{k=1}^{P_i^{k+1}<=n}S(\frac{n}{P_i^k},i+1)+f(P_i^{k+1})$
于是我们得到了$S(n,j)$的表达式:$S(n,j)=g(n,|P|)-\sum\limits_{i=1}^{j-1}f(P_i)+\sum\limits_{i=j}^{P_i^2<=n}\sum\limits_{k=1}^{P_i^{k+1}<=n}S(\frac{n}{P_i^k},i+1)+f(P_i^{k+1})$
时间复杂度依然是$O(\frac{n^\frac{3}{4}}{log n})$。(别问我,我还是不会证φ(..) )
接下来看道Min_25模板题(被自己的各种智障操作弄到怀疑人生)
loj6053 简单的函数
有一个积性函数$f(x)$,满足$f(1)=1$,$f(p^c)=p⊕c$(⊕表示异或,p为质数,$c>0$),求$f(x)$的前n项和。
我们发现,大于2的质数p均满足$f(p)=p-1$,而$f(2)=3=1+2$。我们可令$f(2)=1$,统计答案时再加上2。
我们需要求出g数组,但$f'(x)=x-1$并不满足积性函数的性质,怎么办呢?
我们可以分别计算出$g(x,|P|)=\sum\limits_{i=1}^{x}[i∈P]i$和$h(x,|P|)=\sum\limits_{i=1}^{x}[i∈P]$,计算时相减即可。
1 #include<cstdio> 2 #include<cmath> 3 #include<cstring> 4 #include<algorithm> 5 using namespace std; 6 typedef long long ll; 7 inline ll rd(){ 8 char c=getchar();ll x=0,flag=1; 9 for(;c<'0'||c>'9';c=getchar())if(c=='-')flag=-1; 10 for(;c>='0'&&c<='9';c=getchar())x=x*10+c-'0'; 11 return x*flag; 12 } 13 const int mod=1000000007; 14 const int N=1000000; 15 namespace min_25{ 16 int sqr,cnt=0,tot=0,pr[N+10],id1[N+10],id2[N+10],pre[N+10]={0},f[N+10],g[N+10]; 17 bool vis[N+10]={false};ll w[N+10]; 18 void get_prime(int n){ 19 for(int i=2;i<=n;i++){ 20 if(!vis[i]){ 21 pr[++cnt]=i; 22 pre[cnt]=(pre[cnt-1]+i)%mod; 23 } 24 for(int j=1;j<=cnt&&i*pr[j]<=n;j++){ 25 vis[i*pr[j]]=true; 26 if(i%pr[j]==0) 27 break; 28 } 29 } 30 return; 31 } 32 int S(ll n,ll x,int j){ 33 if(x<pr[j]) 34 return 0; 35 int id=(x<=sqr)?id1[x]:id2[n/x]; 36 int res=(f[id]-pre[j-1]-g[id]+j-1)%mod; 37 (res+=mod)%=mod; 38 for(int i=j;i<=cnt&&(ll)pr[i]*pr[i]<=x;i++){ 39 ll val=pr[i]; 40 for(int k=1;val*pr[i]<=x;k++,val*=pr[i]) 41 res=(res+((ll)S(n,x/val,i+1)*(pr[i]^k)%mod+(pr[i]^(k+1)))%mod)%mod; 42 } 43 return res; 44 } 45 int solve(ll n){ 46 if(n==1)return 1; 47 get_prime(sqr=(int)sqrt(n)); 48 for(ll i=1,j;i<=n;i=j+1){ 49 j=n/(n/i);w[++tot]=n/i; 50 if(w[tot]<=sqr) 51 id1[w[tot]]=tot; 52 else 53 id2[n/w[tot]]=tot; 54 f[tot]=((w[tot]+2)%mod)*((w[tot]-1)%mod)%mod; 55 if(f[tot]&1)f[tot]+=mod;f[tot]>>=1; 56 g[tot]=(w[tot]-1)%mod; 57 } 58 for(int j=1;j<=cnt;j++){ 59 for(int i=1;i<=tot&&(ll)pr[j]*pr[j]<=w[i];i++){ 60 int id=(w[i]/pr[j]<=sqr)?id1[w[i]/pr[j]]:id2[n/(w[i]/pr[j])]; 61 f[i]=(f[i]-(ll)pr[j]*(f[id]-pre[j-1]+mod)%mod+mod)%mod; 62 g[i]=(g[i]-g[id]+j-1)%mod;(g[i]+=mod)%=mod; 63 } 64 } 65 return (S(n,n,1)+3)%mod; 66 } 67 } 68 int main(){ 69 printf("%d\n",min_25::solve(rd())); 70 return 0; 71 }
待更……(我才不会说我咕了呢)(其实是自闭了调不出来啦555~)(毒瘤搬题人出题人)