洛谷 P5170 【模板】类欧几里得算法（类欧几里得）

洛谷题面传送门

一个挺毒瘤的算法，感觉对于应对 NOI 及以下的考试完全没有用（？），不过学了就当一个技能防身吧 /cy

类欧几里是一种能够高效求解带下取整符号，且枚举变量在分子上的和式的算法。其本质上与传统欧几里得算法一点关联都没有，之所以取名为类欧几里得，是因为其复杂度分析与欧几里得相同，均为单次 \(\Theta(\log n)\)。

我们先来看一个类欧几里得最基础的应用：求解 \(\sum\limits_{i=0}^n\lfloor\dfrac{ai+b}{c}\rfloor\)，不妨设其为 \(f(n,a,b,c)\)。

首先特判掉 \(n<0\)​ 的情况——结果显然为 \(0\)​，以及 \(a=0\)​ 的情况，结果为 \((n+1)·\lfloor\dfrac{b}{c}\rfloor\)。

考虑剩余的情况，分 \(a\ge c\) 或 \(b\ge c\)，以及 \(a<c\) 且 \(b<c\) 两种情况处理。

对于 \(a\ge c\) 且 \(b\ge c\)​ 的情况，设 \(a’=\lfloor\dfrac{a}{c}\rfloor,b’=\lfloor\dfrac{b}{c}\rfloor,r_a=a\bmod c,r_b=b\bmod c\)，那么

\[\begin{aligned} f(n,a,b,c)=&\sum\limits_{i=0}^n\lfloor\dfrac{a'i·c+b'c+r_a·i+r_b}{c}\rfloor\\ =&\sum\limits_{i=0}^na'i+b'+\lfloor\dfrac{r_a·i+r_b}{c}\rfloor\\ =&a'·\dfrac{n(n+1)}{2}+b'·(n+1)+f(n,r_a,r_b,c) \end{aligned} \]

递归处理即可。

对于 \(a<c\)​ 且 \(b<c\)​ 的情况，显然有 \(\lfloor\dfrac{an+b}{c}\rfloor<n\)​，不妨设其为 \(M\)​。那么我们就考虑交换下标和值域，即将 \(\lfloor\dfrac{A}{B}\rfloor\)​ 变为 \(\sum\limits_{i=0}^{\lfloor\dfrac{A}{B}\rfloor-1}1\)​，这样有：

\[\begin{aligned} f(n,a,b,c)=&\sum\limits_{i=0}^{n}\sum\limits_{j=0}^{M-1}[j<\lfloor\dfrac{ai+b}{c}\rfloor]\\ =&\sum\limits_{j=0}^{M-1}\sum\limits_{i=0}^n[jc+c-1<ai+b]\\ =&\sum\limits_{j=0}^{M-1}\sum\limits_{i=0}^n[i>\lfloor\dfrac{jc+c-b-1}{a}\rfloor]\\ =&\sum\limits_{j=0}^{M-1}n-\lfloor\dfrac{jc+c-b-1}{a}\rfloor\\ =&\ nM-\sum\limits_{j=0}^{M-1}\lfloor\dfrac{jc+c-b-1}{a}\rfloor\\ =&\ nM-f(M-1,c,c-b-1,a) \end{aligned} \]

可以发现，每次遇到 \(a\ge c\)​ 或 \(b\ge c\)​ 的情况，必然存在一个常数 \(\lambda\approx\dfrac{2}{3}\)，满足 \(a+b\) 至多变为原来的 \(\lambda\) 倍，因此前一种情况最多出现 \(\Theta(\log V)\) 次，而遇到 \(a<c\) 且 \(b<c\) 的情况，必然在进行此次操作以后又变为 \(a\ge c\) 或 \(b\ge c\) 的情况，故后一种情况也最多出现 \(\Theta(\log V)\) 次。也就是说递归层数为 \(\Theta(\log V)\)​ 级别的，进而证明了该算法的复杂度。​

接下来着重研究另外两个函数

\[g(n,a,b,c)=\sum\limits_{i=0}^n\lfloor\dfrac{ai+b}{c}\rfloor^2\\ h(n,a,b,c)=\sum\limits_{i=0}^ni\lfloor\dfrac{ai+b}{c}\rfloor \]

