UOJ Logo EntropyIncreaser的博客

博客

标签
暂无

任意模数二项卷积

2020-07-20 08:48:26 By EntropyIncreaser

二项卷积:

$$ c_n = \sum_k \binom n k a_k b_{n-k} $$

当模数 $M$ 有一质因子 $p\le n$ 时,计算二项卷积无法直接变为 EGF 进行计算,因为 $n!$ 此时不可逆。


我们考虑首先如何解决 $M = p^k$ 时的情况,然后可以用 CRT 合并各情况。

我们记 $v_p(n)$ 是 $n!$ 中 $p$ 的质因子次数,$p$-阶乘为 $n!_p = p^{v_p(n)}$,反 $p$-阶乘为 $\overline{n!_p} = \frac{n!}{n!_p}$。

那么由定义显然 $\overline{n!_p}$ 还是同余 $M$ 可逆的。我们先令 $\widehat a_n = a_n \cdot \left( \overline{n!_p} \right)^{-1} \bmod M$,我们可以得到

$$ \widehat c_n \equiv \sum_k \left(\frac{n!_p}{k!_p (n-k)!_p}\right) \widehat a_k \widehat b_{n-k} \equiv \sum_k p^{v_p(n)-v_p(k)-v_p(n-k)} \widehat a_k \widehat b_{n-k} \pmod M $$

而 $d = v_p(n)-v_p(k)-v_p(n-k)$ 可以推出就是在 $p$ 进制下,$n$ 减去 $k$ 时所发生的退位次数。(这个被叫做 Kummer 定理,过程略)

而 $n$ 在 $p$ 进制下最多只有 $\log_p n$ 位可退,因此我们知道 $p^d \le n$,因此我们在不取模的情况下,可以得到 $\widehat c_n \le n \cdot nM^2 = n^2M^2$。

虽然 $p$ 在模 $M$ 下不可逆,但是当 $p\le n$,自然满足在我们选取的 NTT 模数下都可逆!因此,这一涉及除法的卷积式子,因为已经保证了结果是值域在 $n^2M^2$ 内的整数,所以我们只需选取 NTT 模数进行卷积,之后用 CRT 合并即可。取 $n\le 10^6, M\le 10^9$ 的一般情况下,可得 $c_n \le 10^{30}$,使用四个 NTT 模数进行合并足够。

因此对于每个 $M=p^k$ 的情况,我们都可以通过“四模数 NTT”进行计算,那么对于一般的 $M$ 进行 CRT,算法的复杂度为 $\Theta(n\log n \cdot \omega(M))$,或者也可以解释为,结果值域是 $n^{1+\omega(M)}M^2$ 的卷积。


该提交记录 是一个不算精细的实现。在 $\omega(M)$ 取到最大时大概效率是 loj 的 1s 上下波动(其实这个时候每次 CRT 不用四模数了,如果考虑这一点还能更快)。暂时没想到在“四模数 NTT”的时候避免 int128 的方法。


特殊情况:

  • $p$ 很小的时候可以类比子集卷积,加一个占位多项式来帮助计算,但是这个做法的复杂度看起来会带 $p$ 和 $\log_p n$,在 $p$ 大的时候尚未想到能从此想法得到什么优秀的做法。之前思考这个问题的时候也在这里困惑了很久。

HALF-GCD 算法的阐述

2020-03-10 17:48:53 By EntropyIncreaser

翻了各课件,发现……

某个英文网站 是这么说的:

There is a way to calculate the GCD and resultants in $O(n\log^2n)$. To do this you should note that if you consider $a(x)=a_0+x^ka_1$ and $b(x)=b_0+x^kb_1$ where $k=\min(\deg a,\deg b)/2$ then basically the first few operations of Euclidean algorithm on $a(x)$ and $b(x)$ are defined by the Euclidean algorithm on $a_1(x)$ and $b_1(x)$ for which you may also calculate GCD recursively and then somehow memorize linear transforms you made with them and apply it to $a(x)$ and $b(x)$ to lower the degrees of polynomials. Implementation of this algorithm seems pretty tedious and technical thus it's not considered in this article yet.

你如果脑子一冲,傻乎乎地去直接分治一半算出线性变换,DEBUG 一下会让你冷静过来…… 你发现这个线性变换应用到整体多项式上没有任何效果……

