洲阁筛

用上线性筛,我们可以用O(n)O(n)的线性时间求一个积性函数的点值;如果一个函数和另一个容易求前缀和的函数的狄利克雷卷积的前缀和也很好求,用上杜教筛,我们可以在O(n23)O(n^\frac{2}{3})的时间里求这个函数的前缀和。但是在积性这样一个限制下,仍有十分广阔的空间我们没有做研究。洲阁筛就是一个能在O(n34logn)O(\frac{n^\frac{3}{4}}{\log n})时间内求出很多积性函数的前缀和的算法。新人入坑推荐阅读国家集训队2016任之洲的《积性函数求和的几种方法》。

洲阁筛对它要求前缀和的函数f(x)f(x)有一些限制条件。首先f(x)f(x)必须是积性;其次pP,cN,f(pc)\forall p\in P,c\in N,f(p^c)能快速计算出来,pP,f(p)\forall p\in P,f(p)为一个有限次数多项式。满足这些条件的函数一定可以得到nN,f(n)=i=1kf(pici)\forall n\in N^*,f(n)=\prod_ {i=1}^{k}f(p_ i^{c_ i})。另外,由于积性要求,f(1)=0f(1)=011,由于f(1)=0f(1)=0的情况极为困难我完全不会:),本文只讨论f(1)=1f(1)=1的情况。

我们的目标是在低于线性时间复杂度下求:

i=1nf(i)\sum_ {i=1}^nf(i)

考虑把ii分成两类,一类是含有>n>\sqrt{n}的质因子,一类是不含的。

=i=1p>n,pinf(i)(1+p>npPnif(p))=\sum_ {\begin{matrix} i=1\\ \nexists p>\sqrt{n},p|i \end{matrix}}^n f(i)\left(1+\sum_ {\begin{matrix} p>\sqrt{n}\\ p\in P \end{matrix}}^{\left\lfloor\frac{n}{i}\right\rfloor}f(p)\right)

公式能这样写,是因为>n>\sqrt{n}的质因子最多一个。观察发现在i>ni>\sqrt{n}的时候最后一个求和项恒等于00。如果我们把ini\leq\sqrt{n}f(i)f(i)先预处理出来(看最终复杂度的话,这不会成为瓶颈),我们要求的东西就分为两部分:

i=1p>n,pinf(i)\sum_ {\begin{matrix} i=1\\ \nexists p>\sqrt{n},p|i \end{matrix}}^n f(i)

以及对于任意ini\leq\sqrt{n}要求出

p>npPnif(p)\sum_ {\begin{matrix} p>\sqrt{n}\\ p\in P \end{matrix}}^{\left\lfloor\frac{n}{i}\right\rfloor} f(p)

来吧。

算法的回归:动态规划

Gk[i][j]G_ k[i][j][1,j][1,j]中与前ii个质数互质的数的kk次幂和。

Gk[i][j]=Gk[i1][j]pikGk[i1][jpi]G_ k[i][j]=G_ k[i-1][j]-p_ i^kG_ k\left[i-1\right]\left[\left\lfloor\frac{j}{p_ i}\right\rfloor\right]

这个公式渲染得比较奇怪。

对于n<in\sqrt{n}<i\leq nii要么是质数,要么含有不超过n\sqrt{n}的质因数。对比我们要求的第二个式子,发现:

p>npPnif(p)=p>npPnij=0Kajpk=j=0Kajp>npPnipk=j=0KajGk[P(n)][ni]\begin{matrix} &\sum_ {\begin{matrix} p>\sqrt{n}\\ p\in P \end{matrix}}^{\left\lfloor\frac{n}{i}\right\rfloor} f(p)\\ = & \sum_ {\begin{matrix} p>\sqrt{n}\\ p\in P \end{matrix}}^{\left\lfloor\frac{n}{i}\right\rfloor}\sum_ {j=0}^Ka_ j\cdot p^k\\ = & \sum_ {j=0}^Ka_ j\sum_ {\begin{matrix} p>\sqrt{n}\\ p\in P \end{matrix}}^{\left\lfloor\frac{n}{i}\right\rfloor}p^k\\ = & \sum_ {j=0}^{K}a_ j\cdot G_ k\left[P(\sqrt{n})\right]\left[\left\lfloor\frac{n}{i}\right\rfloor\right] \end{matrix}

其中P(n)P(n)表示n\leq n的质数的数量。只要能得到了所有的GkG_ k,我们就可以在O(K)O(K)的时间内完成这一部分的计算。

要想清楚几个显然的结论:

  1. iN,Gk[i][0]=0\forall i\in N,G_ k[i][0]=0
  2. iN,Gk[i][1]=1\forall i\in N,G_ k[i][1]=1
  3. pi+1>jp_ {i+1}>jGk[i][j]=1G_ k[i][j]=1