的计算方式。

类比之前 \(f(n,a,b,c)\) 的算法，我们先特判掉一些特殊情况：

• \(n<0\)，\(g(n,a,b,c)=h(n,a,b,c)=0\)。
• \(a=0\)，\(g(n,a,b,c)=(n+1)\lfloor\dfrac{b}{c}\rfloor^2\)，\(h(n,a,b,c)=\dfrac{n(n+1)}{2}·\lfloor\dfrac{b}{c}\rfloor\)。

还是分为 \(a\ge c\) 或 \(b\ge c\) 和 \(a,b<c\) 的情况。

若 \(a\ge c\) 或 \(b\ge c\)，那么还是设 \(a’=\lfloor\dfrac{a}{c}\rfloor,b’=\lfloor\dfrac{b}{c}\rfloor,r_a=a\bmod c,r_b=b\bmod c\)

\[\begin{aligned} g(n,a,b,c)=&\sum\limits_{i=0}^n\lfloor\dfrac{ai+b}{c}\rfloor^2\\ =&\sum\limits_{i=0}^n\lfloor\dfrac{a'i·c+b'·c+r_a·i+r_b}{c}\rfloor^2\\ =&\sum\limits_{i=0}^n(a'i+b'+\lfloor\dfrac{r_a·i+r_b}{c}\rfloor)^2\\ =&\sum\limits_{i=0}^na'^2i^2+b'^2+\lfloor\dfrac{r_a·i+r_b}{c}\rfloor^2+2a'i·\lfloor\dfrac{r_a·i+r_b}{c}\rfloor+2b'·\lfloor\dfrac{r_a·i+r_b}{c}\rfloor+2a'b'i\\ =&\dfrac{n(n+1)(2n+1)}{6}·a'^2+(n+1)·b'^2+2a'b'·\dfrac{n(n+1)}{2}+g(n,r_a,r_b,c)+2a'·h(n,r_a,r_b,c)+2b'·f(n,r_a,r_b,c) \end{aligned} \]

\[\begin{aligned} h(n,a,b,c)=&\sum\limits_{i=0}^ni\lfloor\dfrac{ai+b}{c}\rfloor\\ =&\sum\limits_{i=0}^ni\lfloor\dfrac{a'i·c+b'·c+r_a·i+r_b}{c}\rfloor\\ =&\sum\limits_{i=0}^ni(a'i+b'+\lfloor\dfrac{r_a·i+r_b}{c}\rfloor)\\ =&\sum\limits_{i=0}^na'i^2+b'i+i·\lfloor\dfrac{r_a·i+r_b}{c}\rfloor\\ =&\ a'·\dfrac{n(n+1)(2n+1)}{6}+b'·\dfrac{n(n+1)}{2}+h(n,r_a,r_b,c) \end{aligned} \]

递归下去即可。

若 \(a<c\) 且 \(b<c\)，那么

\[\begin{aligned} g(n,a,b,c)=&\sum\limits_{i=0}^n\lfloor\dfrac{ai+b}{c}\rfloor^2\\ =&\sum\limits_{i=0}^n\sum\limits_{j=0}^{M-1}[j<\lfloor\dfrac{ai+b}{c}\rfloor](2j+1)\\ =&\sum\limits_{j=0}^{M-1}\sum\limits_{i=0}^n[j<\lfloor\dfrac{ai+b}{c}\rfloor]+2·\sum\limits_{j=0}^{M-1}\sum\limits_{i=0}^n[j<\lfloor\dfrac{ai+b}{c}\rfloor]·j\\ =&\ f(n,a,b,c)+2·\sum\limits_{j=0}^{M-1}\sum\limits_{i=0}^n[jc+c-1<ai+b]·j\\ =&\ f(n,a,b,c)+2·\sum\limits_{j=0}^{M-1}j·(n-\lfloor\dfrac{jc+c-b-1}{a}\rfloor)\\ =&\ f(n,a,b,c)+2nM(M-1)-2·\sum\limits_{j=0}^{M-1}j·\lfloor\dfrac{jc+c-b-1}{a}\rfloor\\ =&\ f(n,a,b,c)+2nM(M-1)-2·h(M-1,c,c-b-1,a) \end{aligned} \]