而 Picks 的 Introduction to Polynomials 也写的比较晦涩,因此在这里稍微详细地谈谈,以正视听……


给多项式 $a, b \in \mathbb F[z]$,我们原本的问题是求 GCD 以及扩展欧几里得的解:$x, y\in \mathbb F[z], ax + by = \gcd (a, b)$

由于多项式取模最坏情况下每次只让一个多项式的度数减小,这种方法的复杂度最好也就是 $\Theta((\deg a) (\deg b))$ 的。

我们考虑求解一个中间情况叫做 HALF-GCD,设 $k = \min(\deg a, \deg b)$,这个东西让我们对于 $a, b$ 能得到一个矩阵 $\begin{pmatrix} x_1 & y_1 \\ x_2 & y_2 \end{pmatrix}$,其中 $\deg x \le \frac k 2 + o(k)$,且 $\deg (a x_1 + b y_1), \deg (a x_2 + b y_2) \le \frac k 2 + o(k)$(除非它们的 $\gcd$ 的 $\deg$ 比这个大)。也就是说这个 HALF-GCD 可以理解成是将欧几里得过程执行到一半的产物。

睿智的你想必已经意识到了这个东西如何得出真正的欧几里得的解。我们只需要调用一次这个算法,然后就规约到了一个 $\frac n2$ 的子问题(可能需要额外支付 $O(M(n))$ 时间用于调整较大的多项式),那么我们就得到了一个 $\Theta(M(n)) + T(n) + T(n/2) + T(n/4) + \cdots$ 的方法。

那么接下来考虑通过分治来解决 HALF-GCD:设 $k=\min(\deg a, \deg b)$,我们首先算出 $a/z^{k/2}, b/z^{k/2}$ 对应的矩阵,这样一个矩阵作用与原多项式有两部分:考虑将原来的表为 $a=a_0 + a_1 z^{k/2}, b=b_0 + b_1 z^{k/2}$,那么 $ax + by$ 中,$\deg ((a_1 x + b_1 y)z^{k/2}) = \frac 34 k + o(k)$,另外 $\deg (a_0 x + b_0 y) = \frac 34 k + o(k)$ 也比较显然。接着我们将得到的两个 $\frac 34 k$ 次多项式除掉前 $z^{k/4}$ 项送进去,这样又再次保证了递归返回的又是一个 $\frac 14 k + o(k)$ 量级的结果,我们将两个矩阵相乘就得到了变换矩阵。而结果正好满足度数在 $\frac k 2 + o(k)$ 以内。时间复杂度 $T(n) = 2T(n/2) + M(n)$,即 $M(n)\log n$。消耗了 $M(n)$ 的一方面是进行矩阵乘法,另一方面是可能在一些情况减少两个多项式次数的差,同时也可能发生剪枝。

综上所述,我们在 $M(n)\log n$ 的时间内完成了两个多项式的欧几里得求算,高精度整数的做法也是相似的。

这里 是一份代码,这里处理的时候是令 $k = \max (\deg a, \deg b)$,理论上是差不多的,目前尚且不知道这个代码是否存在一些 edge case 上的问题,欢迎大家指正。

区间增量最大子段和的 polylog 做法

2019-09-19 22:07:18 By EntropyIncreaser

先抛出目前得到的结果:区间给每个数加 $x (x \ge 0)$,区间查询最大子段和,我们有一个 $\mathcal O(n\log^3 n + m\log^4 n + q\log n)$ 的做法(不排除可能通过一些分析发现这个算法的复杂度实际上更低)。

其中,$n$ 是序列长度,$m$ 是修改次数,$q$ 是询问次数。

我们考虑当前情况下先用传统的方法用线段树维护最大子段和,即维护 $sum, lmax, rmax, totmax$,也就是:

$$ \begin{align} sum &= ls.sum + rs.sum\\ lmax &= \max(ls.lmax, ls.sum+rs.lmax)\\ rmax &= \max(rs.rmax, rs.sum+ls.rmax)\\ totmax &= \max(ls.totmax, rs.totmax, ls.rmax + rs.lmax) \end{align} $$

