关于min_25筛的一些简单笔记
90nwyn
2020-10-07 13:32:32
设$f(i)$是一个积性函数,满足$f(p),p \in prime$可以用多项式表示,同时$f(p^k)$可以快速计算出来
考虑如何求解$\sum _{i=1}^nf(i)(1<n<10^{11})$
------------
另一个新的问题:对于积性函数$f(i)$,求$\sum_{i=1}^n[i \in prime] f(i)$
为了方便理解,问题变为求$\sum_{i=1}^n[i \in prime] i^k$
定义$prime_i$为第$i$个质数,$minp(i)$为$i$的最小质因子,令$minp(1)=inf$
对于$\sqrt n$范围内的前$i$个质数的$k$次方前缀和$sum_i=\sum _{j=1}^i prime_j^k$可以线性筛解决
设$g(n,i)=\sum_{j=1}^n[j \in prime $ $or$ $ minp(j)>prime_i]j^k$
直观地来看,$g(n,i)$就是在 $n$ 范围内的埃式筛算法进行第$i$轮后尚未被筛去的数的 $k$ 次方和
显然,$g(n,0)=\sum_{i=2}^ni^k$
思考从$g(n,i-1)$到$g(n,i)$的转移
- $prime_i^2>n$,在埃氏筛的第 $i$ 轮没有数字被筛掉,$g(n,i)=g(n,i-1)$
- $prime_i^2 \leq n$,考虑$prime_i$在$g(n,i-1)$中的贡献,$g(n,i)=g(n,i-1)-prime_i^k(g(\left \lfloor \frac{n}{prime_i}\right \rfloor,i-1)-sum_{i-1})$
令$cnt$为$\sqrt n$范围内的质数个数,有
$\sum_{i=1}^n[i \in prime] i^k=sum_{cnt}+g(n,cnt)$
时间复杂度 $O(\frac{n^{3/4}}{logn})$
------------
现在问题回归到求解$\sum _{i=1}^nf(i)$
定义$s(n,i)=\sum_{j=1}^n[minp(j) \geq prime_i] f(j)$
那么,$\sum _{i=1}^nf(i)=s(n,1)+f(1)$
问题最终化为如何求解$s(n,i)$
因为$\sum_{i=1}^n[i \in prime] f(i)$已经可以快速计算,那么不难得到答案中质数的贡献$\sum_{j=1}^n[j \in prime] f(j)-\sum_{j=1}^{i-1}f(prime_j)$
考虑合数的贡献,枚举每个合数的最小质因子以及最小质因子的次幂,由$f$为积性函数得
$\sum_{j=1}^{prime_j^2 \leq n}\sum_{k=1}^{prime_j^{k+1} \leq n} ( s( \left \lfloor \frac{n}{prime_j^k}\right \rfloor,j+1)f(prime_j^k)+f(prime_j^{k+1}))$
------------
[题目链接](https://vjudge.net/problem/LibreOJ-6053)
考虑$f(p),p \in prime$的值
若$p=2$,$f(p)=p-1$,反之$f(p)=p+1$
问题相当于min_25筛求出$n$以内的质数和与质数个数
------------
```cpp
#include<bits/stdc++.h>
using namespace std;
typedef long long int ll;
const int M=1e6+5,mod=1e9+7,inv2=500000004;
int prime[M],tot,vis[M],sqr,sum[M],id1[M],id2[M],m,g[M],h[M];
ll n,w[M];
void init(int siz)
{
for(int i=2;i<=siz;i++)
{
if(!vis[i])prime[++tot]=i,sum[tot]=(sum[tot-1]+i)%mod;
for(int j=1;j<=tot&&i*prime[j]<=siz;j++)
{
vis[i*prime[j]]=1;
if(i%prime[j]==0)break;
}
}
}
int getid(ll x){return (x<=sqr)?id1[x]:id2[n/x];}
ll s(ll x,int y)
{
if(x<=1||prime[y]>x)return 0;
int pos=getid(x);
ll res=(g[pos]-sum[y-1]-h[pos]+y-1+mod)%mod;
if(y==1)res+=2;
for(int i=y;i<=tot&&(ll)prime[i]*prime[i]<=x;i++)
{
ll p1=prime[i];
ll p2=p1*p1;
for(int e=1;p2<=x;e++,p1=p2,p2*=prime[i])
res=(res+s(x/p1,i+1)*(prime[i]^e)%mod+(prime[i]^(e+1))%mod)%mod;
}
return res;
}
int main()
{
scanf("%lld",&n);
sqr=sqrt(n);
init(sqr);
for(ll l=1,r=0;l<=n;l=r+1)
{
w[++m]=n/l;
r=n/w[m];
if(w[m]<=sqr)id1[w[m]]=m;
else id2[n/w[m]]=m;
g[m]=((2+w[m])%mod)*((w[m]-1)%mod)%mod*inv2%mod;
h[m]=(w[m]-1)%mod;
}
for(int j=1;j<=tot;j++)
for(int i=1;i<=m&&(ll)prime[j]*prime[j]<=w[i];i++)
{
int k=getid(w[i]/prime[j]);
g[i]=(g[i]-(ll)prime[j]*((g[k]-sum[j-1]+mod)%mod)+mod)%mod;
h[i]=(h[i]-h[k]+j-1+mod)%mod;
}
printf("%lld\n",(s(n,1)+1)%mod);
return 0;
}
```