min25筛

发表于 更新于 阅读次数：

min25筛是一种针对特殊的积性函数的筛法，可以在O(\frac{n^{\frac{3}{4}}}{logn})的复杂度内筛出某函数的前缀和，在OI中应用不及杜教筛，但是仍然有学习价值。

min25筛的要求：

1.对于f(p^k)(p\in prim)有明确定义，并且比较好求。

2.该函数在f(p)(p \in prim)时，为完全积性函数，或者可以拆成几个完全积性函数。

由此可见，上述要求较为苛刻，这种函数除非特意构造，不然很难自然的出现在题目之中。

但是min25筛相比杜教筛的优势是，不用去特意构造狄利克雷卷积，也就不依赖数论函数(\mu,\phi等)。

实现:

Part1:处理质数部分对f的贡献

因为f(p)为完全积性的，我们构造一个g，使得g在任何地方的定义都是与f在质数时的定义相同，我们求出g的前缀和，然后依次筛掉g合数部分的值，剩下的部分和\sum_{p \in prim} f(p)是等价的。

我们利用类似于埃氏筛的思想，依次筛掉最小质因子为p_i的贡献，因为最小质因子一定<=\sqrt{n}，我们只要筛出<=\sqrt{n}的质数就可以了。

此处比较神仙：

我们设

G(n,j)=\sum_{i=1}^{n}g(i)[i \in prim ||minp(i)>p_j]

考虑把G(n,j-1)变成G(n,j)，其实就是要筛掉p_j作为最小质因子的数的贡献，而g是完全积性的，我们可以直接提出一个p_j，枚举其他因子：

p_j*G(\frac{n}{p_j},j-1)

但是这样会算多，因为\frac{n}{p_j}>=p_j，而G定义中会包含i\in prim的g_i,所以p_{1}到p_j-1也会被算进去，在这部分中，p_j就不是最小质因子了，和定义不符，需要减去。

所以我们要预先处理出<=\sqrt {n}的\sum_{p \in prim} f(p)，记为sp(i)，表示前i小的质数的f值之和。

于是转移就可以写成

G(n,j)=G(n,j-1)-p_j*[G(\frac{n}{p_j},j-1)-sp(j-1)]

可以发现这个东西可以滚动数组，我们只要预先求出所有形如G(\frac{n}{i},0)的函数值(就是g前\frac{n}{i}项的前缀和)，然后枚举n从大到小，直接更新即可。

Part2 依次加入合数对f的贡献

我们记G(n,|P|)表示筛完了合数后的g(i)的和，显然这东西就等于质数部分的f(i)之和。

然后这里再次神仙一波：

我们设

F(n,j)=\sum_{i=1}^{n}f(i)[minp(i)>=p_j]

注意：此处的定义于上面不同，质数将不再无脑合法，且>变成了>=。

然后我们枚举最小质因子及其次数，有转移:

F(n,j)=G(n,|P|)-sp(j-1)+\sum_{k=j}^{cnt}\sum_{e}^{p_k^{e+1}<=n}[f(p_k^e)*F(\frac{n}{p_k^e} , j+1)+f(p_k^{e+1})]

看起来蛮恐怖的，其实联系实际意义想还好，p_k^1的贡献已经在前面了，只是把单含p_k这个质因子的数和含其他因子的数分开算了而已，当p_k^{e+1}>n的时候p_k^e*p_{k+1}显然也>n，就没有必要枚举了。

注意：不管是筛g还是f，1的贡献最好单独计算。

代码实现

写一下处理所有形如\frac{n}{i}的g的前缀和这部分。

这里因为开不下数组，我们离散化一下，对于\frac{n}{i}<=\sqrt{n}的直接记录下标，对于\frac{n}{i}>=\sqrt{n}，我们记录这个i即可。

注意：我们处理这个时使用数论分块，对于\frac{n}{i}>=\sqrt{n}我们记录的i为分块时的r，因为我们引用下标的时候用\frac{n}{\frac{n}{i}}就可以得到这个r。

