万能欧几里得算法 笔记

\[ \]

类欧几里得算法

类欧几里得算法常用于解决形如 \(\sum_i \left\lfloor \frac{ai+b}c \right\rfloor\)（\(a, c > 0, b \ge 0\)）的问题。

代数推导，当 \(i\) 从 \(0\) 求和到 \(n\) 的答案为 \(f(n, a, b, c)\)：

  • \(a>c \lor b > c\)，变为 \(f(n, a \bmod c, b \bmod c, c)\)。

  • 考虑交换求和顺序，令 \(m=\left\lfloor \frac{an+b}c \right\rfloor\)：

    \[\sum_{i=0}^n\left\lfloor \frac{ai+b}c \right\rfloor = \sum_{j=0}^{m}\sum_{i=0} [j < \left\lfloor \frac{ai+b}c \right\rfloor] \]

    对于 \(\left\lfloor \frac{ai+b}c \right\rfloor\)，它等于 \(\left\lceil \frac{ai+b+1}{c} \right\rceil - 1\)。则有 \(j+1 < \left\lceil \frac{ai+b+1}c \right\rceil\)，变为 \(cj+c-b-1 < ai\)，可以得到 \(i > \left\lfloor \frac{cj+c-b-1}{a} \right\rfloor\)。因此 \(f(n, a, b, c)\) 可以变为 \(f(m, c, c-b-1, a)\)。

\(n\) 这一维是单调减的，\(a\) 这一维每两步至少除以二，因此时间复杂度是 \(O(\log n)\) 的。

它的几何意义：

  • 对第一步是把斜率大于等于 \(1\) 的都减去 \(\left\lfloor k \right\rfloor\)。
  • 对第二步则是把横纵坐标反转。

万能欧几里得算法

我们对类欧几里得算法进行了扩展，考虑一下它的几何意义。

考虑用折线来刻画一下这个求值。我们的直线每穿过一条横线的时候记为 \(U\)，穿过一条竖线记为 \(R\)。特别地，当它同时穿过横线和竖线的时候先 \(U\) 后 \(R\)。我们要求这个操作序列最后的位置为 \(R\)。

图例为 \(y=\frac{12.2x+2}{9}\)：

euclidean-graph1.png

此时 \(x\) 轴为 \(i\)，我们维护 \(\left\lfloor \frac{ai+b}{c} \right\rfloor\) 和 \(\sum\left\lfloor \frac{ai+b}{c} \right\rfloor\) 的值。当 \(x \to x+1\) 时，\(\left\lfloor \frac{ax+b}c \right\rfloor\) 会加 \(t\) 次（向上走 \(t\) 步），而 \(\sum\left\lfloor \frac{ai+b}{c} \right\rfloor\) 会加上 \(\left\lfloor \frac{ax+b}c \right\rfloor\)。我们使用矩阵维护（\(\begin{bmatrix}\left\lfloor \frac{ai+b}c \right\rfloor \\ \sum\left\lfloor \frac{ai+b}c \right\rfloor\end{bmatrix}\)），则每次向上的时候让将其乘上 \(U\)，向右乘上 \(R\)。令 \(f(n, a, b, c, U, R)\) 为上述问题的答案（即有 \(n\) 个 \(R\)，第 \(i\) 个 \(R\) 前有 \(\left\lfloor \frac{ai+b}c \right\rfloor\) 个 \(U\)，最后一个是 \(R\) 的答案），则：

  • 当 \(b \ge c\) 的时候，答案为 \(U^{\left\lfloor b/c \right\rfloor}f(n, a, b \bmod c, c, U, R)\)。

  • 当 \(a \ge c\) 的时候，答案为 \(f(n, a \bmod c, b, c, U, U^{\left\lfloor a/c \right\rfloor}R)\)。

    即我们变为 euclidean-graph2.png

  • 此时 \(a < c \land b < c\)，我们关于 \(y=x\) 对称，即我们对每个 \(U\) 统计其之前的 \(R\) 的个数（这样才能变为一个类似于欧几里得算法的形式）。（几何意义：对于整点的特殊情况，我们平移 \(\frac 1c\) 解决）

    euclidean-graph3.png

  • 对于第 \(i\) 个 \(U\)，求 \(\sum_{j \ge 1} \left[\left\lfloor \frac{aj+b}{c} \right\rfloor < i\right]\)。对 \(i > \left\lfloor \frac{aj+b}c \right\rfloor\)，有之前的类似转化，变为：

    \[i > \frac{aj+b}c \]

    \[j < \frac{ci-b}a \]

    \[j < \left\lceil \frac{ci-b}a \right\rceil \]

    \[j \le \left\lfloor \frac{ci-b-1}a \right\rfloor \]

    因此第 \(i\) 个 \(U\) 前会有 \(\left\lfloor \frac{ci-b-1}a \right\rfloor\) 个 \(R\)，令 \(m = \left\lfloor \frac{an+b}c \right\rfloor\)。我们可以把 \((n, a, b, c)\) 变为 \((m, c, -b-1, a)\)。但是 \(-b-1\) 是负数，我们初始要求 \(a, c \ge 0, b > 0\)。这一个可以直接把 \(i=1\) 的单独计算来解决。还有一个问题是我们要求最后一个是 \(R'\)（即 \(U\)），但是最后若干个变为了 \(U'\)（\(R\)）。因此我们也需要提出最后一个连续段。

    在提出 \(U^{\left\lfloor \frac{c-b-1}{a} \right\rfloor}\) 后，第 \(x\) 个 \(R\) 前应有 \(\left\lfloor \frac{cx+c-b-1}{a} \right\rfloor - \left\lfloor \frac{c-b-1}{a} \right\rfloor = \left\lfloor \frac {cx+(c-b-1)\bmod a}{a} \right\rfloor\)。最后应有 \(n-\left\lfloor \frac{cm - b - 1}{a} \right\rfloor\) 个 \(R\)。

    即 \(f(n, a, b, c, U, R) = R^{\left\lfloor (c-b-1)/a \right\rfloor}U\cdot f(m, c, (c-b-1)\bmod a, a, R, U) \cdot R^{n - \left\lfloor \frac{cm-b-1}{a} \right\rfloor}\)。

    特别地，由于我们需要掐头去尾，我们特判 \(m=0\) 的情况，此时答案为 \(R^n\)。

