广义积性函数的新筛法
Naszt
·
·
算法·理论
前言
就思想上来说，这个算法或许可以称作「欧拉乘积块筛筛法」。
※ 块筛指的是按整除分块求和，后文会再详细说明。
我觉得整个技巧比较新颖，而且用途较广，我想它命名为「Naszt 筛」。
其可以在 O\left(N^{2/3}\right) 或更快的时间复杂度内算出一类数论函数的前缀和问题。
这类数论函数需要满足可以写成以下 DGF 形式：
\tilde{F}(x)=\sum_{n\ge1}\frac{f(n)}{n^x}=\prod_{i\ge1}\left(1+A_i(x)\right)
其中 A_i(x) 是形如 A_i(x)=\sum_n a_{i,n}n^{-x} 的 DGF。
对于给定的 A_i(x)，下文提出了一种分析时间复杂度的方法。
要达到这个时间复杂度的条件比较宽松，积性函数是完全满足的。
不过对于大部分的积性函数，还是要计算「质数前缀块筛」，没必要用这个筛法。
所以这个筛法主要对能写成上述形式的广义积性函数有用。
例题
我们先从一个具体问题引入，给定 N，求：
\sum_x[\varphi(x)\le N]
我就是从这开始思考的，这题是我由上次的 C 题 P12470 改的。
但这题我找到了 原题，它的官方题解 用 min-25 筛做到了 O\left(\frac{N^{3/4}}{\log N}\right)。
用我的这个方法可以做到 O\left(N^{2/3}\right)。
推式子
该问题其实就是 欧拉函数的重数函数 前缀和。
\sum_{x}[\varphi(x)\le N]=\sum_{i=1}^{N}\sum_{x}[\varphi(x)=i]=\sum_{i=1}^{N}\psi(i)
要把这个函数写成 DGF，可以看看我上次的 C 题的 官方题解 的推导：
\tilde{\Psi}(x)=\sum_{n\ge1}\frac{\psi(n)}{n^x}=\sum_{n\ge1}\frac{1}{\varphi^x(n)}=\zeta(x)\tilde{F}(x)
\tilde{F}(x)=\sum_{n\ge1}\frac{f(n)}{n^x}=\prod_{p\in\mathcal{P}}\left(1+(p-1)^{-x}-p^{-x}\right)
于是我们只需计算 \tilde{F}(x) 的块筛：
\sum_{i=1}^{N}\psi(i)=\sum_{i=1}^{N}\sum_{d\mid i}f(d)=\sum_{i=1}^{N}S_f\left(\left\lfloor\frac{N}{i}\right\rfloor\right)
对于序列 a，其块筛指的是对于所有不同的 x=\left\lfloor\frac{N}{k}\right\rfloor，求出 S_a(x)={\rm sum}[1,n]=\sum_{i=1}^xa(i)。
根据 整除分块 的结论，我们可以得到只有 2\sqrt N 个块。
根据 杜教筛，我们可以在 O(N^{3/4}) 的时间复杂度计算出块筛卷积。
这里杜教筛做不到 \widetilde{O}(n^{2/3}) 是因为块筛不能提供足够的信息预处理。
不同于上面一般的定义，下文中我们都把块筛（块）定义为：
{\rm sum}\left[\left\lfloor\frac{N}{x+1}\right\rfloor+1,\left\lfloor\frac{N}{x}\right\rfloor\right]=S_f\left(\left\lfloor\frac{N}{x}\right\rfloor\right)-S_f\left(\left\lfloor\frac{N}{x+1}\right\rfloor\right)
相比于一般定义，这样做差分的好处是更新某一处的值只会影响一个块。
筛法
记第 i 个质数为 p_i，\tilde{F}(x) 可以拆成这么两部分：
\begin{aligned}
\tilde{F}(x)&=\tilde{A}(x)\times\tilde{B}(x)\\
&=\prod_{i=1}^{\pi(\sqrt N)}\left(1+(p_i-1)^{-x}-p_i^{-x}\right)\times\prod_{i>\pi(\sqrt N)}\left(1+(p_i-1)^{-x}-p_i^{-x}\right)
\end{aligned}
其中 \tilde{B}(x) 可以直接截断截断，于是可以改写为求和：
\tilde{B}(x)=1+\sum_{i>\pi(\sqrt N)}\left((p_i-1)^{-x}-p_i^{-x}\right)
于是对于 \tilde{B}(x)，当 x>\sqrt N 时，有 {\rm sum}[1,x]=1+[x+1\in\mathcal{P}].
检验 \sqrt N 次质数即可得出其块筛，用 Miller Rabin 素性检验即可。
假设已经得出了 \tilde{A}(x) 的块筛，由于 \tilde{B}(x) 在 [2,\sqrt N] 上没值，
显然可以在 O(\sqrt N\log N) 内完成二者的块筛卷积。
我们现在只用考虑计算 \tilde{A}(x) 的块筛。
每次卷上 (1+(p_i-1)^{-x}-p_i^{-x}) 时不用关心 i 的顺序，
于是我们考虑 p_{\pi(\sqrt N)}\to p_{1}，即为反向遍历质数。
设此时卷到了 p_i 项，我们需要先更新 f(p-1),f(p)。
由于此时区间 [2,p_{i+1}-2] 上系数为零，
剩下的需要更新的区间的范围是 [(p_i-1)(p_{i+1}-1),N]。
我们需要做的转移的区间可以近似的看作：
[p_i,N/p_i]\xrightarrow{\times p_i}[p_i^2,N]
从数字密度上来看，我们要均匀的将 [p_i,N/p_i] 转移到 [p_i^2,N]。
但块筛是前密后疏的，极其不均匀的，这启发我们把要转移的区间拆为：
[p_i,B]\cup[B,N/p_i]
它会转移到：
[p_i^2,p_iB]\cup[p_iB,N]
对前半部分的项做主动转移：
要考虑 (B-p_i) 个元素的转移。
对后半部分的项做被动转移：
相当于要考虑最后 \left(\frac{N}{p_iB}\right) 个块筛会被更新。
例如我们考虑倒数第 x 块，会被这样更新（被动转移）：
{\rm sum}\left[\left\lfloor\frac{N}{x+1}\right\rfloor+1,\left\lfloor\frac{N}{x}\right\rfloor\right] \leftarrow {\rm sum}\left[\left\lfloor\frac{N}{(x+1)p_i}\right\rfloor+1,\left\lfloor\frac{N}{xp_i}\right\rfloor\right]
因为每次转移时会被 p_i 个块贡献，所以被动转移要用树状数组维护。
这两部分的时间复杂度和为：
O\left(B-p_i+\frac{N}{p_iB}\right)\log N
平衡两部分：
当 p_i\ge N^{1/3} 时，取 B=p_i
时间复杂度为：
\log N\sum_{i=\pi(N^{1/3})}^{\pi(N)}\frac{N}{p_i^2}=O\left(N^{2/3}\right)
当 p_i\le N^{1/3} 时，取 B=\sqrt{\frac{N}{p_i}}
时间复杂度为：
\log N\sum_{i=1}^{\pi(N^{1/3})}\left(2\sqrt{\frac{N}{p_i}}-p_i\right)=O\left(N^{2/3}\right)
因此总时间复杂度为 O\left(N^{2/3}\right)。
推广
更一般的，我们考虑推广：
\tilde{F}(x)=\sum_{n\ge1}\frac{f(n)}{n^x}=\prod_{i\ge1}\left(1+A_i(x)\right)
其中 A_i(x) 是形如 A_i(x)=\sum_n a_{i,n}n^{-x} 的 DGF。
我们定义 t_i=\min\{k:a_{i,k}\neq0\}
不失一般性，令 t_i<t_{i+1}，
我们仍旧可以把 \tilde{F}(x) 拆成两类：
\begin{aligned}
\tilde{F}(x)&=\tilde{A}(x)\times\tilde{B}(x)\\
&=\prod_{i=1}^{M}(1+A_i(x))\times\prod_{i>M}(1+A_i(x))
\end{aligned}
其中 M=\max\{i:t_i<\sqrt N\}。
于是 \tilde{B}(x) 可以改写为求和：
\tilde{B}(x)=1+\sum_{i>M}A_i(x)
因此一个条件是满足此时 \tilde{B}(x) 的块筛好求。
一般来说这是容易满足的，可能需要用到质数前缀块筛。
这是一种 \widetilde O(N^{2/3}) 的质数前缀块筛。
我们考虑反向遍历 i 为 M\to1。
对于 t_i 的转移可以看作 [t_i,N/t_i]\xrightarrow{\times t_i}[t_i^2,N]。
如果 t_i 较为均匀，按照先前的根号分治，可以做到 O(N^{2/3}\log N)(M/\sqrt N)。
对于常见的 M=\pi(\sqrt N) 即可直接得出这里的时间复杂度为 O(N^{2/3})。
由于 A_i(x) 可能有好几项，设另一项 v_i>t_i 且满足 a_{i,v_i}\neq0，
对于 v_i 的转移可以看作 [t_i,N/v_i]\xrightarrow{\times v_i}[t_iv_i,N]。
我们可以这么估计时间复杂度：
f_i(x)=x-t_i+\frac{N}{v_ix}\qquad(x>t_i)
\log N\sum_{i=1}^Mf_i(x)_{min}
根据这个方法，仍然假设 t_i 较为均匀，我们考虑两个特殊值：
若 v_i=kt_i+o(t_i)，可以做到 O((N/k)^{2/3}\log N)(M/\sqrt N)。
若 v_i=kt_i^2+o(t_i^2)，可以做到 O((N/k)^{1/2}\log N)(M/\sqrt N)。
因此对于大部分情况都可以做到 O(N^{2/3})，具体情况需要具体分析。
展望与对比
按照推广后的方法，对于先前的例题直接筛 \tilde{\Psi} 也是一样的时间复杂度。
不过这要多求一次质数前缀块筛。
若一个积性函数有重数函数，那么其一定可以写成上述 DGF 形式，
那么对于大部分的数论函数的重数函数前缀和问题都可以做到 O(N^{2/3})。
例如 \sum_x[\sigma_k(x)\le N] 等。
若该积性函数仅在 PN 处有值，那么它的重数函数前缀和问题可以做到 \widetilde O(N^{1/2})。
其实这个算法还有不少可以优化的地方，但是很麻烦：
例如主动转移时有不少值是 0，不用转移，于是记录一张指向后面跳的表。
我们再重新平衡一下 B，时间复杂度应该能快一个 \log^{1/3}N。
我和 _fewq 讨论后发现可以用 ln-exp 优化：
\begin{aligned}
\tilde{A}(x)&=\prod_{i=1}^M\left(1+A_i(x)\right)\\
&=\exp\left(\sum_{i=1}^M\ln(1+A_i(x))\right)\\
&=\exp\left(\sum_{i=1}^M\sum_{j\ge1}\frac{(-1)^{j-1}}{j}A_i(x)^j\right)
\end{aligned}
于是时间复杂度是 T(N)+\sqrt Npoly(n) 的。
其中 T(n) 是求 块筛 \exp 的时间复杂度，n 是 a_{i,n} 非零项个数。
我们可以使用 O(\log N) 次的 O(N^{2/3}\log^{1/3}N) 杜教筛来做 块筛 \exp。
周康阳的集训队论文把块筛卷积做到了 O(\sqrt N\log^{1.5}N\log\log N) 的时间复杂度，
但是常数巨大，在 OI 范围没用（跑的更慢）。
这算是另一种方法了。
对比这两种方法，\exp 方法实现可能还更简洁，但时间复杂度大了 4/3 个 \log。
希望哪位大佬可以实现对比一下，不知道能不能抵上这个算法的常数。
但正如前文的展望，我的算法能在某些要求上做到更快的 \widetilde O(N^{1/2})。
这也算一些优势罢，希望我的方法还是有用武之地的。
后话
我最开始想到的方法是根号分治套根号分治的，还有巨量的分类讨论。
但我再不断提取核心思想，改到了现在这个版本，算法逻辑已经很简单了。
关于块筛还有很多科技：
我 和 飞雨烟雁 的这两篇文章探讨了块筛的牛顿迭代，各位可以看看。
周康阳的集训队论文也探讨了一种理论时间复杂度很低的块筛卷积，不过常数也很大。