题解 P6021 【【模板】质数前缀统计】

· · 题解

我们考虑借鉴埃式筛法的思想,逐渐使用质数将不合法的数筛去。

因此,实现埃式筛的时候,只需要利用\sqrt{N}以内的质数来筛除,就能正确地得到N以内所有的质数。

令h(n,k)表示埃式筛法枚举了前k个质数后,不超过n的,还剩下的数的c次方和,即

\sum\limits_{i=1}^n\large i^c\ \small[\text{i不含}\leq p_k\text{的素因子}]

我们以k从小到大的方式,筛去“最小素因子为p_k的数”,注意只筛到\sqrt{N}即可。

考虑每次减去被筛掉的数的贡献,可以得到递推式 :

h(n,k)=h(n,k-1)-p_k^c\begin{cases}{\big(h(\lfloor n/p_k\rfloor,k-1)-h(p_k-1,k-1)\big)}&(p^k\leq\sqrt{n})\\0&(p^k>\sqrt{n})\end{cases}

边界 : h(n,0)=\sum\limits_{i=2}^ni^k,此时只有1被筛去。

至于如何求解自然数幂和,不在本文的范围,可见 P4593 & CF622F , 使用拉格朗日插值这部分复杂度为O(k\sqrt{n})

我们必须不重不漏地选出最小素因子为p_k的数。

考虑埃式筛是怎么筛的 : 每次从已知的数范围内枚举x,如果其没有被筛掉,那么就能得到x没有小于p_k的素因子。

那么就能得到x*p_k最小素因子为p_k(不重),我们筛去x*p_k,这里能得到x\leq\lfloor n/p_k\rfloor。

不漏的性质显然。

这样子,我们去除h(\lfloor n/p_k\rfloor,k-1)的贡献。

注意到h(p_k-1,k-1)中的贡献都是质数,明显拥有小于p_k的素因子,这部分不应减去。

当p_k>\sqrt{n}时,有h(n,k)=h(n,k-1),即没有数被筛除,结合引理容易理解。

形式上就是 : \lfloor n/p_k\rfloor<p_k。

注意到我们在递推中需要利用\lfloor n/p_k\rfloor处的取值,有引理:

结合整除分块的结论,我们从$N$出发,除上一系列整数,只会得到$O(\sqrt{N})$种结果。 我们只用维护任意的$\lfloor N/d\rfloor$处的答案即可,这一共$O(\sqrt{N})$个。 接下来我们讨论实现方法。 - ① 如果采用暴力实现,在每个质数处使用上述递推式,复杂度为: 质数个数$O(\frac{\sqrt{N}}{logN})$,维护的值$O(\sqrt{N})$,总复杂度为$O(\frac{N}{logN})$,无法通过。 - ② 注意到某些值在递推中是不会改变的,在递推式中不断`+0`。 有$O(\frac{\sqrt{m}}{logm})$个质数在转移时会影响到$h(m,?)$地值,取最大的$\sqrt{n}$个$\lfloor N/d\rfloor$积分, 总复杂度 : $O(\sum\limits_{d=1}^{\sqrt{N}}\frac{\sqrt{N/d}}{log(N/d)})=O(\frac{N^{3/4}}{logN})$,这也就是洲阁筛第一部分的经典复杂度。 大多数题目里,做到这样已经够了,由于代码简单,常数小,**实战中最为常用**。 - ③ 当然,可以考虑利用朴素筛法预处理一部分,假设我们预处理$[1,K]$。 每次新加质数的时候都把$K$以内暴力筛除一遍,需要**树状数组**维护动态前缀和,复杂度$O(KlogK)

递推的运算次数O(\sum\limits_{d=1}^{N/K}\frac{\sqrt{N/d}}{logN})=O(\frac{N}{logN\sqrt{K}}),不过注意到树状数组需要一个O(logK),这部分复杂度为O(\frac{N}{\sqrt{K}})

总的复杂度就是O(\frac{N}{\sqrt{K}}+KlogK),取K=O((\frac{N}{logN})^{2/3})可得复杂度为O(N^{2/3}log^{1/3}N)\rightarrow O(N^{2/3+e})。

这还不够优秀,并且由于常数问题,跑过暴力需要精细实现。

考虑省掉递推运算中的部分树状数组操作。

注意到h(n,k)中如果p_k>\sqrt{n}的话函数值不再改变,此时我们直接存下对应的函数值即可,没必要再使用树状数组。

