首页 > 技术文章 > [整理]min_25 筛

juruoajh 2021-11-06 22:15 原文

大概是一种牛逼的筛法,本质思想是 dp。
方便起见先规定本文的一些符号:
\(p_i\) 表示第 \(i\) 个质数,\(p\) 表示任意一个质数;
\(\text{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}\) 表示 \(\sum\limits_{i=1}^n[i\in\text{P}\lor\min_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\rfloor,j-1}\)
吗?我们多减了一些质数,应该加回来一个 \(g_{p_{j-1},j-1}\)

求出了 \(g\),我们尝试寻找 \(g\)\(S\) 之间的关系。
显然,\(S(n,j)\) 可以按照类似的处理方式提取质因数后递归地计算:\(S(n,j)=g_{n,\text{P}}-\sum\limits_{i=1}^{j-1}f(p_i)+\sum\limits_{i>j,p_i^x\le n}f(p_i^x)(S(\lfloor n/p_i^x\rfloor,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\n",(S(n,0)+1)%p);
  KafuuChino HotoKokoa
}

推荐阅读