另外一个需要记得的是，枚举最小质因子及其次数，要到\sqrt{n}就停止。

题目

[Loj6053]简单的函数

题目链接

solution:此处注意到f(p) =p-1并不是完全积性的，我们令g(i)=i,h(i)=1，并令sp(i)表示前i个素数之和。

我们按照上面的方法设出G(n,j)和H(n,j)等到筛完了G和H后，我们发现他们的定义是:

G(n,|P|)=\text{n以内的质数和}

H(n,|P|)=\text{n以内的质数个数}

有了这两个我们就可以筛f了：

F(n,j)=G(n,|P|)-H(n,|P|)+sp(j-1)-(j-1)+\sum_{k=j}^{cnt}\sum_{e}^{p_k^{e+1}<=n}[f(p_k^e)*F(\frac{n}{p_k^e} , j+1)+f(p_k^{e+1})]

code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
#define gc getchar()
inline int read(){
    register int x(0),f(1);register char c(gc);
    while(c>'9'||c<'0')f=c=='-'?-1:1,c=gc;
    while(c>='0'&&c<='9')x=x*10+c-48,c=gc;
    return f*x;
}
#define ll long long
const int N=210000,P=1e9+7,inv2=(P+1)/2;
ll n,block,w[N];
ll prim[N],id1[N],id2[N],h[N],g[N],m;
int vis[N],cnt,sp[N];
void pre(int n){
    register int i,j;
    for(i=2;i<=n;i++){
        if(!vis[i])prim[++cnt]=i,sp[cnt]=(sp[cnt-1]+i)%P;
        for(j=1;j<=cnt&&i*prim[j]<=n;j++){
            vis[i*prim[j]]=1;
            if(i%prim[j]==0)break;
        }
    } 
}
int S(ll x,ll y){
    if(x<=1||prim[y]>x)return 0;
    int k=(x<=block)?id1[x]:id2[n/x],ret=(g[k]-h[k]+P-sp[y-1]+P+y-1)%P;
    if(y==1)ret+=2;
    for(int i=y;i<=cnt&&1ll*prim[i]*prim[i]<=x;i++){
        ll t1=prim[i],t2=1ll*prim[i]*prim[i];
        for(int e=1;t2<=x;e++,t1=t2,t2*=prim[i])
            (ret+=((1ll*S(x/t1,i+1))*(prim[i]^e)%P+(prim[i]^(e+1))%P))%=P;
    }
    return ret;
}
int main(){
    cin>>n;
    block=sqrt(n);
    pre(block);
    ll l,r;
    for(l=1;l<=n;l=r+1){
        r=n/(n/l);
        w[++m]=n/l;
        h[m]=(w[m]-1)%P;
        g[m]=w[m]%P*(w[m]+1)%P*inv2%P-1;
        if(w[m]<=block)id1[w[m]]=m;
        else id2[r]=m;
    }
    register int i,j;
    for(j=1;j<=cnt;j++)
        for(i=1;i<=m&&1ll*prim[j]*prim[j]<=w[i];i++){
            int k=(w[i]/prim[j]<=block)?id1[w[i]/prim[j]]:id2[n/(w[i]/prim[j])];
          	(g[i]-=1ll*prim[j]*(g[k]-sp[j-1])%P)%=P;
            h[i]=(h[i]-(h[k]-(j-1))+P)%P;
        } 
    ll ans=S(n,1)+1;
    cout<<ans;
    return 0;
}
[luogu5325]min25筛模板

题目链接

题意:求\sum_{i=1}^{n}f(i),f(p^k)=p^k*(p^k-1),f为积性函数。

solution:考虑把f(p^k)拆成g(p^k)=p^{2k}和h(p^k)=p^k,其实就是id和id^2。

这样当k=1的时候他们就是完全积性的，和上面类似的方法筛就可以了。