如果暴力递推,我们枚举第二维jj。对于j<nj<\sqrt{n}的部分,不考虑pi>jp_ i>j的情况(根据第三条结论显然为11),用质数密度O(1logn)O(\frac{1}{\log n})估计,复杂度为O(jlogj)O(\frac{j}{\log j});对于j>nj>\sqrt{n}的部分,第一维全部转移,复杂度O(nlogn)O(\frac{\sqrt{n}}{\log\sqrt{n}})。算上一开始枚举的ii,近似估计复杂度为:

j=1nO(jlogj)+i=1nO(nlogn)O(nlogn)\sum_ {j=1}^{\sqrt{n}}O(\frac{j}{\log j})+\sum_ {i=1}^{\sqrt{n}}O(\frac{\sqrt{n}}{\log\sqrt{n}})\approx O(\frac{n}{\log n})

虽然已经比线性快,但是常数巨大,而且logn\log n是很小的东西,所以我们还要继续考虑优化。

从上面显然结论的三条里面推导,其实我们可以得到更多的简化。联系递推式里的Gk[i1][jpi]G_ k\left[i-1\right]\left[\left\lfloor\frac{j}{p_ i}\right\rfloor\right],如果p(i1)+1>jpip_ {(i-1)+1}>\left\lfloor\frac{j}{p_ i}\right\rfloorj<pi2j<p_ i^2,这个项会变成一。归纳起来:

Gk[i][j]=Gk[i1][j]+{pikGk[i1][jpi](jpi2)pik(pij<pi2)0(j<pi)G_ k[i][j]=G_ k[i-1][j]+\left\{\begin{matrix} -p_ i^k\cdot G_ k\left[i-1\right]\left[\left\lfloor\frac{j}{p_ i}\right\rfloor\right] & (j\geq p_ i^2)\\ -p_ i^k & (p_ i\leq j<p_ i^2)\\ 0 & (j<p_ i) \end{matrix}\right.

因为pinp_ i\leq\sqrt{n},所以pikp_ i^k的前缀和很容易O(kn)O(k\sqrt{n})预处理出来,O(1)O(1)查询。我们计算量堆积在了j>pi2j>p_ i^2这种情况上,其他情况都可以在调用Gk[i][n]G_ k[i][n]的时候累加几个质数幂的和,再归为第一种情况。而第一种情况的条件为jpi2j\geq p_ i^2,也就是对于第二维jj,我们只需要考虑第一维j\leq\sqrt{j}O(jlogj)O(\frac{\sqrt{j}}{\log\sqrt{j}})个质数。分块后总时间复杂度约为

j=1nO(jlogj)+i=1nO(nilogni)i=1nO(nilogni)O(0nnxdxlogn)=O(n120nx12dxlogn)=O(n12(n)12logn)=O(n34logn)\begin{matrix} & \sum_ {j=1}^{\sqrt{n}}O(\frac{\sqrt{j}}{\log j})+\sum_ {i=1}^{\sqrt{n}}O(\frac{\sqrt{\frac{n}{i}}}{\log\frac{n}{i}})\\ \leq & \sum_ {i=1}^{\sqrt{n}}O(\frac{\sqrt{\frac{n}{i}}}{\log\frac{n}{i}})\\ \approx & O(\frac{\int_{0}^{\sqrt{n}}\sqrt{\frac{n}{x}}\mathrm{d} x}{\log n})\\ = & O(\frac{n^{\frac{1}{2}}\int_{0}^{\sqrt{n}}x^{-\frac{1}{2}}\mathrm{d} x}{\log n})\\ = & O(\frac{n^{\frac{1}{2}}(\sqrt{n})^{\frac{1}{2}}}{\log n})\\ = & O(\frac{n^\frac{3}{4}}{\log n}) \end{matrix}

神奇的复杂度出现了!上界卡的紧紧的!有没有像分块大小取最优的感觉?原来除的一个log\log是因为质数密度的问题。

在计算微积分的过程中,因为OO记号省了很多常数的运算,甚是方便。

故技重施

在第一项中,我们要求的是

i=1p>n,pinf(i)\sum_ {\begin{matrix} i=1\\ \nexists p>\sqrt{n},p|i \end{matrix}}^n f(i)

F[i][j]F[i][j][1,j][1,j]中仅由后ii个质数(相对于n\sqrt{n})组成的数的ff函数值的和。每次新增一个质数pip_ i,枚举它在每个数中出现的次数。

F[i][j]=c=0logpijf(pic)F[i1][jpic]F[i][j]=\sum_ {c=0}^{\left\lfloor\log_ {p_ i}j\right\rfloor}f(p_ i^c)\cdot F[i-1]\left[\left\lfloor\frac{j}{p_ i^c}\right\rfloor\right]

这里终于用了一次积性。

类似的显然结论:

  • iN,F[i][0]=0\forall i\in N,F[i][0]=0
  • jN,F[0][j]=1\forall j\in N,F[0][j]=1
  • pi>jp_ i>j,则F[i][j]=1F[i][j]=1

结合第三条显然结论,暴力复杂度一样是O(nlogn)O(\frac{n}{\log n})的。那就来一波毛爷爷的神优化:若pi2>jp_ i^2>j,因为pip_ i是当前最小的质数,所以在[pi,min(j,n)][p_ i,\min(j,\sqrt{n})]内只有质数是要被统计进来的。

或者用类似GkG_ k的方法更形式化地证明一下,当j<pi2j<p_ i^2时,jpi<pi<pi1\left\lfloor\frac{j}{p_ i}\right\rfloor<p_ i<p_ {i-1},依然由显然结论三可得

F[i][j]=c=01f(pic)F[i1][jpic]=F[i1][j]+f(pi)F[i1][jpi]=f(1)F[i1][j]+f(pi)\begin{matrix} F[i][j] & = & \sum_ {c=0}^{1}f(p_ i^c)\cdot F[i-1]\left[\left\lfloor\frac{j}{p_ i^c}\right\rfloor\right]\\ & = & F[i-1][j]+f(p_ i)\cdot F[i-1]\left[\left\lfloor\frac{j}{p_ i}\right\rfloor\right]\\ & = & f(1)\cdot F[i-1][j]+f(p_ i) \end{matrix}

因此,

F[i][j]={c=0logpijf(pic)F[i1][jpic](jpi2)F[i1][j]+f(pi)(pij<pi2)1(j<pi)F[i][j]=\left\{\begin{matrix} \sum_ {c=0}^{\left\lfloor\log_ {p_ i}j\right\rfloor}f(p_ i^c)\cdot F[i-1][\frac{j}{p_ i^c}] & (j\geq p_ i^2)\\ F[i-1][j]+f(p_ i) & (p_ i\leq j<p_ i^2)\\ 1 & (j<p_ i) \end{matrix}\right.

预处理出质数的函数值的前缀和,很容易就能搞出j<pi2j<p_ i^2的情况。和GkG_ k一样十分类似,再推一波复杂度,也是O(n34logn)O(\frac{n^\frac{3}{4}}{\log n}),我就省略了。因为要枚举cc,常数估计会大一些。

于是

i=1p>n,pinf(i)=F[P(n)][n]\sum_ {\begin{matrix} i=1\\ \nexists p>\sqrt{n},p|i \end{matrix}}^n f(i)=F[P(\sqrt{n})][n]

实现细节

  1. 保证n\sqrt{n}的取值rr满足r2nr^2\geq n
  2. 线性筛出前P(n)P(\sqrt{n})个质数,然后计算他们这些质数的kk次幂前缀和,以及他们的kk次幂函数值的前缀和。
  3. 计算GkG_ k
    • 虽然我们在时间复杂度分析中第一层枚举的是jj,但是实际操作的时候可以第一层枚举ii,时间复杂度不变;第二层类似背包倒序枚举jj,在j<pi2j<p_ i^2的时候停止。这样方便开滚动第一维,空间复杂度O(n)O(\sqrt{n})
    • GkG_ k滚动数组的值分两部分存储:一部分是n\leq\sqrt{n}的,记录j[1,n]j\in[1,\sqrt{n}];一部分是>n>\sqrt{n}的,记录j{ni1i<n}j\in\left\{\left\lfloor\frac{n}{i}\right\rfloor|1\leq i<\sqrt{n}\right\}。由向下取整的性质,我们能保证所有可能出现的值都在这些范围内。可以开一个struct来封装一下map操作。
      • 这一步我请教了同学才想起来nab=nab\left\lfloor\frac{\left\lfloor\frac{n}{a}\right\rfloor}{b}\right\rfloor=\left\lfloor\frac{n}{ab}\right\rfloor会使得无论nn除了多少次数,只要没有加减,都能一定保证在上述范围之内。我原先的思路要求预处理完还要DFS展开,复杂度远不止这个数量级。这样比我原先想的方案不知道简单到哪里去了!
    • 我们需要储存一个最后转移到该jj的时间ii'',在调用Gk[i][j]G_ k[i][j]的时候,找到第一个pijp_ {i'}\leq j,把Gk[??][j]G_ k[??][j]的值i=i+1ipik-\sum_ {i=i''+1}^{i'}p_ i^k就是正确答案了。????是滚动的意思。
  4. 计算FF:和GkG_ k几乎一模一样的。
  5. 合并答案的时候,注意GkG_ k要乘上f(x)f(x)中的多项式系数和f(i)f(i)

总结

洲阁筛主要通过积性和质数的一些性质,将[1,n][1,n]中的数分类讨论,由于它们共同的性质大大简化了动态规划的转移和最终答案的计算,思路具有“一定的启发性”。但是我比较弱,没有能想到什么推广。

问我从这个过程里面学习到了什么?

  • 分类讨论后强行简化的思想;
  • 分块后讨论的数的范围通常是O(n)O(\sqrt{n})的;
  • Latex数学公式等号可以用\begin{matrix}\end{matrix}对齐。
avatar
Kerry Su