观察h(n,k)递推式中使用的函数值:

h(\lfloor n/p_k\rfloor,k-1)

此时,这个函数值仍在变化的条件是\sqrt{n/p_k}>p_k\rightarrow n>p_k^3

也就是说,在h(m,?)的求值中有O(\frac{\sqrt[3]{m}}{logm})个树状数组操作,复杂度即为O(\sqrt[3]{m})。

这部分总复杂度是O(\sum\limits_{d=1}^{N/K}\sqrt[3]{N/d})=O(\frac{N}{K^{2/3}})<O(\frac{N}{logN\sqrt{K}}),可不计。

h(p_k-1,k-1)

这时该函数值早已确定,前缀和即可。

这样子的总复杂度是O(\frac{N}{logN\sqrt{K}}+KlogK)取K=O(\frac{N^{2/3}}{log^{4/3}N})可得复杂度为O(\frac{N^{2/3}}{log^{1/3}N})\rightarrow O(N^{2/3-e})。

想要更快的话请见 : zzt巨神 :《一些特殊的数论函数求和问题》(2018国家候选队队论文)

里面提出了一种O((\frac{N}{logN})^{2/3})的做法,但是实现较为复杂。

这里只给出O(\frac{n^{3/4}}{logn}+k\sqrt{n})的代码实现。

使用了一个除法优化。

#include<algorithm>
#include<cstdio>
#include<cmath>
#define Limit 320000
typedef long long ll;
using namespace std;
int mod;
ll powM(ll a,int t)
{
  ll ans=1;
  while(t){
    if (t&1)ans=ans*a%mod;
    a=a*a%mod;
    t>>=1;
  }return ans;
}
int k,m;
ll lf[15],rf[15],ifac[15],y[15];
void Init()
{
  m=k+2;
  for (int i=1;i<=m;i++)
    y[i]=(y[i-1]+powM(i,k))%mod;
  ifac[0]=ifac[1]=1;
  for (int i=2;i<=m;i++)
    ifac[i]=ifac[mod%i]*(mod-mod/i)%mod;
  for (int i=2;i<=m;i++)
    ifac[i]=ifac[i]*ifac[i-1]%mod;
  lf[0]=1;
}
ll gPows(ll N)
{
  N%=mod;
  for (int i=1;i<=m;i++)
    lf[i]=lf[i-1]*(N-i)%mod;
  rf[m+1]=1;
  for (int i=m;i;i--)
    rf[i]=rf[i+1]*(N-i)%mod;
  ll ans=0;
  for (int i=1;i<=m;i++)
    ans=(ans+y[i]*lf[i-1]%mod*
         rf[i+1]%mod*ifac[i-1]%mod*
         ((m-i)&1 ? mod-ifac[m-i] : ifac[m-i]))%mod;
  return ans;
}
ll N,h0[Limit],h1[Limit];
int lim;
int main()
{
  scanf("%lld%d%d",&N,&k,&mod);
  Init();
  lim=sqrtl(N);
  for(int i=1;i<=lim;++i){ 
    h1[i]=gPows(N/i)-1; 
    h0[i]=gPows(i)-1;
  }
  for(int i=2;i<=lim;++i){
    h0[i]=(h0[i]+mod)%mod;
    if(h0[i]==h0[i-1])continue;
    ll x0=h0[i-1],r=(ll)i*i,p0=powM(i,k);
    int u=min((ll)lim,N/((ll)i*i)),uu=min(u,lim/i);
    for(int j=1;j<=uu;++j)
      h1[j]=(h1[j]-p0*(h1[j*i]-x0))%mod;
    ll t=N/i;
    if (t<(1ll<<31))
      for(int tt=t,j=uu+1;j<=u;++j)
        h1[j]=(h1[j]-p0*(h0[tt/j]-x0))%mod;
    else 
      for(int j=uu+1;j<=u;++j)
        h1[j]=(h1[j]-p0*(h0[t/j]-x0))%mod;
    for(int j=lim;j>=r;--j)
      h0[j]=(h0[j]-p0*(h0[j/i]-x0))%mod;
  }
  ll ans=0;
  for (int i=1;i<=lim;i++)
    ans=(ans+h1[i]*i%mod*i)%mod;
  printf("%lld\n",(ans+mod)%mod);
}