万能欧几里得算法

发现自己在算法方面还有很多欠缺的，趁着还有时间赶紧补一下。

Intro

万能欧几里得算法解决的是出现 \(\lfloor\frac{ai+b}{c}\rfloor\) 的求和式，其中 \(i\) 是求和指标。

几何意义转化一下，发现 \(\lfloor\frac{ai+b}{c}\rfloor\) 表示的是 \(y=\frac{ax+b}{c}\) 这条直线在 \(x=i\) 处的下方的整点数。

我们只关心直线下方的整点数，可以把这条直线画在网格图坐标系上，从左到右考虑这条直线每次穿过网格边的情况。维护一个变量 \(t\) 表示当前横坐标下面有多少个整点。

• 首先，在 \(x=0\) 的地方，要有 \(t\leftarrow \lfloor\frac{b}{c}\rfloor\)。
• 当这条直线向上穿过一条横着的网格边，有 \(t\leftarrow t+1\)。
• 当这条直线向右穿过一条横着的网格边，让 \(t\) 按照某种方式贡献进答案中，贡献方式取决于要求的具体式子。
• 当这条直线向右上穿过一个整点，让 \(t\leftarrow t+1\) 然后再贡献。

如果把每次 \(t\leftarrow t+1\) 记为 \(U\) 操作，贡献答案记为 \(R\) 操作，那么这条直线在定义域 \([0,n]\) 的与答案贡献有关的所有信息可以由一个由 \(U,R\) 组成的序列来表示。而且一般可以用矩阵乘法表示。

例如，假如我们要求这个式子：

\[\sum_{i=0}^n\lfloor\frac{ai+b}{c}\rfloor \]

那么我们可以用一个行向量 \((\sum t,t,1)\) 来表示当前的状态。 \(U\) 操作相当于右乘上矩阵 \(\begin{bmatrix}1,\ 0,\ 0\\ 0,\ 1,\ 0\\ 0,\ 1,\ 1\end{bmatrix}\)，\(R\) 操作相当于 \(\begin{bmatrix}1,\ 0,\ 0\\ 1,\ 1,\ 0\\ 0,\ 0,\ 1\end{bmatrix}\)。最后答案就是 \(\sum t\)。

于是这个东西总算有了点用处，可以让我们在 \(O(n+\frac{an+b}{c})\) 的时间内求出这个式子的值，但这不是比暴力还慢吗。

考虑优化。把问题描述为矩阵乘法最大的好处就是，这个运算有结合律，可以快速幂计算。设 \(f(a,b,c,n,U,R)\) 为 \(y=\frac{ax+b}{c}\) 在 \([0,n]\) 上，每向上穿过一次横线就乘上一个 \(U\)，向右一次就乘上一个 \(R\) 最终得到的矩阵。

根据式子可以知道，直线为 \(y=\frac{ax+b}{c}\) 与对于任意的 \(i\)，第 \(i\) 个 \(R\) 前面有恰好 \(\lfloor\frac{ai+b}{c}\rfloor\) 个 \(U\) 是等价的。

如果 \(b\geq c\)，那么我们可以先把开头的 \(\lfloor\frac{b}{c}\rfloor\) 个 \(U\) 算了，有：

\[f(a,b,c,n,U,R)=U^{\lfloor\frac{b}{c}\rfloor}f(a,b\bmod c,c,n,U,R) \]

如果 \(a\geq c\)，那么中间的每个 \(R\) 之前必定会出现至少 \(\lfloor\frac{a}{c}\rfloor\) 个 \(U\)，考虑把这个东西合并起来。合并之后第 \(i\) 个 \(R\) 前面 \(U\) 的个数为：

\[\lfloor\frac{ai+b}{c}\rfloor-i\lfloor\frac{a}{c}\rfloor=\lfloor\frac{(a\bmod c)i+b}{c}\rfloor \]

所以：

\[f(a,b,c,n,U,R)=f(a\bmod c,b,c,n,U,U^{\lfloor\frac{a}{c}\rfloor}R) \]

