UOJ Logo negiizhao的博客

博客

OI中常用数论函数求和法的简化陈述

2021-07-11 06:32:00 By negiizhao

如果还不会,或许也可以先从wc2018的原始课件或者 zzt blog(里面有一些小bug) 看一点基础知识.(upd:csdn 现在需要登录并关注才能看了)


upd:本文于 2024.01.30 凌晨 review 了一遍,干掉了很多坏的表述(感觉后半基本是重写了),原来写的感觉自己都要看不懂了!!1


从最基础的知识开始介绍吧,首先是 Dirichlet 双曲线法,对于数论函数 $f, g, h : \mathbb Z^+ \to R,f=g*h$,其中 $*$ 为 Dirichlet 卷积,$R$ 为交换环(并且有时假设较小的正整数在 $R$ 中可逆),那么有 $$ \sum _{n\leq x}f(n)=\sum _{n\leq x}\sum _{ab=n}g(a)h(b)=\sum _{a\leq {\sqrt {x}}}\sum _{b\leq {\frac {x}{a}}}g(a)h(b)+\sum _{b\leq {\sqrt {x}}}\sum _{a\leq {\frac {x}{b}}}g(a)h(b)-\sum _{a\leq {\sqrt {x}}}\sum _{b\leq {\sqrt {x}}}g(a)h(b) $$

从图像上理解,相当于是平面直角坐标系第一象限有一条双曲线 $\displaystyle y=\frac{n}{x}$ 也就是 $xy=n$,我们要对其下方(包括线上)的所有整点对 $(x, y)$,求 $g(x)h(y)$ 的和.我们对 $x \leq \sqrt n$ 的部分、$y \leq \sqrt n$ 的部分分别求和,再减去两者的公共部分.

james2.png

也许会有人更熟悉“整除分块”的做法,它和上式只相差一个 Abel 变换,我想这个做法或许会方便操作,常数更小一些,不过无论如何二者能力是差不多的.

无论是用哪个式子,我们都需要知道 $g$ 和 $h$ 在一些位置上的前缀和(这里当然认为我们对 $g$ 和 $h$ 的单点值并没有更多了解).这个位置的集合是我们非常熟悉的整除集合 $\displaystyle D_n=\left\{\left\lfloor\frac nm \right\rfloor \;\middle|\; m\in\mathbb Z^+ \right\}=\left\{0, 1, 2, \dots, \left\lfloor\sqrt n\right\rfloor-1, \left\lfloor\sqrt n\right\rfloor, \left\lfloor\frac n {\left\lfloor\sqrt n\right\rfloor}\right\rfloor, \left\lfloor\frac n {\left\lfloor\sqrt n\right\rfloor-1}\right\rfloor, \dots, \left\lfloor\frac n 2\right\rfloor, n\right\}$,这个集合的一些性质就不赘述了,可以参考 zzt blog.

可以看出来,这些位置的前缀和值是很重要的.我们记 $S(f) : D_n \to R$ 表示 $f$ 在整除集合位置的前缀和函数(符号中省去 $n$,因为我们这里考虑的 $n$ 总是固定的),那么上文说明我们可以由 $S(g)$ 和 $S(h)$ 得到 $g*h$ 的前 $n$ 项和.我们对此并不满意:也许 $g*h$ 只是计算过程的一小步,我们希望可以得到 $S(g*h)$.而这也不难做到,因为整除位置的整除集合是整除集合的子集,只需要直接对每个位置求一遍即可.经过计算容易得出这一过程的时间为 $\Theta(n^{3/4})$.

我们仍然不满意:这一过程还是稍微慢了一点.考虑对小于 $B$ 的位置暴力计算 Dirichlet 卷积求解,这样时间为 $O\left(B\log B+\sqrt n\sqrt{\frac nB}\right)$,取 $B=\Theta\left(\left(\frac{n}{\log n}\right)^{2/3}\right)$ 时为 $\Theta(n^{2/3}\log^{1/3} n)$.若数论函数为积性函数,可以通过记录 $B$ 以内函数的完整信息,增大空间开销使计算卷积的时间降为 $\Theta(B)$,此时取 $B=\Theta(n^{2/3})$ 时总时间为 $\Theta(n^{2/3})$.此外注意到这种方法亦可用于已知 $S(f*g)$ 和 $S(g)$ 反推 $S(f)$,即计算除法.由于一些历史原因,这种用法常常被称为“杜教筛”.

