自动搬运

查看原文

来自洛谷，原作者为

	command_block 众水皆昂首，饮月唯我一。

搬运于2025-08-24 22:11:19，当前版本为作者最后更新于2020-01-27 16:30:38，作者可能在搬运后再次修改，您可在原文处查看最新版

自动搬运只会搬运当前题目点赞数最高的题解，您可前往洛谷题解查看更多

以下是正文

---

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

• 引理 : 合数$m$必有一个$\sqrt{m}$以内的因子。

因此,实现埃式筛的时候,只需要利用$\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$处的取值,有引理: $\lfloor\frac{N}{ab}\rfloor=\left\lfloor\frac{\lfloor\frac{N}{a}\rfloor}{b}\right\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^{2/3}}{logn})$的做法。

---

这里只给出$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);
}