如果 \(a<c\)，我们仿照类欧几里得算法，对称坐标系，相当于交换 \(U,R\)。新的线段的定义域 \(n'\) 变成原来的值域 \(\lfloor\frac{an+b}{c}\rfloor\)，参数要满足第 \(i\) 个 \(R\) 之前 \(U\) 的个数等于之前线段第 \(i\) 个 \(U\) 之前 \(R\) 的个数，而后者等于最大的 \(j\) 满足：

\[\begin{aligned} \lfloor\frac{aj+b}{c}\rfloor&<i\\ \frac{aj+b}{c}&<i\\ j&<\frac{ci-b}{a}\\ j&<\lceil\frac{ci-b}{a}\rceil\\ j&\leq\lfloor\frac{ci-b-1}{a}\rfloor \end{aligned} \]

但是这样直接做会出现一些问题。

首先截距 \(\frac{-b-1}{a}\) 是负数，注意如果我们把变化后的线段往左平移一格之后截距就是正数了，因为原线段截距 \(<1\)，原线段向下平移一格之后与 \(x\) 轴的交点肯定来到正半轴。这相当于我们把前面的 \(R^{\lfloor\frac{c-b-1}{a}\rfloor}U\) 的提前拿出来然后再翻转。注意能平移一格的条件是原线段至少有一个 \(U\)，即 \(n'>0\)。如果 \(n'=0\)，那么直接返回 \(R^n\)。

如果成功平移了一格，那么变换后的线段第 \(i\) 个 \(R\) 之前的 \(U\) 的个数变为：

\[\lfloor\frac{c(i+1)-b-1}{a}\rfloor-\lfloor\frac{c-b-1}{a}\rfloor=\lfloor\frac{ci+(c-b-1)\bmod a}{a}\rfloor \]

其次变换之后的线段后面会多出来一些 \(U\)，这就是让我们把原线段最后面 \(n-\lfloor\frac{cn'-b-1}{a}\rfloor\) 个 \(R\) 提前拿出来算再翻转。

综上所述，有：

\[f(a,b,c,n,U,R)=R^{\lfloor\frac{c-b-1}{a}\rfloor}Uf(c,(c-b-1)\bmod a,a,n'-1,R,U)R^{n-\lfloor\frac{cn'-b-1}{a}\rfloor} \]

于是这样就可以 \(O(\log\max\{a,c\})\) 计算一次了。非常厉害。关键在于只要你能设计出一个用 \(U,R\) 表示答案的矩阵，就能直接套这个模板，非常厉害。

现在让我们来做一做模板题 P5170。如果你真的拿矩阵维护这一坨东西，你就会发现根本过不去。注意只要信息能表示成一个幺半群，那么我们就能快速计算，而矩阵慢的原因是其维护了太多无用信息。

先来看这个式子：

\[\sum_{i=0}^n\lfloor\frac{ai+b}{c}\rfloor \]

首先我们的模板不能计算 \(i=0\) 的贡献，所以最后要加上 \(\lfloor\frac{b}{c}\rfloor\)。

设记号：

\[\sum_S f(x,y)=\sum_{i}[S_i=R]f(\sum_{j\leq i}[S_j=R],\sum_{j\leq i}[S_j=U]) \]

其中 \(S\) 是一个由 \(U,R\) 构成的操作序列，\(f(x,y)\) 是一个关于 \(x,y\) 的函数。

那么设这道题最终求出来操作序列为 \(S\)，我们要求的就是 \(\sum_Sy\)。

设一段操作序列 \(S\) 的信息为 \((x,y,s)\)，其中 \(x\) 表示 \(\sum_{j\leq i}[S_j=R]\)，\(y\) 表示 \(\sum_{j\leq i}[S_j=U]\)。

考虑合并两个操作序列 \(S_1,S_2\)，信息会如何合并：

\[(x_1,y_1,\sum_{S_1}y)\times (x_2,y_2,\sum_{S_2}y)=(x_1+x_2,y_1+y_2,?) \]

推导最后一项如何合并：

\[\sum_{S_1+S_2}y=\sum_{S_1}y+\sum_{S_2}(y+y_1)=\sum_{S_1}y+\sum_{S_2}y+x_2y_1 \]

于是合并结果为：

\[(x_1+x_2,y_1+y_2,\sum_{S_1}y+\sum_{S_2}y+x_2y_1) \]

显然 \(U=(0,1,0),R=(1,0,0)\)，而且有幺元 \((0,0,0)\)。并且满足结合律。所以可以直接做了。

来推剩下的几个式子。第二个式子要求的就是 \(\sum_S y^2\)。直接推：

\[\begin{aligned} &\quad\sum_{S_1+S_2}y^2\\ &=\sum_{S_1}y^2+\sum_{S_2}(y_1+y)^2\\ &=\sum_{S_1}y^2+\sum_{S_2}y_1^2+2y_1y+y^2\\ &=\sum_{S_1}y^2+x_2y_1^2+2y_1\sum_{S_2}y+\sum_{S_2}y^2 \end{aligned} \]

第三个式子就是：

\[\begin{aligned} &\quad\sum_{S_1+S_2}xy\\ &=\sum_{S_1}xy+\sum_{S_2}(x+x_1)(y+y_1)\\ &=\sum_{S_1}xy+\sum_{S_2}xy+y_1\sum_{S_2}x+x_1\sum_{S_2}y+x_1y_1\sum_{S_2}1\\ &=\sum_{S_1}xy+\sum_{S_2}xy+y_1\frac{(1+x_2)x_2}{2}+x_1\sum_{S_2}y+x_1y_1x_2 \end{aligned} \]

至此我们只需要每个信息维护 \((x,y,\sum_Sy,\sum_Sxy,\sum_Sy^2)\) 即可合并。

其他的一些例题

上面说过，万能欧几里得的好处在于，你可以直接将工作投入到推贡献式子上，而不用每次重新推导一坨互相转移的式子。

LOJ 138

模板题。读者自推用于检验学习成果。

提示：二项式定理展开。

LOJ 6440

模板题。读者自推用于检验学习成果。

提示：把矩阵这个整体当作半群里的一个信息来维护。

P5172

首先 \(\sqrt r\) 为整数的情况是平凡的，可以直接特判处理。现在 \(\sqrt r\) 为无理数。

此时直接应用万能欧几里得算法：

信息可以设为 \((x,s)\)。\(x\) 表示 \(-1\) 的当前直线下方点的个数次幂，\(s\) 则表示答案。此时有：

\[U= \begin{bmatrix} -1&, 0\\ 0&, 1 \end{bmatrix}， R= \begin{bmatrix} 1,\ 0\\ 1,\ 1 \end{bmatrix} \]

那么考虑推到万欧函数。注意这道题 \(\sqrt r\) 是无理数，所以不能直接沿用之前的算法。设 \(f(n,k,U,R)\) 表示正在处理 \(y=kx\) 这条直线，定义域 \([1,n]\) 的答案。还是类似对称坐标系的想法：

• 若 \(k\geq 1\)，则每个 \(R\) 前面至少出现 \(\lfloor k \rfloor\) 个 \(U\)：

\[f(n,k,U,R)=f(n,k-\lfloor k \rfloor,U,U^{\lfloor k \rfloor}R) \]

• 若 \(k<1\)，则对称坐标系，与上面同样需要注意的点在于，末尾的 \(R\) 的特殊处理：

\[m=\lfloor nk\rfloor\\ f(n,k,U,R)=f(m,\frac{1}{k},R,U)R^{n-\lfloor\frac{m}{k}\rfloor}\]

分析一下复杂度。

• 对于 \(k\geq 1\)，会花费 \(1\) 步变成 \(k<1\)。
• 对于 \(k<1\)，若 \(k<\frac{1}{2}\)，则 \(n\) 折半。否则，下一轮 \(k<1\) 时 \(k\) 会变为 \(\frac{1}{k}-1\)，此时 \(k(\frac{1}{k}-1)=1-k<\frac{1}{2}\)，\(n\) 会折半。

所以复杂度是 \(O(\log n)\) 的。

注意如果直接用 long double 实现，精度会爆掉。有两种解决方法：

• 使用 __float128。
• 注意到 \(k\) 的迭代过程就是求连分数展开，可以直接求连分数展开。