min-25筛学习笔记

Min_25筛简介

( ext{min_25})筛是一种处理一类积性函数前缀和的算法。

其中这类函数(f(x))要满足(sum_{i=1}^{n}[iin prime]cdot f(i))可以被(sum_{i=1}^{n}[iin prime]cdot i^k)简单表示或者快速计算，其中(k)为较小的常数。

时间复杂度好像是(O(frac{n^{0.75}}{log n}))，不过据说被证伪了...也有人说是(O(n^{1-epsilon}))，反正可以做(n=1e10)，貌似跑的比杜教筛快。

空间复杂度(O(sqrt{n}))。

part 1

首先我们要处理的问题是(sum_{i=1}^{n}[iin prime]cdot f(i))，由于函数特性转化为了求(sum_{i=1}^{n}[iin prime]cdot i^k)。

设(g(n,i)=sum_{i=1}^{n}[iin prime or t(i)>p_i]cdot i^k)，其中(t(x))表示(x)的最小质因子，(p_i)表示第(i)个质数。

直观理解就是埃拉托斯特尼筛法进行了(i)轮之后(1sim n)剩下的数的(k)次方之和。

显然当(p_i^2>n)时(g(n,i)=g(n,i-1))，下面只考虑(p_i^2leqslant n)。

我们考虑(g(n,i))如何从(g(*,i-1))转移过来，也就是考虑第(i)轮筛掉了那些数。

那么可以得到式子：

[g(n,i)=g(n,i-1)-p_i^kcdot left(g(lfloorfrac{n}{p_i} floor,i-1)-s_{i-1} ight) ]

其中(s_i=sum_{j=1}^{i}p_i^k)。

解释一下这个式子，我们考虑筛掉了哪些数，然后减掉就好了，筛掉的那些数最小的质因子显然是(p_i)，而正好(g(lfloorfrac{n}{p_i} floor,i-1))代表最小质因子(geqslant p_i)的数以及质数，我们把质数减掉，拼起来就是(1sim n)。

那么写在一起就是：

[g(n,i)= egin{cases} g(n,i-1)&p_i^2>n\ g(n,i-1)-p_i^kcdot left(g(lfloordfrac{n}{p_i} floor,i-1)-s_{i-1} ight)&p_i^2leqslant n end{cases} ]

根据常识可以知道第一维只有(2sqrt n)个取值，离散化一下然后递推就好了。

注意到这里我们顺便求出了对于每个(x),(sum_{i=1}^{lfloor n/x floor}[iin prime]cdot i^k)的值。

part 2

我们想算出(f)的前缀和，同样的我们设(h(n,i))如下：

[h(n,i)=sum_{j=1}^{n}[t(j)geqslant p_i]cdot f(j) ]

那么(h(n,1)+f(1))就是答案，注意我们(h)不会把(1)算进去。

首先还是很显然的，当(p_i>n)时，(h(n,i)=0)，同样下面只考虑(p_ileqslant n)的情况。

我们可以通过(g)很简单的算出质数的贡献，即：

[sum_{j=1}^{n}[jin prime]cdot f(j)-sum_{j=1}^{i-1}f(p_j) ]

这些都可以预处理。

然后我们考虑暴力地枚举最小的质因子是多少以及用了多少个，即：

[sum_{j=i}^{p_j^2leqslant n}sum_{k=1}^{p_j^{k+1}leqslant n}left(f(p_j^k)cdot h(lfloorfrac{n}{p_j^k} floor,j+1)+f(p_j^{k+1}) ight) ]

解释下后面的(f(p_j^{k+1}))，这是因为(h)没有把(1)算进去，所以会漏算质数的次幂，直接加上就行了。

综合起来就是：

