狄利克雷前缀和

Aleph_Drawer · 2024-02-26 22:17:34 · 算法·理论

0 前言

膜拜 GhostLX 佬。膜拜 zak！

耗时三天。效率极为低下。

1 狄利克雷前缀和

狄利克雷（Dirichlet）前缀和通常指的是如下问题：

给定 \{a_i\}, 1 \leq i \leq n，求解 \{b_i\}, 1 \leq i \leq n，其中有
b_i = \sum_{d \mid i} a_d
越快越好。

首先我们知道这个问题可以简单的 \mathcal O(n \log n) 解决，方法不再赘述。不妨假设 n = 10^7，此时这个时间复杂度就相当勉强了。而现在介绍的就是一种能在 \mathcal O(n \log \log n) 解决这种问题的方法。

1.1 从狄利克雷卷积分析本质

还记得狄利克雷卷积吗？可以回去看看。不妨令 A(x) = a_x，B(x) = b_x，那么此时 B = A * \mathbf 1。

我们已经知道，两个积性函数卷起来还是积性函数。不妨针对 \mathbf 1 下手。

不妨令 f_i(p^k) = [p = p_i]，其中 p_i 是第 i 小的质数。特别的，我们规定 f(1) = 1，且 f 为积性函数，这就意味着，如果 x 不是特定质数的特定次幂，那么 f 的值恒定为 0。那么 \mathbf 1 = \prod_i f_i。这里卷积显然是狄利克雷卷积。

怎么说明这一点呢？首先我们知道唯一分解定理，知道每个数都会被唯一分解成若干质因数的幂的乘积。

不妨常规地设 x = \prod_i p_i^{c_i}，不妨从小到大的去把 f 给卷起来，观察到卷到 x 最大的质因子的时候，得到的函数在 x 出取值恒定为 1。这个可以用归纳法来说明。对于 f_1，根据定义很容易说明所有形如 2^x 形式的式子值均为 1。

不妨假设形如 2^{c_1} \cdot 3^{c_2} \cdot \dots p_{i - 1}^{c_{i - 1}} 类的值均在对应位置取到了 1，此时准备把这个卷起来的庞然大物（不妨设为 F(x)）和 f_i 卷。对于 2^{c_1} \cdot 3^{c_2} \cdot \dots p_i^{c_i} 之类的数，我们观察到 F(2^{c_1} \cdot 3^{c_2} \cdot \dots p_{i - 1}^{c_{i - 1}}) = 1，而且 f_i(p_i^{c_i}) = 1，又根据唯一分解定理，所以有 (F * f_i)(^{c_1} \cdot 3^{c_2} \cdot \dots p_i^{c_i}) = 1。那么我们就证明了 \mathbf 1 = \prod_i f_i，证毕。

那么我们原来不是要求 A * \mathbf 1 吗？那么现在就变成了 A * f_1 * f_2 * \dots。我们怎么去做 A * f_i 呢？

首先观察 x \nmid p_i 的值，观察到就是原来的值，所以我们不需要去动他们。那么剩下的就是 x \mid p_i 的值了。

观察到对于 x \mid p_i，我们都要加上 A_{x / p_i} 的值。

对于 x \mid p_i^2，我们都要加上 A_{x / p_i^2} 的值。

以此类推，我们发现我们可以利用后效性去维护这个东西。具体的，以 p = 2, x = 12 为例，对于 A(12) 我们要加上 A(3), A(6)，但是如果我们先做了 A(6)，已经加上了 A(3)，那么在处理到 A(12) 的时候我们只需要加上 A(6) 即可。

所以最后代码并不长。

for(int i = 1; p[i] <= n; i++)
    for(int d = 1; d * p[i] <= n; d++)
        a[p[i] * d] += a[d];

分析一下时间复杂度，发现是 \mathcal O(\sum_{p} \dfrac np)。这个是埃筛的经典时间复杂度，就是我们要的 \mathcal O(n \log \log n)。

可能有人会问，这个时间复杂度是怎么出来的。

其实可以简单研究一下。

References：https://www.zhihu.com/question/35112789

不严谨的，我们知道 \pi(x) \sim \dfrac x{\ln x}，那么对于函数 a(x) = [x \in P]，我们有 a(x) = \pi(x) - \pi(x - 1) \sim \dfrac{1}{\ln x}。

则

\begin{aligned} \sum_{p \leq n} \dfrac 1p &= \sum_{i = 2}^{n} \dfrac 1i a(i) \\ &= \sum_{i = 2}^n \dfrac 1{i \ln i} \\ &\sim \int_2^n \dfrac 1{x \ln x} \mathrm dx \end{aligned}

经常做微积分的朋友都知道 f(x) = \dfrac 1{x \ln x} 的反微分是 F(x) = \ln \ln x + C（注意限定 x > 1），那么上面那个式子就是 F(n) - F(2) 就是 \ln \ln n - \ln \ln 2。所以这个式子就近似于 \ln \ln n。那么时间复杂度就是 \mathcal O(n \ln \ln n)。