但是我们同时还维护一个额外的信息,就是当前这个值在本区间 $+x$ 后的情况,也就是说我们存储的是一个一次函数 $kx + b$,我们将每个量都看成一个一次函数

在这一情况下,当所加的 $x$ 极小,我们每个 max 符号选取的位置不会发生变化,那么我们还想记录最小的 $x$,表示将这个子树内的数都 $+x$ 后,某个 max 改为选取另一条直线,我们称这是一次“击败”事件。在这种情况下我们向下递归到该节点并重新选取当前 $x$ 下值最大的直线。

接下来证明本算法的时间复杂度。

我们记一个包含 max 比较的变量的 rank 值 $r(v, para)$:

  • 对于 $r(v, lmax)$ 或 $r(v, rmax)$:对于该变量的决定式 $f(x) = \max(a(x), b(x))$,rank 值表示当前 $a, b$ 中斜率大于 $f$ 的斜率的直线的数量乘以 $d^2$,$d$ 为该节点到根节点距离的节点数。
  • 对于 $r(v, totmax)$:对于该变量的决定式 $f(x) = \max(a(x), b(x), c(x))$,rank 值表示当前 $a, b, c$ 中斜率大于 $f$ 的斜率的直线的数量乘以 $d$。

我们直接定义势能 $\Phi = \sum_v \sum_{para} r(v, para)$。

我们现在考虑最坏实现,那就是每次击败都重新从根节点向下找到该节点修改。分析一次击败事件的摊还代价,我们记单位时间为一次线段树修改,因此这里算出的代价换算到时间复杂度后要乘以一个 log:

  • 发生一次 $lmax$ 或 $rmax$ 的击败(此处以 $lmax$ 举例):
    • 当前节点的 $\Delta r(v, lmax)$ 显然是 $-d^2$,因为替换上来的函数显然斜率更大。
    • 父节点的 $\Delta r(prt, lmax) \le (d-1)^2$,因为替换上来的函数有可能使得在上方的节点增加一次被击败的可能性。
    • 父节点的 $\Delta r(prt, totmax) \le d - 1$,同理。
    • 其他的 rank 值不应该有改变,注意不要混淆一个概念:如果此时上面有的函数因为下面的更新以及引发了修改,那应当视作是另一次击败事件。
    • 综上,我们计算这次击败的势能变化:

$$ \Delta \Phi \le 1 -d^2 + (d-1)^2 + d-1 = 1-d \le 0 $$

  • 发生一次 $totmax$ 的击败:
    • 当前节点的 $\Delta r(v, totmax)$ 是 $-d$。
    • 父节点的 $\Delta r(prt, totmax) \le d - 1$。
    • 其他的 rank 值不应该有改变。
    • 综上,我们计算这次击败的势能变化:

$$ \Delta \Phi \le 1 -d + d-1 = 0 $$

可以看出,我们完全消解了击败操作的摊还代价!这意味着所有的时间复杂度由初始势能以及修改造成的额外势能改变贡献。

$$ \Phi_0 = \sum_v \sum_{para} r(v, para) \le \sum_v d + 2d^2 = \Theta(n\log^2 n) $$

显然一个节点 $v$ 的 $\sum_{para} r(v, para) \le d + 2d^2$,而一次区间修改会使得 $\mathcal O(\log n)$ 个节点的 rank 变大,因此 $\Delta \Phi = \mathcal O(\log^3 n)$。

别忘了这个计量的是线段树上的操作次数……因此实际复杂度是 $\mathcal O(n\log^3n + m\log^4n)$,询问显然是 $\Theta(q\log n)$ 没话说。


UPD

我们考虑每个节点的 $lmax$ 和 $rmax$ 的决策发生变化的次数。注意到 $lmax$ 的决策点只会递增,而在一个节点发生斜率切换的必要条件是决策位置从 $[l, mid]$ 切换到了 $[mid+1, r]$。每个节点最多触发一次,故 $lmax, rmax$ 的总共额外复杂度是 $\Theta(n\log n)$。

因此我们现在知道 KTT 维护 $lmax$ 和 $rmax$ 的总共 dfs 到的节点数是 $\Theta((n + q)\log n)$ 的。我们沿用势能分析,是 $\Theta((n+q)\log^3 n)$。

共 3 篇博客