[h(n,i)=egin{cases} 0&p_i>n\ sum_{j=1}^{n}[jin prime]cdot f(j)-sum_{j=1}^{i-1}f(p_j)+sum_{j=i}^{p_j^2leqslant n}sum_{k=1}^{p_j^{k+1}leqslant n}left(f(p_j^k)cdot h(lfloorfrac{n}{p_j^k} floor,j+1)+f(p_j^{k+1}) ight)&p_ileqslant n end{cases} ]

写成这样求和符号好像有点丑...不管它

直接暴力递归算就好了，递推好像还慢些。

code

[LOJ6053] 简单的函数

(f(p)=p-1)，(p=2)的情况特殊处理下就好了。

#include<bits/stdc++.h>
using namespace std;

void read(int &x) {
    x=0;int f=1;char ch=getchar();
    for(;!isdigit(ch);ch=getchar()) if(ch=='-') f=-f;
    for(;isdigit(ch);ch=getchar()) x=x*10+ch-'0';x*=f;
}
 
void print(int x) {
    if(x<0) putchar('-'),x=-x;
    if(!x) return ;print(x/10),putchar(x%10+48);
}
void write(int x) {if(!x) putchar('0');else print(x);putchar('
');}

#define lf double
#define ll long long 

#define pii pair<int,int >
#define vec vector<int >

#define pb push_back
#define mp make_pair
#define fr first
#define sc second

#define FOR(i,l,r) for(int i=l,i##_r=r;i<=i##_r;i++) 

const int maxn = 1e6+10;
const int inf = 1e9;
const lf eps = 1e-8;
const int mod = 1e9+7;
const int inv2 = 5e8+4;

int add(int x,int y) {return x+y>=mod?x+y-mod:x+y;}
int del(int x,int y) {return x-y<0?x-y+mod:x-y;}
int mul(int x,int y) {return 1ll*x*y-1ll*x*y/mod*mod;}

ll n,w[maxn];
int b,tot,pri[maxn],vis[maxn],id[maxn],id2[maxn],m,g[maxn],h[maxn],f[maxn],p[maxn];

void sieve() {
    for(int i=2;i<=b<<1;i++) {
        if(!vis[i]) pri[++tot]=i,p[tot]=add(p[tot-1],pri[tot]);
        for(int j=1;j<=tot&&i*pri[j]<=b<<1;j++) {
            vis[i*pri[j]]=1;
            if(i%pri[j]==0) break;
        }
    }
}

int get_id(ll x) {return x<=b?id[x]:id2[n/x];}

void get_g() {
    for(ll l=1;l<=n;) {
        ll v=n/l,r=n/v;
        if(v<=b) id[v]=++m;else id2[n/v]=++m;  //利用整除分块来离散化
        w[m]=v,v%=mod;g[m]=mul(v+2,mul(v-1,inv2)),h[m]=del(v,1);l=r+1;  
    }
    for(int i=1;pri[i]<b;i++)
        for(int j=1;j<=m;j++) {
            if(1ll*pri[i]*pri[i]>w[j]) break;   //剪枝
            g[j]=del(g[j],mul(pri[i],del(g[get_id(w[j]/pri[i])],p[i-1]))); //g函数记录的是sum [iin pri] i，其中第二位滚动掉了
            h[j]=del(h[j],del(h[get_id(w[j]/pri[i])],i-1));     //h记录的是sum [iin pri] 1
        }
    for(int i=1;i<=m;i++) f[i]=del(g[i],h[i]);
}

int calc(ll x,int i) {
    if(x<=1||pri[i]>x) return 0;
    int res=del(f[get_id(x)],del(p[i-1],i-1));if(i==1) res=add(res,2);  //这里处理了f(2)=2+1的情况
    for(int j=i;1ll*pri[j]*pri[j]<=x;j++)
        for(ll k=1,e=pri[j];e*pri[j]<=x;e*=pri[j],k++)
            res=add(res,add(mul(pri[j]^k,calc(x/e,j+1)),pri[j]^(k+1))); //直接暴力算
    return res;
}

int main() {
    scanf("%lld",&n);b=ceil(sqrt(n));
    sieve();get_g();write(add(calc(n,1),1));  //注意加 1
    return 0;
}

---

洛谷P5325 【模板】Min_25筛

这和上题也是一样的，注意别爆( ext{long long})。

#include<bits/stdc++.h>
using namespace std;

void read(int &x) {
    x=0;int f=1;char ch=getchar();
    for(;!isdigit(ch);ch=getchar()) if(ch=='-') f=-f;
    for(;isdigit(ch);ch=getchar()) x=x*10+ch-'0';x*=f;
}
 
void print(int x) {
    if(x<0) putchar('-'),x=-x;
    if(!x) return ;print(x/10),putchar(x%10+48);
}
void write(int x) {if(!x) putchar('0');else print(x);putchar('
');}

#define lf double
#define ll long long 

#define pii pair<int,int >
#define vec vector<int >

#define pb push_back
#define mp make_pair
#define fr first
#define sc second

#define FOR(i,l,r) for(int i=l,i##_r=r;i<=i##_r;i++) 

const int maxn = 1e6+10;
const int inf = 1e9;
const lf eps = 1e-8;
const int mod = 1e9+7;
const int inv2 = 5e8+4;
const int inv6 = 166666668;

int add(int x,int y) {return x+y>=mod?x+y-mod:x+y;}
int del(int x,int y) {return x-y<0?x-y+mod:x-y;}
int mul(int x,int y) {return 1ll*x*y-1ll*x*y/mod*mod;}

ll n,w[maxn];
int b,tot,pri[maxn],vis[maxn],id[maxn],id2[maxn],m,g[maxn],h[maxn],f[maxn],p[maxn],p2[maxn];

void sieve() {
    for(int i=2;i<=b<<1;i++) {
        if(!vis[i]) pri[++tot]=i,p[tot]=add(p[tot-1],pri[tot]),p2[tot]=add(p2[tot-1],mul(i,i));
        for(int j=1;j<=tot&&i*pri[j]<=b<<1;j++) {
            vis[i*pri[j]]=1;
            if(i%pri[j]==0) break;
        }
    }
}

int get_id(ll x) {return x<=b?id[x]:id2[n/x];}

void get_g() {
    for(ll l=1;l<=n;) {
        ll v=n/l,r=n/v;
        if(v<=b) id[v]=++m;else id2[n/v]=++m;
        w[m]=v,v%=mod;
        g[m]=del(mul(mul(v,mul(v+1,(2*v+1)%mod)),inv6),1);
        h[m]=del(mul(v,mul(v+1,inv2)),1);l=r+1;
    }
    for(int i=1;pri[i]<b;i++)
        for(int j=1;j<=m;j++) {
            if(1ll*pri[i]*pri[i]>w[j]) break;
            g[j]=del(g[j],mul(mul(pri[i],pri[i]),del(g[get_id(w[j]/pri[i])],p2[i-1])));  //sum of square of prime
            h[j]=del(h[j],mul(pri[i],del(h[get_id(w[j]/pri[i])],p[i-1]))); //sum of prime
        }
    for(int i=1;i<=m;i++) f[i]=del(g[i],h[i]); 
}

int calc(ll x,int i) {
    if(x<=1||pri[i]>x) return 0;
    int res=del(f[get_id(x)],del(p2[i-1],p[i-1]));
    for(int j=i;1ll*pri[j]*pri[j]<=x;j++)
        for(ll k=1,e=pri[j],ee=e;e*pri[j]<=x;e*=pri[j],k++,ee=e%mod) //注意e过大可能会爆，先模一下就好了
            res=add(res,add(mul(mul(ee,ee-1),calc(x/e,j+1)),mul(mul(ee,pri[j]),del(mul(ee,pri[j]),1))));
    return res;
}

int main() {
    scanf("%lld",&n);b=ceil(sqrt(n));
    sieve();get_g();write(add(calc(n,1),1));
    return 0;
}