2 一些拓展

2.1 从特例到一般

上面展示了 A * \mathbf 1 的快速方法，那么对于数论函数 f 和积性函数 g，是否也可以类似的去 \mathcal O(n \ln \ln n) 的计算 A * g 呢？答案是可以的！

首先我们要 O(n) 的处理出 g 函数的值。这个可以用线性筛做到。

接下来模仿上面的思路，构建 h_{i}(p^x) = \begin{cases} g(p), & p = p_i \lor p^x = 1, \\ 0, &\text{otherwise} \end{cases}，然后我们可以用同样的方法说明 g = \prod_i h_i，最后我们利用类似的方法来做 f * h_i，具体的，对于 x \nmid p_i，我们不需要做任何事情，而对于 x \mid p_i，我们就顺序去做，令 A_{dp_i} \leftarrow A_{dp_i} + A_{d} \cdot g(p_i)。

但是这会有问题！诸如 n = 4 在 p = 2 进行到这一步时会算出 g^2(2) A_1 + g(2) A_2 + A_4，不太能忍受。

所以我们退而求其次，我们放弃后效性的利用，转而分别对 g(2), g(4) 等相关的值分批做操作。具体的，枚举质数 p，枚举 k 使得 p^k \leq n，再枚举 d 使得 dp^k \leq n。我们使用临时数组 tmp 来存储额外加上的值，具体的，tmp_{\normalsize dp^k} \leftarrow tmp_{\normalsize dp^k} + A_d \cdot g(p^k)（临时使用 \normalsize）。

然后最后把 tmp 转移到 a 里面。

有人可能问，这么做复杂度不会爆炸吗？分析一波。我们的时间消耗是 \sum_{p \leq n} \sum_{p^k \leq n} \dfrac n{p^k}。

我们尝试证明：

\mathcal O(\sum_{p^k \leq n} \dfrac n{p^k}) = \mathcal O(\dfrac np)

观察到后面那个东西可以等差数列求和，首先这个限制是毫无用处的，因为当 p^k > n 的时候后面这个值趋近于 0，我们可以简单的忽略。

令 S = \sum_{k} \dfrac n{p^k}，则 \dfrac {S}{p} = \sum_{k} \dfrac n{p^{k + 1}}，那么 \dfrac {S(p - 1)}p = \dfrac n p，S = \dfrac n{p - 1}，显然趋近于 \dfrac np，证毕。

所以最后时间复杂度依然很神奇的是 \mathcal O(n \log \log n)！

2.2 改前缀变后缀

修改问题变为如下形式：

给定 \{a_i\}, 1 \leq i \leq n，求解 \{b_i\}, 1 \leq i \leq n，其中有
b_i = \sum_{i \mid d} a_d
越快越好。

可以观察到此时只是从约数变为了倍数。这个怎么做呢？

仍然还是尝试借用上面的思路。

我们分开从小到大枚举 p，此时我们发现，把循环顺序 reverse 一下就是我们要的东西。

具体的，以 p = 2, A_2 为例，我们发现它会从 A_4 转移过来，而 A_4 又会从 A_8 转移来，实现了与上面类似的“后效性传递”。

原理其实和上面是一样的。

直接放代码吧。

for(int i = 1; p[i] <= n; i++)
    for(int d = n / p[i]; d >= 1; d--)
        a[d] += a[p[i] * d];

2.3 改前缀为差分

现在我们的问题又变了。

给定 \{b_i\}, 1 \leq i \leq n，求解 \{a_i\}, 1 \leq i \leq n，其中有
b_i = \sum_{d \mid i} a_d
越快越好。

注意到给定的东西与求解的东西互换了位置。

首先有一个朴素的思路，我们不是知道 B = A * \mathbf 1 吗？那么显然我们有 A = B * \mathbf 1^{-1} = B * \mu，先线性筛 \mu 然后直接按照 2.1 去做就行，但是常数比较大。

当然我们还有更好的思路，观察到我们用到的只是 \mu(p^k)，这个不是 1 就是 -1，而且非常有规律。

观察到 \mu(p^a) \times \mu(p^b) = \mu(p^{a + b})，这个是成立的，也就是说，在我们用到的范围内，\mu 是一个完全积性函数，那么首先我们就没有必要枚举 p^k 了，因为枚举的意义就是因为可能积性函数不是完全积性的，\mu(p^a) \times \mu(p^b) 就可能不等于 \mu(p^{a + b})。那么此时与一般的狄利克雷前缀和不同的，只剩下了：

式子要额外乘上 \mu(p)，即 -1；
这个式子不需要利用“后效性”。

这两个不同点都非常好解决，第一个改加法为减法，第二个就直接改变枚举顺序即可。

代码如下：

for(int i = 1; p[i] <= n; i++)
    for(int d = n / p[i]; d >= 1; d--)
        a[d] -= a[p[i] * d];

而对于狄利克雷后缀差分，就是二者的结合体，我们依然还是再次改变枚举顺序即可。

3 后记

就这样吧。好想睡觉。