万能欧几里得算法

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

万能欧几里得算法解决的是出现 \(\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\) 表示答案的矩阵，就能直接套这个模板，非常厉害。