在一些特殊情况下(对,我的意思是下文会用到),两个函数的非零项密度为 $O(\log^{-1} n)$,此时时间为 $O\left(B\log B\log^{-2} n+\sqrt n\sqrt{\frac nB}\log^{-1} n\right)$,取 $B=\Theta(n^{2/3})$ 时总时间为 $O(n^{2/3}\log^{-1}n)$.


现在我们考虑适用于一类积性函数的求和法,它要求函数在素数处的取值是关于素数的常数次多项式.

为了方便以下不区分积性函数和它的 Dirichlet 生成函数,$f(x)$ 意思仍然和上文一样,在本文语境下对 Dirichlet 级数进行求值显然没有任何意义.

用 $f_p$ 表示对 $f$ 仅保留 $p$ 的幂处的值,其余值设为 0,积性函数无非是意味着 $\displaystyle f=\sum_{n=1}^\infty f(n)n^{-s}=\prod_p\left(1+\sum_{k=1}^\infty f(p^k)(p^k)^{-s}\right)=\prod_p f_p$.

我们希望对一个中间结果附加或取消 $f_p$ 的贡献,即由 $S(g)$ 得出 $S(g*f_p)$ 或 $S(g / f_p)$,这可以用 $O(\sqrt n\log_pn)$ 时间做到(对每个点的前缀和考虑如何贡献到它即可)

(比较重要的情况为仅有 $p$ 处有贡献,而更高幂次处没有贡献.直接计算乘法显然是 $O(\sqrt n)$ 时间,而除法自然是按乘法的过程倒推回去.)

(另一个常见的情况是对完全积性函数附加/取消贡献,完全积性意味着 $f_p$ 的逆很简单,为 $1-f_p(p)p^{-s}$,此时将乘或除以 $f_p$ 变为除以或乘 $f_p^{-1}$ 即可更轻松地进行操作.)

仅考虑大于 $\sqrt n$ 的素数的贡献 $\displaystyle S\left(\prod_{p>n^{1/2}}f_p\right)$(以下简称“大素数贡献”),此时任何 $pq$ 都会大于 $n$,因此只会在素数处非 0,于是可以将任意素数处取值为多项式的积性函数的大素数贡献拆解为 $id_k$ 的大素数贡献的线性组合.

于是我们有了个做法:先通过对 $id_k$ 取消小素数贡献得到大素数贡献,然后组合出想要求的积性函数的大素数贡献,然后再附加上小素数贡献即可.这样我们就得到了一个 $\displaystyle \sum_{p\le n^{1/2}}O(\sqrt n\log_pn)=\int_1^{n^{1/2}}O(\sqrt n\log_pn)\log^{-1}p\,\mathrm{d}p=O(n\log^{-1}n)$ 的做法.

但我们真的只能做这么好吗?考虑取消贡献时从小到大取消,附加贡献则相反,从大到小附加.当我们准备取消素数 $p$ 的贡献时,$p^2$ 以下只有 $1$ 和大于等于 $p$ 的素数处有值,这意味着这部分的值可以快速获取而无需进行维护,我们只需要考虑大于等于 $p^2$ 的部分受到了什么样的影响即可.既然如此,我们不如取消贡献时仍然保留素数处的值,这样 $p^2$ 以下的前缀和都无需再进行维护.以 $n^{1/4}$ 为分界,分别计算可得 $\displaystyle \sum_{p\le n^{1/4}}O(\sqrt n\log_pn)+\sum_{n^{1/4}< p\le n^{1/2}}O\left(\frac{n}{p^2}\log_pn\right)=O(n^{3/4}\log^{-1}n)$