\[\begin{aligned} h(n,a,b,c)=&\sum\limits_{i=0}^n\lfloor\dfrac{ai+b}{c}\rfloor·i\\ =&\sum\limits_{i=0}^n\sum\limits_{j=0}^{M-1}i·[jc+c-1<ai+b]\\ =&\sum\limits_{j=0}^{M-1}\sum\limits_{i=0}^ni·[i>\lfloor\dfrac{jc+c-b-1}{a}\rfloor]\\ =&\sum\limits_{j=0}^{M-1}\dfrac{n(n+1)}{2}-\dfrac{1}{2}·\lfloor\dfrac{jc+c-b-1}{a}\rfloor·(\lfloor\dfrac{jc+c-b-1}{a}\rfloor-1)\\ =&\dfrac{n(n+1)}{2}·M-\dfrac{1}{2}\sum\limits_{j=0}^{M-1}(\lfloor\dfrac{jc+c-b-1}{a}\rfloor^2-\lfloor\dfrac{jc+c-b-1}{a}\rfloor)\\ =&\dfrac{n(n+1)}{2}·M-\dfrac{1}{2}(g(M-1,c,c-b-1,a)-f(M-1,c,c-b-1,a)) \end{aligned} \]

找到了当年推 P5518 [MtOI2019]幽灵乐团 / 莫比乌斯反演基础练习题 的感觉（

直接递归是指数级别的，即 \(2^{\log n}=\mathcal O(n)\)​，会爆，可以使用记忆化搜索将其优化到时空复杂度均 \(\mathcal O(\log n)\)，不过对于此题而言会 MLE。不过注意到三个函数可以并行计算，因此考虑定义一个返回值为 tuple<int, int, int> 的函数一次项返回三个函数的值，这样空间即可做到 \(\mathcal O(1)\)，即可通过。

附：MLE 了的代码：

const int MOD = 998244353;
const int INV2 = 499122177;
const int INV6 = (MOD + 1) / 6;
const int BS1 = 1145141919;
const int BS2 = 1004535809;
const int BS3 = 1000000007;
const int BS4 = 924844033;
int n, a, b, c;
int get1(int x) {return 1ll * x * (x + 1) % MOD * INV2 % MOD;}
int get2(int x) {return 1ll * x * (x + 1) % MOD * (x << 1 | 1) % MOD * INV6 % MOD;}
int clcF(int n, int a, int b, int c);
int clcG(int n, int a, int b, int c);
int clcH(int n, int a, int b, int c);
unordered_map<ll, int> F, G, H;
ll gethash(int n, int a, int b, int c) {
	return (1ll * n * BS1 + 1ll * a * BS2 + 1ll * b * BS3 + 1ll * c * BS4);
}
int clcF(int n, int a, int b, int c) {
	if (n < 0) return 0;
	if (!a) return 1ll * (n + 1) * (b / c) % MOD;
	if (F.count(gethash(n, a, b, c))) return F[gethash(n, a, b, c)];
	if (a >= c || b >= c) return F[gethash(n, a, b, c)] = (1ll * (a / c) * get1(n) + 1ll * (n + 1) * (b / c) + clcF(n, a % c, b % c, c)) % MOD;
	else {
		int M = (1ll * a * n + b) / c;
		return F[gethash(n, a, b, c)] = (1ll * n * M % MOD - clcF(M - 1, c, c - b - 1, a) + MOD) % MOD;
	}
}
int clcG(int n, int a, int b, int c) {
	if (n < 0) return 0;
	if (!a) return 1ll * (n + 1) * (b / c) % MOD * (b / c) % MOD;
	if (G.count(gethash(n, a, b, c))) return G[gethash(n, a, b, c)];
//	printf("G(%d, %d, %d, %d)\n", n, a, b, c);
	if (a >= c || b >= c) {
		int _a = a / c, _b = b / c, ra = a % c, rb = b % c;
		return G[gethash(n, a, b, c)] =
		(1ll * _a * _a % MOD * get2(n) + 1ll * _b * _b % MOD * (n + 1) +
		clcG(n, ra, rb, c) + 2ll * _a * clcH(n, ra, rb, c) % MOD + 2ll * _b * clcF(n, ra, rb, c) % MOD
		+ 2ll * _a * _b % MOD * get1(n)) % MOD;
	} else {
		int M = (1ll * a * n + b) / c;
		return G[gethash(n, a, b, c)] = (2ll * get1(M - 1) * n % MOD - 2ll * clcH(M - 1, c, c - b - 1, a) % MOD + clcF(n, a, b, c) + MOD) % MOD;
	}
}
int clcH(int n, int a, int b, int c) {
	if (n < 0) return 0;
	if (!a) return 1ll * get1(n) * (b / c) % MOD;
	if (H.count(gethash(n, a, b, c))) return H[gethash(n, a, b, c)];
//	printf("H(%d, %d, %d, %d)\n", n, a, b, c);
	if (a >= c || b >= c) {
		int _a = a / c, _b = b / c, ra = a % c, rb = b % c;
		return H[gethash(n, a, b, c)] = (1ll * _a * get2(n) + 1ll * _b * get1(n) + clcH(n, ra, rb, c)) % MOD;
	} else {
		int M = (1ll * a * n + b) / c;
		return H[gethash(n, a, b, c)] = (1ll * M * get1(n) % MOD - 1ll * INV2 * clcG(M - 1, c, c - b - 1, a) % MOD -
		1ll * INV2 * clcF(M - 1, c, c - b - 1, a) % MOD + MOD + MOD) % MOD;
	}
}
void solve() {
	scanf("%d%d%d%d", &n, &a, &b, &c);
	printf("%d %d %d\n", clcF(n, a, b, c), clcG(n, a, b, c), clcH(n, a, b, c));
}
int main() {
	int qu; scanf("%d", &qu);
	while (qu--) solve();
	return 0;
}