练习1：P5170 【模板】类欧几里德算法
Problem

多组询问，每次给定 \(a, b, c, n\)。

求 \(\sum_{i\le n} \left\lfloor \frac{ai+b}c \right\rfloor\)，\(\sum_{i\le n} i\left\lfloor \frac{ai+b}c \right\rfloor\)，\(\sum_{i\le n} \left\lfloor \frac{ai+b}c \right\rfloor^2\)。

Sol

令 \(y = \left\lfloor \frac{ax+b}c \right\rfloor\)。需要维护 \(x, y, \sum x, \sum y, \sum y^2, \sum xy\) 的值。

遇到 \(U\) 时，\(y \gets y + 1\)。

遇到 \(R\) 时，\(x \gets x+1\)，\(sx \gets sx + x\)，\(sy \gets sy + y\)，\(sy2 \gets sy2 + y^2\)，\(sxy \gets sxy + xy\)。

考虑合并两个区间的信息。则有 \(x = x_l + x_r\)，\(y = y_l + y_r\)，\(sx = sx_l + sx_r + x_l\cdot x_r\)，\(sy = sy_l + sy_r + y_l \cdot x_r\)，\(sy2 = sy2_l + sy2_r + x_ry_l^2 + 2y_l\cdot sy_r\)，\(sxy = sxy_l + sxy_r + x_rx_ly_l + x_lsy_r + y_lsx_r\)。

Code
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int P = 998244353;
struct Node {
	ll x, y, sx, sy, sy2, sxy;
	Node() : x(0), y(0), sx(0), sy(0), sy2(0), sxy(0) {}
	Node(ll _x, ll _y, ll _sx, ll _sy, ll _sy2, ll _sxy) : x(_x), y(_y), sx(_sx), sy(_sy), sy2(_sy2), sxy(_sxy) {}
	friend Node operator*(const Node &l, const Node &r) {
		Node res;
		res.x = (l.x + r.x) % P;
		res.y = (l.y + r.y) % P;
		res.sx = (l.sx + r.sx + l.x * r.x) % P;
		res.sy = (l.sy + r.sy + l.y * r.x) % P;
		res.sy2 = (l.sy2 + r.sy2 + r.x * l.y % P * l.y + 2 * l.y * r.sy) % P;
		res.sxy = (l.sxy + r.sxy + r.x * l.x % P * l.y % P + l.x * r.sy + l.y * r.sx) % P;
		return res;
	}
};
Node QPow(Node a, ll b) {
	Node res;
	for (; b; b >>= 1, a = a * a)
		if (b & 1)
			res = res * a;
	return res;
}
Node F(ll n, ll a, ll b, ll c, Node U, Node R) {
	if (!n) return Node();
	if (b >= c) return QPow(U, b / c) * F(n, a, b % c, c, U, R);
	if (a >= c) return F(n, a % c, b, c, U, QPow(U, a / c) * R);
	ll m = (a * n + b) / c;
	if (!m) return QPow(R, n);
	return QPow(R, (c - b - 1) / a) * U * F(m - 1, c, (c - b - 1) % a, a, R, U) * QPow(R, n - (c * m - b - 1) / a);
}
int main() {
	Node U, R;
	U.y = 1;
	R.x = 1, R.sx = 1;
	int T;
	scanf("%d", &T);
	while (T--) {
		ll n, a, b, c;
		scanf("%lld%lld%lld%lld", &n, &a, &b, &c);
		Node ret = F(n, a, b, c, U, R);
		ret.sy += b / c;
		ret.sy2 += (b / c) * (b / c);
		printf("%lld %lld %lld\n", ret.sy % P, ret.sy2 % P, ret.sxy);
	}
	return 0;
}
练习2：
Problem