(上面的分析里有很多 $\log_pn$ 项似乎使得分析变得不方便,其实容易发现它根本不影响结果.可以选定足够小的 $c$ 在 $n^c$ 处进行分段,小于等于 $n^c$ 的部分很小,可以直接放缩成 $\log n$ 而不影响最终结果,大于 $n^c$ 的部分 $\log_pn$ 为常数.亦可直接在算法中用 $O(n^{1/2}\log n)$ 时间(upd:我这部分写得比较晚,这里的 $O(n^{1/2}\log n)$ 的方法是注意到仅 powerful 部分非 0 的积性函数非常稀疏,直接枚举每个非 0 值,使用我下一篇 blog 所说“最优”乘法的做法中一行的贡献方法暴力算贡献即可)修改 powerful 部分的贡献(比如我们要从积性函数 $f$ 得到积性函数 $g$,若它们的素数处值都相同,则 $g/f$ 为仅有 powerful 部分非 0 的积性函数,可暴力算出 $S(g/f)$ 乘到 $S(f)$ 上.利用 powerful number 的想法,我最早能找到 利用powerful number求积性函数前缀和 - fjzzq2002(如果仅求单点值的话,修性质较好的 powerful 贡献其实能做到 $O(n^{2/5} \log^{-8/5} n)$,方法基本是下一篇 blog 思路的延续,从而根本无需乘上 $\log_pn$)

这种做法常被称为“洲阁筛”.其另有一个较小数据范围内速度较快的搜索变种,常被称为“min_25筛”,根据 OEIS 所说,“min_25筛”的复杂度大致为 $n \exp\left(-(1 + o(1))\sqrt{\log n \log \log n}\right)$.


这样的方法还是不够令人满意.我们以 $O(n\log^{-1}n)$ 的做法为框架,考虑再将小素数按 $n^{1/6}$ 分为 2 类,不妨重新将 $p\le n^{1/6}$ 称为小素数,而将 $n^{1/6}< p\le n^{1/2}$ 称为中素数.

假使我们已经得到中素数的贡献 $\displaystyle S\left(\prod_{n^{1/6}< p\le n^{1/2}}f_p\right)$,然后直接将小素数贡献附加上去得到 $\displaystyle S\left(\prod_{p\le n^{1/2}}f_p\right)$,仿照上文分析易知附加小素数需 $O(n^{2/3}\log^{-1}n)$ 时间.而计算出小中素数单独的贡献后,从整体中取消贡献或与大素数贡献合并仅需 $O(n^{1/2}\log n)$ 时间.

于是只需要考虑如何计算中素数的贡献.注意到 $\displaystyle \prod_{n^{1/6}< p\le n^{1/2}}f_p=\exp\sum_{n^{1/6}< p\le n^{1/2}}\log f_p$,而 $\exp$ 内部式子最低项 $6$ 次方已超出 $n$,对答案无影响.于是只需要计算其 2、3、4、5 次方.注意到涉及的函数的非零项密度皆为 $O(\log^{-1} n)$(证明可参考zzt blog),所以计算这部分总时间为 $O(n^{2/3}\log^{-1}n)$.

以上算法大概就是简化后的zzt求和法了()可以看出,通过调整小素数分类位置选取,可以进一步减少附加小素数的开销,而对计算中素数贡献的开销只有常数影响,除非将乘法时间降到 $\tilde O(n^{1/2})$.因此这种求和法所需时间实际上只和做乘法的时间有关.

暂时还没有进一步研究,也没考虑怎么卡常,因为↓

写了一份代码,虽然过了模板题,但好像还有些没弄明白的bug(我家乘法不结合.png),其中282行~294行是与题目本身有关的内容,此外都是模板,感觉代码量还是有点大的.

现在修好了,拆解了一些过程,现在287行~299行是与题目本身有关的内容.

修了一下,现在303行~315行是与题目本身有关的内容.

upd:仿佛整出来了个 $O(n^{2/3}\log^{-4/3}n)$ 甚至 $O(n^{2/3}\log^{-2}n)$(upd:很久之后发现算错了,其实是 $O(n^{2/3}\log^{-3/2}n)$)的做法,还有更快的算 $\pi(n)$ 的做法(x)(upd:被我鸽了 2.5 年后写在了这里

评论

暂无评论

发表评论

可以用@mike来提到mike这个用户,mike会被高亮显示。如果你真的想打“@”这个字符,请用“@@”。