万能欧几里得算法

类欧几里得算法

类欧几里得算法常用于解决形如

\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}

：

此时

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)$ 。

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

CPP

#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

CPP

#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;
}

万能欧几里得算法

文章操作

类欧几里得算法

万能欧几里得算法

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

Problem

Sol

Code

练习2：

Problem

Code

相关推荐

评论

万能欧几里得算法