[整理]min_25 筛

大概是一种牛逼的筛法,本质思想是 dp。
方便起见先规定本文的一些符号:
(p_i) 表示第 (i) 个质数,(p) 表示任意一个质数;
( ext{P}) 表示全体质数;
(min_i) 表示 (i) 的最小质因数。
现在我们要求一个积性函数 (f) 的前缀和,这个 (f) 需要满足在 (p^k) 处为关于 (p) 的多项式,这样我们只考虑 (f(p)=p^k) 的情况(是一个完全积性函数)。
容易发现如果设 (S(n,j)) 为最小质因数大于 (p_j)(f) 之和,我们要求的就是 (S(n,0)),而下面可以看到这个 (S) 是更加好求的。

接下来就是比较牛逼的了,我们设一个 dp 数组 (g_{n,j}) 表示 (sumlimits_{i=1}^n[iin ext{P}lormin_i>p_j]i^k)
显然 (p_{j-1}^2>n) 时这个东西就等价于 (g_{n,j-1}),下面考虑其他情况。
容易发现比起 (g_{n,j-1}) 少的部分就是最小质因数恰好为 (p_j) 的合数。利用完全积性函数的性质,可以提一个 (p_j^k) 出来,剩下的就是 (g_{lfloor n/p_j floor,j-1})
吗?我们多减了一些质数,应该加回来一个 (g_{p_{j-1},j-1})

求出了 (g),我们尝试寻找 (g)(S) 之间的关系。
显然,(S(n,j)) 可以按照类似的处理方式提取质因数后递归地计算:(S(n,j)=g_{n, ext{P}}-sumlimits_{i=1}^{j-1}f(p_i)+sumlimits_{i>j,p_i^xle n}f(p_i^x)(S(lfloor n/p_i^x floor,i)+[x>1]))

然后因为变量都是 (n) 整除一个数,可以离散化后做到空间 (sqrt{n}) 级别。
非常丑的洛谷模板题代码:

const int N=400010,p=1e9+7;
int n,B,g1[N],g2[N],id1[N],id2[N],inv3;
bool np[N];int pri[N],sp1[N],sp2[N],pri0;
il int Pow(int a,int b=p-2){
  int res=1;
  for(;b;a=(LL)a*a%p,b>>=1)if(b&1)res=(LL)res*a%p;
  return res;
}
il void Init(){
  np[1]=1;
  for(int i=2;i<N;i++){
    if(!np[i]){
      pri[++pri0]=i;
      sp1[pri0]=(sp1[pri0-1]+i)%p;
      sp2[pri0]=(sp2[pri0-1]+i*i%p)%p;
    }
    for(int j=1;j<=pri0&&i*pri[j]<N;j++){
      np[i*pri[j]]=1;
      if(i%pri[j]==0)break;
    }
  }
}
int w[N],tot;
il void DP(){
  //j=0 时就是求所有数的和
  for(int l=1,r;l<=n;l=r+1){
    r=n/(n/l);
    w[++tot]=n/l,g1[tot]=w[tot]%p;
    g2[tot]=g1[tot]*(g1[tot]+1)/2%p*(2*g1[tot]+1)%p*inv3%p;
    g1[tot]=(g1[tot]*(g1[tot]+1)/2%p-1+p)%p;
    g2[tot]=(g2[tot]-1+p)%p;
    if(n/l<=B)id1[n/l]=tot;
    else id2[n/(n/l)]=tot;
  }
  for(int i=1;i<=pri0;i++){
    for(int j=1;j<=tot&&pri[i]*pri[i]<=w[j];j++){
      int qwq=w[j]/pri[i];
      if(qwq<=B)qwq=id1[qwq];
      else qwq=id2[n/qwq];
      g1[j]=(g1[j]-pri[i]*(g1[qwq]-sp1[i-1]+p)%p+p)%p;
      g2[j]=(g2[j]-pri[i]*pri[i]%p*
      (g2[qwq]-sp2[i-1]+p)%p+p)%p;
    }
  }
}
int S(int x,int j){
  if(pri[j]>=x)return 0;
  int qwq=x<=B?id1[x]:id2[n/x];
  int res=(g2[qwq]-g1[qwq]-(sp2[j]-sp1[j]+p)%p+p+p)%p;
  for(int i=j+1;i<=pri0&&pri[i]*pri[i]<=x;i++){
    int pw=pri[i];
    for(int k=1;pw<=x;k++,pw=pw*pri[i]){
      int fuck=pw%p;
      (res+=fuck*(fuck-1)%p*(S(x/pw,i)+(k>1))%p)%=p;
    }
  }
  return res;
}
signed main(){
  Read(n),inv3=Pow(3),B=sqrt(n);
  Init(),DP();
  printf("%lld
",(S(n,0)+1)%p);
  KafuuChino HotoKokoa
}
内容来自_ajthreac_的博客(https://www.cnblogs.com/juruoajh/),未经允许,不得转载。
原文地址:https://www.cnblogs.com/juruoajh/p/15518626.html