求：

\[C=\sum\limits_{x=1}^LA^xB^{\left\lfloor\frac{Px+R}{Q}\right\rfloor} \]

其中 \(A,B\) 是 \(N\) 行 \(N\) 列的矩阵。\(N \le 20, L,\left\lfloor \frac{PL}{Q} \right\rfloor \le 10^{18}\)。

我们维护 \(\prod A, \prod B,\sum A^iB^j\)。

遇到 \(U\) 时：\(pb \gets pb \cdot B\)。

遇到 \(R\) 时：\(pa \gets pa \cdot A\)，\(sab \gets sab + pa\cdot pb\)。

考虑合并左右区间的信息。有 \(pa = pa_l \cdot pa_r\)，\(pb = pb_l\cdot pb_r\)，\(sab = sab_l + pa_l\cdot sab_r\cdot b_l\)，这个是因为 \(pb_l, B^k\) 都是 \(B\) 的次幂，所以能交换。

Code
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef __int128_t i128;
const int P = 998244353;
struct Matrix {
	int n, m;
	ll a[25][25];
	Matrix() : n(0), m(0) { memset(a, 0, sizeof (a)); }
	Matrix(int _n) : n(_n), m(_n) {
		memset(a, 0, sizeof (a));
		for (int i = 1; i <= n; i++) a[i][i] = 1;
	}
	Matrix(int _n, int _m) : n(_n), m(_m) { memset(a, 0, sizeof (a)); }
	ll *operator[](int x) { return a[x]; }
	const ll *operator[](int x) const { return a[x]; }
	friend Matrix operator+(const Matrix &a, const Matrix &b) {
		Matrix c(a.n, a.m);
		for (int i = 1; i <= a.n; i++)
			for (int j = 1; j <= a.m; j++)
				c[i][j] = (a[i][j] + b[i][j]) % P;
		return c;
	}
	friend Matrix operator*(const Matrix &a, const Matrix &b) {
		Matrix c(a.n, b.m);
		for (int i = 1; i <= a.n; i++)
			for (int k = 1; k <= a.m; k++)
				for (int j = 1; j <= b.m; j++)
					(c[i][j] += a[i][k] * b[k][j]) %= P;
		return c;
	}
};
int m;
struct Node {
	Matrix x, y, s;
	Node() : x(m), y(m), s(m, m) {}
	friend Node operator*(const Node &l, const Node &r) {
		Node res;
		res.x = l.x * r.x;
		res.y = l.y * r.y;
		res.s = l.s + l.x * r.s * l.y;
		return res;
	}
};
Node QPow(Node a, ll b) {
	Node res;
	for (; b; b >>= 1, a = a * a)
		if (b & 1)
			res = res * a;
	return res;
}
Node F(ll n, ll a, ll b, ll c, Node U, Node R) {
	if (!n) return Node();
	if (b >= c) return QPow(U, b / c) * F(n, a, b % c, c, U, R);
	if (a >= c) return F(n, a % c, b, c, U, QPow(U, a / c) * R);
	ll m = ((i128) a * n + b) / c;
	if (!m) return QPow(R, n);
	return QPow(R, (c - b - 1) / a) * U * F(m - 1, c, (c - b - 1) % a, a, R, U) * QPow(R, n - ((i128) c * m - b - 1) / a);
}
ll a, b, c, n;
int main() {
	scanf("%lld%lld%lld%lld%d", &a, &c, &b, &n, &m);
	Node U, R;
	for (int i = 1; i <= m; i++)
		for (int j = 1; j <= m; j++)
			scanf("%lld", &R.x[i][j]), R.s[i][j] = R.x[i][j];
	for (int i = 1; i <= m; i++)
		for (int j = 1; j <= m; j++)
			scanf("%lld", &U.y[i][j]);
	Node ans = F(n, a, b, c, U, R);
	for (int i = 1; i <= m; i++)
		for (int j = 1; j <= m; j++)
			printf("%lld%c", ans.s[i][j], " \n"[j == m]);
	return 0;
}
posted @ 2025-11-26 15:30  Pengzt  阅读(1)  评论(0)    收藏  举报