但是要注意在处理S的时候p_k^e*p_k^e可能会爆longlong，得记得取模。

code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
#include<cstdio>
#include<iostream>
#include<cmath>
#include<cstring>
#include<algorithm>
using namespace std;
#define gc getchar()
#define ll long long
inline ll read(){
    register ll x(0),f(1);register char c(gc);
    while(c>'9'||c<'0')f=c=='-'?-1:1,c=gc;
    while(c>='0'&&c<='9')x=x*10+c-48,c=gc;
    return f*x;
}
const int N=410000,P=1e9+7,inv2=(P+1)/2,inv6=166666668;
ll w[N],m,h[N],g[N],id1[N],id2[N],n,sp[N],prim[N],sp2[N];
int cnt,block,vis[N];
void pre(int n){
    register int i,j;
    for(i=2;i<=n;i++){
        if(!vis[i])prim[++cnt]=i,sp[cnt]=(sp[cnt-1]+1ll*i*i%P)%P,sp2[cnt]=(sp2[cnt-1]+1ll*i)%P;
        for(j=1;j<=cnt&&1ll*i*prim[j]<=n;j++){
            vis[i*prim[j]]=1;
            if(i%prim[j]==0)break;
        }
    }
}
ll S(ll x,ll y){
    if(x<=1||prim[y]>x)return 0;
    int k=x<=block?id1[x]:id2[n/x];
    ll ans=(g[k]-sp[y-1]+P-h[k]+P+sp2[y-1])%P;
    for(int i=y;i<=cnt&&prim[i]*prim[i]<=x;i++)
        for(ll t1=prim[i],t2=prim[i]*prim[i];t2<=x;t1*=prim[i],t2*=prim[i]){
            (ans+=t1%P*(t1%P-1)%P*S(x/t1,i+1)%P+t2%P*(t2%P-1)%P)%=P;
        }			
    return ans;
}
signed main(){ 
    n=read();block=sqrt(n);
    pre(block);
    register int i,j;
    register ll l,r;
    for(l=1;l<=n;l=r+1){
        r=n/(n/l);
        w[++m]=n/l;
        g[m]=w[m]%P*(w[m]+1)%P*(2ll*w[m]%P+1)%P*inv6%P-1;
        h[m]=w[m]%P*(w[m]+1)%P*inv2%P-1;
        if(w[m]<=block)id1[w[m]]=m;
        else id2[r]=m;
    }
    for(j=1;j<=cnt;j++)
        for(i=1;i<=m&&prim[j]*prim[j]<=w[i];i++){
            ll k=(w[i]/prim[j]<=block)?id1[w[i]/prim[j]]:id2[n/(w[i]/prim[j])];
            g[i]=(g[i]-prim[j]*prim[j]%P*(g[k]-sp[j-1]+P)%P+P)%P;
            h[i]=(h[i]-prim[j]*(h[k]-sp2[j-1]+P)%P+P)%P;
        }
    cout<<S(n,1)+1;
    return 0;
}
[SPOJ20173]DIVCNT2 - Counting Divisors (square)

题目链接

题意：求

\sum_{i=1}^{n}d(i^2)

solution1:暴力min25筛。考虑f(p)=3,f(p^k)=2k+1,我们只要设g(x)=1最后再乘上个3即可，复杂度O(T\frac{n^{\frac{3}{4}}}{logn})。

solution2:大概就是拆成

\sum_{i=1}^{n}{d*\mu^2}

然后对于\mu^2和d，大于n^{\frac{2}{3}}的现场O(\sqrt{n})暴力，小于n^{\frac{2}{3}}的事先预处理出来，复杂度就是O( n^\frac{2}{3})，证明可以参考杜教筛。

但是我们会发现因为某种玄妙原因，不管是否记忆化这个方法都跑不过min25筛，并且这种筛法空间是O(n^ \frac{2}{3})的，全面被吊打。

0%