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~)(毒瘤搬题人出题人)

posted @ 2019-04-03 14:20  乖巧的小团子QwQ  阅读(261)  评论(0编辑  收藏  举报