[模板]Min_25筛

时间:2023-03-08 21:57:29

用途

快速($O(\frac{n^{3/4}}{logn})$)地计算一些函数f的前缀和,以及(作为中间结果的)只计算质数的前缀和

一般要求$f(p)$是积性函数,$f(p)$是多项式的形式,且$f(p^k)$可以快速计算

做法

首先考虑求出范围内的质数的取值的和

如果有$f(p)=\sum{a_ip^i}$

那么我们构造$h_i(x)=x^i$,不难发现$h_i$是完全积性的

就是说,把f在质数的时候的式子拆开,然后让它在不是质数的时候也成立

考虑求其中的一个h,接下来设$pri_j$是第j个质数

设$g(n,j)=\sum\limits_{i=1}^{n}[i \in Prime OR min\_divisor(i)>pri_j]h(i)$,即n范围内质数或者最小质因数$>pri_j$的h值之和

可以当成是埃氏筛法的过程,已经筛掉了前j个质数

那么$g(n,inf)$就是要求的质数的和,g(n,0)就是所有数(1以外的)的h的和

(关于计算g(n,0),次数要是小的话可以直接代公式,大的话就要斯特林数或者别的什么来求自然数幂和了..不过别忘了减掉1)

考虑转移,有:

$g(n,j)=g(n,j-1) , pri_j*pri_j>n$

$g(n,j)=g(n,j-1)-h(pri_j)(g(\lfloor\frac{n}{pri_j}\rfloor,j-1)-\sum\limits_{i=1}^{j-1}h(pri_i)) , pri_j*pri_j<=n$

考虑$pri_j*pri_j<=n$的情况,它在从$g(n,j-1)$转移过来的时候,需要减掉那些最小质因子是$pri_j$的,而h又是完全积性的,所以可以直接乘。但这样减会多减掉那些比$pri_j$小的质数乘$pri_j$,需要再加回来

于是我们就求出来质数的h的和啦!然而这有什么卵用呢..

设$S(n,j)=\sum\limits_{i=1}^{n}[min\_divisor(i)>=pri_j]f(i)$,那么$S(N,1)+f(1)$就是要求的答案

考虑转移,有$S(n,j)=\sum a_ig_i(n,inf) - \sum\limits_{i=1}^{j-1}f(pri_i) + \sum\limits_{k=j}^{pri_k*pri_k<=n}\sum\limits_{e=1}^{pri_k^{e+1}<=n}f(pri_k^e)S(\lfloor \frac{n}{pri_k^e} \rfloor,k+1)+f(pri_k^{e+1})$

就是说,既然我们已经求出了质数的情况,那就先把符合要求的质数加上,然后再加上合数,枚举它的质因子然后都除下去就好了(因为这里已经没有完全积性了,所以一次要把某个质因子都除掉),因为S里是没有1的,所以还要额外加上$p^{e+1}$的情况

例题

loj6053 简单的函数

注意到除了f(2)=3以外,f(p)=p-1

于是拆成p和1去做即可。最后算S的时候,如果j<=1,那再加个2

 #include<bits/stdc++.h>
#include<tr1/unordered_map>
#define CLR(a,x) memset(a,x,sizeof(a))
#define MP make_pair
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> pa;
const int maxn=1e6+,P=1e9+; inline ll rd(){
ll x=;char c=getchar();int neg=;
while(c<''||c>''){if(c=='-') neg=-;c=getchar();}
while(c>=''&&c<='') x=x*+c-'',c=getchar();
return x*neg;
} int pri[maxn],sump[maxn],cnt,tot,sqrtn;
bool np[maxn];
ll N,w[maxn];
tr1::unordered_map<ll,int> id;
int g[maxn],h[maxn]; inline int S(ll n,int j){
if(n<=||pri[j]>n) return ;
int re=(0ll+g[id[n]]-h[id[n]]-sump[j-]+j-)%P;
if(j==) re=(re+)%P;
for(int i=j;i<=cnt&&1ll*pri[i]*pri[i]<=n;i++){
ll tmp=pri[i];
for(int k=;tmp*pri[i]<=n;k++){
re=(re+1ll*S(n/tmp,i+)*(pri[i]^k)+(pri[i]^(k+)))%P;
tmp=tmp*pri[i];
}
}return re;
} int main(){
//freopen("","r",stdin);
ll i,j;
N=rd();sqrtn=sqrt(N);
for(i=;i<=sqrtn;i++){
if(!np[i]){
pri[++cnt]=i;
sump[cnt]=(sump[cnt-]+i)%P;
}
for(j=;j<=cnt&&pri[j]*i<=sqrtn;j++){
np[pri[j]*i]=;
if(i%pri[j]==) break;
}
}
for(i=;i<=N;i=j+){
j=N/(N/i);
w[++tot]=N/i;id[w[tot]]=tot;
g[tot]=1ll*(w[tot]+)%P*((w[tot]-)%P)%P;
if(g[tot]&) g[tot]+=P;
g[tot]=g[tot]/%P;
h[tot]=(w[tot]-)%P;
}
for(j=;j<=cnt;j++){
for(i=;i<=tot&&1ll*pri[j]*pri[j]<=w[i];i++){
g[i]=(g[i]-1ll*pri[j]*(g[id[w[i]/pri[j]]]-sump[j-]))%P;
h[i]=(0ll+h[i]-(h[id[w[i]/pri[j]]]-j+))%P;
}
}
printf("%d\n",((S(N,)+)%P+P)%P); return ;
}