AC 代码：

const int MOD = 998244353;
const int INV2 = 499122177;
const int INV6 = (MOD + 1) / 6;
int n, a, b, c;
int get1(int x) {return 1ll * x * (x + 1) % MOD * INV2 % MOD;}
int get2(int x) {return 1ll * x * (x + 1) % MOD * (x << 1 | 1) % MOD * INV6 % MOD;}
tuple<int, int, int> calc(int n, int a, int b, int c) {
	if (n < 0) return mt(0, 0, 0);
	if (!a) return mt(
	1ll * (n + 1) * (b / c) % MOD, 1ll * (n + 1) * (b / c) % MOD * (b / c) % MOD,
	1ll * get1(n) * (b / c) % MOD);
	int S1 = get1(n), S2 = get2(n);
	if (a >= c || b >= c) {
		int _a = a / c, _b = b / c, ra = a % c, rb = b % c;
		tuple<int, int, int> T = calc(n, ra, rb, c);
		return mt((1ll * _a * S1 + 1ll * (n + 1) * _b + get<0>(T)) % MOD,
		(1ll * _a * _a % MOD * S2 + 1ll * _b * _b % MOD * (n + 1) +
		get<1>(T) + 2ll * _a * get<2>(T) % MOD + 2ll * _b * get<0>(T) % MOD
		+ 2ll * _a * _b % MOD * S1) % MOD,
		(1ll * _a * S2 + 1ll * _b * S1 + get<2>(T)) % MOD);
	} else {
		int M = (1ll * a * n + b) / c;
		tuple<int, int, int> T = calc(M - 1, c, c - b - 1, a);
		int F = (1ll * n * M % MOD - get<0>(T) + MOD) % MOD;
		int G = (2ll * get1(M - 1) * n % MOD - 2ll * get<2>(T) % MOD + F + MOD) % MOD;
		int H = (1ll * M * S1 % MOD - 1ll * INV2 * get<1>(T) % MOD -
		1ll * INV2 * get<0>(T) % MOD + MOD + MOD) % MOD;
		return mt(F, G, H);
	}
}
void solve() {
	scanf("%d%d%d%d", &n, &a, &b, &c);
	tuple<int, int, int> T = calc(n, a, b, c);
	printf("%d %d %d\n", get<0>(T), get<1>(T), get<2>(T));
}
int main() {
	int qu; scanf("%d", &qu);
	while (qu--) solve();
	return 0;
}