欧几里得
欧几里得
全文 绝大多数 内容是对 [0] 中讲述的 粗略抄写 和 胡乱加工
\(\gcd\)，\(\operatorname {lcm}\)
1. 欧几里得算法
即 \(\gcd\)，又称 辗转相除法，用于 递归 得到 最大公因数，时间复杂度 \(O (\log V)\)
当 \(b = 0\) 时 返回 \(a\)（或 \(a = b\) 时 返回）
证明是 简单的，不妨设 \(a > b\)（否则 递归一层后 也会保证 \(a > b\)）
设 \(\gcd (a, b) = c\)，则显然 \(c \mid a, c \mid b\)，有 \(a = kb + r \rightarrow r = a - kb\)
于是同除 \(c\) 有 \(\dfrac {r} {c} = \dfrac {a} {c} - \dfrac {kb} {c}\)，显然有 \(\dfrac {a} {c}, \dfrac {kb} {c} \in \mathbb {Z}\)，于是必有 \(\dfrac {r} {c} \in \mathbb {Z}\)
即 \(c \mid r\)，故一定有 \(\gcd (b, r) = pc, p \ge 1\)
而若 \(p > 1\)，则有 \(pc \mid r, pc \mid b\)，又 \(a = kb + r\)，故 \(pc \mid a\)，与 \(\gcd (a, b) = c\) 矛盾
故 \(p = 1\)，即 \(\gcd (b, r) = c\)，得证
时间复杂度的证明我们 分情况讨论
- \(a > 2b\)，那么 \((b, a \bmod b)\) ，\(b\) 减半，至多递归 \(O (\log V)\) 次
- \(a < 2b\)，两次得到 \((a \bmod b, b \bmod (a \bmod b))\) ，\(b\) 减半，至多递归 \(O (2 \log V)\) 次
- \(a = 2b\)，一次结束
故 时间复杂度 \(O (\log V)\)
代码有 下面两个版本，区别在于 是否使用递归
inline int GCD (const int x, const int y) {
return y ? GCD (y, x % y) : x;
}
上面这个是 递归版，也是 十分经典的版本
inline int Gcd (int x, int y) {
while (y ^= x ^= y ^= x %= y);
return x;
}
上面这个是 非递归版，前面的 \(3\) 个 ^= 其实就是起到 swap 的作用，这样就好理解了
但这个东西 不够快，考虑 二进制位运算 来优化这个问题，这样的方法又称 \(Binary ~ GCD\)
具体而言，我们可以把 \(a, b\) 按 奇偶 分成 如下 三种情况（以下默认 \(b < a\)）
-
\(a, b\) 均为 偶数，显然有 公因数 \(2\)，故而可以变成 \(\gcd (\dfrac {a} {2}, \dfrac {b} {2})\)，最后答案 乘上 \(2\)
-
\(a, b\) 均为 奇数，由于 \(\gcd (a, b) = \gcd (b, a \bmod b)\)，故显然 \(\gcd (a, b) = \gcd (|a - b|, b)\)
-
\(a, b\) 一奇 一偶，其不可能有 偶公因数，故把 偶数 \(\div 2\) 即可，即 \(\gcd (\dfrac {a} {2}, b)\) 或 \(\gcd (a, \dfrac {b} {2})\)
进而容易发现，每次 \(\div 2\) 时，我们会把 对应偶数 的 \(2\) 除尽
故我们可以使用 __builtin_ctz，\(O(1)\) 找到 末尾 \(0\) 的个数，来优化常数
inline int BGCD (int x, int y) {
#define ctz __builtin_ctz
int a = ctz (x), b = ctz (y), c = min (a, b), d;
y >>= b;
while (x) x >>= a, d = x - y, a = ctz (d), y = min (x, y), x = abs (d);
return y << c;
}
实现的较为精细的话应当能通过 Luogu P5435 基于值域预处理的快速 GCD 这个题
CF1806F2 GCD Master (hard version)
前置任务 CF1806F1 GCD Master (easy version)
先考虑 所有数不同的情况
根据题目，做一次操作就是 将两数合并成它们的 \(\gcd\)，容易发现这个东西 可以扩展
即对 \(x, y, z\) 依次操作（\(\gcd (\gcd (x, y), z)\)）相当于 将三个数合并成它们的 \(\gcd\)
若存在集合 \(S\)，以下定义 \(\gcd (S)\) 为 整个集合的最大公因数
继续推广，那么原题意实际上求的是 一种分组方案，将原序列分成 若干 (\(R\)) 组（\(P_i\)）
使得 \(\sum \limits_{i = 1} ^ {R - 1} |P_i| - 1 = k\)，同时使 \(\sum \limits_{i = 1} ^ {R - 1} \gcd (P_i) + \sum \limits_{p \in P_R} p\) 尽量大
第 \(R\) 组即 原序列 没有被拿去求 \(\gcd\) 的数（剩下来的）
不难看出，当 分的组数越多，第 \(R\) 组的 大小就会越小（每组会 “损失” 一个）
而第 \(R\) 组的 元素和 会 直接贡献到最终答案，贡献很大
故我们可以猜想，需要 分的组数尽量少，比如 只取一组（大小为 \(k + 1\)），然后尝试证明
考虑 反证法
若存在两组 \(S, T\)，显然其 区间和 贡献是 \(\gcd (S) + \gcd (T) - \sum \limits_{i \in S} i - \sum \limits_{i \in T} i\)
若将其 合并成一组，则显然 \(\gcd\) 部分贡献变成 \(\gcd (\gcd (S), \gcd (T))\)
而由于要求 操作次数恒定，故 合并成一组之后 大小应当减少 \(1\)
即我们可以 丢一个数去 ”第 \(R\) 组“，直接加上这个数的贡献，可以想到我们要丢 较大的数
显然，\(\gcd (S)\) 随 \(|S|\) 的 减小（\(|S| > 3\)）而 单调不降
不用担心 丢了一个数 会对 \(\gcd\) 做 负贡献
于是贡献变成 \(\gcd (\gcd (S), \gcd (T)) + E - \sum \limits_{i \in S} i - \sum \limits_{i \in T} i\)（\(E\) 是 丢的数）
不妨设 \(\gcd (S) \ge \gcd (T)\)，由于 \(S, T\) 中 均没有重复的数
故有 \(\operatorname {max\_element} (S) \ge 2 \gcd (S)\)，于是令 \(E = \operatorname {max\_element} (S)\)
容易发现 \(\gcd (\gcd (S), \gcd (T)) + E - \sum \limits_{i \in S} i - \sum \limits_{i \in T} i > \gcd (S) + \gcd (T) - \sum \limits_{i \in S} i - \sum \limits_{i \in T} i\)
即 合并成一组一定更优
原题就变成了 在原序列中 找到 \(k + 1\) 个数，使得 这些数的 \(\gcd\) 加上 没有被选的数 和 最大
但是上面都是建立在 原序列不重复 的前提上的，但实际上 原序列可以重复
考虑 如何处理重复元素，事实上这种东西 性质很好
因为 对重复元素操作 完全等价于 删去其中一个数（不用考虑 \(\gcd\)）
于是我们可以把 重复的数 单独拿出来（一种数有 \(P\) 个，拿出 \(P - 1\) 个），从小到大 排序
容易求得 前缀和，设这样的数 一共 有 \(Cnt\) 个
于是后续 在不可重序列 上的操作只需要做 \([k - Cnt, k]\) 次，再补上这部分 前缀和 即可
对于 \(\sum M < 10 ^ 6\) 的 煎蛋版本，有一个 \(O (M \log M)\) 的 经典做法
我们 枚举 \(\gcd\) 的值，然后 枚举其倍数（在值域 \([1, M]\) 内），以此统计 能整除该 \(\gcd\) 的数个数
当 这个个数 \(\in [k - Cnt, k]\) 时，我们即可 得到可行解，最终答案 即是 所有可行解 的 \(\max\)
枚举倍数 这个过程代价是 调和级数（\(\ln\)），通常我们可以也视作 \(\log\)
#include <bits/stdc++.h>
const int MAXN = 1000005;
using namespace std;
int T;
int N, M, K;
long long Ret, Ans, Sum;
long long P[MAXN];
int B[MAXN], C[MAXN];
int Cnt, tot;
inline void Solve () {
Cnt = tot = 0, Ans = Sum = 0;
cin >> N >> M >> K, ++ K;
for (int i = 1; i <= N; ++ i) {
cin >> B[i];
if (C[B[i]]) P[++ Cnt] = B[i];
C[B[i]] ++, Sum += B[i];
}
sort (P + 1, P + Cnt + 1);
for (int i = 2; i <= Cnt; ++ i) P[i] += P[i - 1];
if (Cnt >= K - 1) Ans = max (Ans, Sum - P[Cnt]);
for (int i = 1; i <= M; ++ i) {
Ret = i, tot = 0;
for (int p = i; p <= M && tot < K; p += i)
if (C[p])
++ tot, Ret -= p, (tot >= K - Cnt) ? Ans = max (Ans, (Sum + Ret - P[K - tot])) : 0;
}
for (int i = 1; i <= N; ++ i) C[B[i]] = 0, B[i] = 0;
cout << Ans << '\n';
}
int main () {
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
cin >> T;
while (T --) Solve ();
return 0;
}
考虑 坚硬的 版本，显然我们不能枚举 \(\gcd\) 了
但是在刚刚枚举 \(\gcd\) 的 思路 中，我们不难？发现一些事实
比如，当 \(\gcd\) 确定的时候，选取 最小的几个数字 一定是最优的
虽然当 \(\gcd\) 不固定时，直接选取 可行最小的数 可能导致 \(\gcd\) 值 大幅下降
但这不妨碍我们将其 作为一种 调整答案的思路，于是我们可以考虑 构造调优 的 思想
先来一组 可行解，再 调整得到 最优解
设 当前我们选中 用来求 \(\gcd\) 的集合 \(S\)，剩余数集 \(P_R\)，\(S\) 最大公约数 \(x = \gcd (S)\)
考虑 如何让 \(\gcd (S) + \sum \limits_{i \in P_R} i\) 变大，根据上面的思路，将 \(P_R\) 中 最小的数 \(p\) 加入 \(S\)
这会使得 \(\gcd (S') = \gcd (x, p)\)，此时我们需要从 \(S'\) 中 拿出一个数（需要保证 \(|S|\) 不变）
由于 删数不会使 \(gcd\) 变小，故为了使 \(\sum \limits_{i \in P_R} i\) 部分更大，我们应拿出 \(S'\) 中 最大的数 \(q\)
容易发现，这样对 最终和 的 贡献 为 \(q - p + gcd (x, p) - x\)
此时当 \(S'\) 中存在另一个数 \(r > p ~ (r < q)\) 时，\(q - p + gcd (x, p) - x > 0\)
由于 \(\gcd (S) = x, r < q\)，故 \(q - r \ge x\)，则 \(q - p + gcd (x, p) - x > r - p + gcd (x, p)\)
又 \(r > p\)，故 得证
即 满足上述条件（\(S\) 中 存在两个数 大于 剩余最小值）时，这样调整一定更优
整理就可以得到结论，即 选择的 \(k + 1\) 个数 中 必定包含前 \(k\) 小的数
进而可以知道，最后选择的数必然是 前面 \(t\) 种数 + 重复的 \(k - t - 1\) 个数 + \(1\) 个后面的数
一个后面的数 来源于 调整条件 需要 \(S'\) 中存在 另一个数 大于 剩余最小值
我们对 剩余的数 排序，然后做 前缀 \(\gcd\)，可知前缀 \(\gcd\) 单调递减 且 只有 \(\log V\) 个 拐点
于是先 枚举 \(t\) 在哪一段，然后枚举 后面的数是哪个
最后枚举 这一段选了多少个数 \(\to\) 选了多少重复的数 即可，时间复杂度 \(O (N \log M)\)
快读是用来 输出
__int128的
#include <bits/stdc++.h>
const int MAXN = 1000005;
using namespace std;
namespace Fastio {
#define USE_FASTIO 1
#define IN_LEN 450000
#define OUT_LEN 450000
char ch, c; int len;
short f, top, s;
inline char Getchar() {
static char buf[IN_LEN], *l = buf, *r = buf;
if (l == r) r = (l = buf) + fread(buf, 1, IN_LEN, stdin);
return (l == r) ? EOF : *l++;
}
char obuf[OUT_LEN], *ooh = obuf;
inline void Putchar(char c) {
if (ooh == obuf + OUT_LEN) fwrite(obuf, 1, OUT_LEN, stdout), ooh = obuf;
*ooh++ = c;
}
inline void flush() { fwrite(obuf, 1, ooh - obuf, stdout); }
#undef IN_LEN
#undef OUT_LEN
struct Reader {
template <typename T> Reader& operator >> (T &x) {
x = 0, f = 1, c = Getchar();
while (!isdigit(c)) { if (c == '-') f *= -1; c = Getchar(); }
while ( isdigit(c)) x = (x << 3) + (x << 1) + (c ^ 48), c = Getchar();
x *= f;
return *this;
}
Reader() {}
} cin;
const char endl = '\n';
struct Writer {
typedef long long mxdouble;
template <typename T> Writer& operator << (T x) {
if (x == 0) { Putchar('0'); return *this; }
if (x < 0) Putchar('-'), x = -x;
static short sta[40];
top = 0;
while (x > 0) sta[++top] = x % 10, x /= 10;
while (top > 0) Putchar(sta[top] + '0'), top--;
return *this;
}
Writer& operator << (const char *str) {
int cur = 0;
while (str[cur]) Putchar(str[cur++]);
return *this;
}
inline Writer& operator << (char c) {Putchar(c); return *this;}
Writer() {}
~ Writer () {flush();}
} cout;
#define cin Fastio::cin
#define cout Fastio::cout
#define endl Fastio::endl
}
inline long long BGCD (long long x, long long y) {
#define ctz __builtin_ctzll
int a = ctz (x), b = ctz (y), c = min (a, b);
long long d;
y >>= b;
while (x) x >>= a, d = x - y, a = ctz (d), y = min (x, y), x = abs (d);
return y << c;
}
map <long long, int> Mp;
int E[MAXN];
long long A[MAXN], B[MAXN], G[MAXN];
__int128 P[MAXN], S[MAXN];
__int128 Sum = 0, Ret = 0, Ans = 0;
long long M;
int N, K;
int tot, Cnt, tmp;
inline void Solve () {
Sum = 0, Ans = 0, tot = 0, Cnt = 0, tmp = 0, Mp.clear ();
cin >> N >> M >> K, ++ K;
for (int i = 1; i <= N; ++ i) {
cin >> B[i];
if (Mp.count (B[i])) P[++ Cnt] = B[i];
Mp[B[i]] = 1, Sum += B[i];
}
sort (P + 1, P + Cnt + 1);
for (auto i : Mp) A[++ tot] = i.first;
for (int i = 2; i <= Cnt; ++ i) P[i] += P[i - 1];
for (int i = 1; i <= tot; ++ i) S[i] = S[i - 1] + A[i];
if (Cnt >= K - 1) Ans = max (Ans, Sum - P[K - 1]);
for (int i = 1; i <= tot; ++ i) {
G[i] = BGCD (G[i - 1], A[i]);
if (G[i] != G[i - 1]) E[++ tmp] = i;
if (K - Cnt <= i && i <= K) Ans = max (Ans, Sum - S[i] - P[K - i] + G[i]);
}
E[tmp + 1] = tot + 1;
for (int i = 1; i <= tmp; ++ i) {
Ret = - Sum;
for (int p = E[i + 1]; p <= tot; ++ p) Ret = max (Ret, (__int128) BGCD (G[E[i]], A[p]) - A[p]);
for (int p = E[i]; p < min (E[i + 1], K); ++ p)
if (K - Cnt - 1 <= p) Ans = max (Ans, Sum + Ret - S[p] - P[K - p - 1]);
}
cout << Ans << '\n';
}
int T;
int main () {
cin >> T;
while (T --) Solve ();
return 0;
}
2. 扩展欧几里得算法
即 \(\operatorname {exgcd}\)，用于求 \(ax + by = \gcd (a, b)\) 的 不定方程 的 一组整数解 \(x, y\)
容易发现，当 \(\gcd\) 进行到 最后返回（\(b = 0\)） 的时候
对于当时的 \(a, b\)，上述方程一定有解 \(x = 1, y = 0\)（当时 \(a = \gcd (a, b), b = 0\)）
显然在 \(\gcd\) 整个过程中，\(\gcd (a, b)\) 均不变，即 方程右边恒等
于是有
然后我们就可以在 \(\gcd\) 的过程中 递归求解 答案了
inline int EXGCD (const int a, const int b, int &x, int &y) {
if (b == 0) return x = 1, y = 0, a;
int d = EXGCD (b, a % b, y, x);
y -= a / b * x; return d;
}
由于需要保证 方程右边 \(\gcd\) 不变，故此处 二进制优化 不适用（无意义）
此处我们得到的是 \(ax + by = \gcd (a, b)\) 的一组 特解 \(x = x_0, y = y_0\)
其 不具有什么特殊性质，如果需要得到一些 特殊的解，可以使用下面的 通解 进一步计算
其中 \(x_0, y_0\) 就是我们得到的 特解，且 \(a' = \dfrac {a} {\gcd (a, b)}, b'\ = \dfrac {b} {\gcd (a, b)}\)
关于 \(ax + by = c\) 解的 存在性 讨论
裴蜀定理：如果 \(a, b \in \mathbb {Z}\)，则有 \(x, y \in \mathbb {Z}\) 满足 \(ax + by = \gcd (a, b)\)
\(\operatorname {exgcd}\) 的算法流程 即在 构造上 给出了该定理的证明
故当 \(\gcd (a, b) | c\) 时，该不定方程 有解，反之则 无解
Luogu P5656 【模板】二元一次不定方程 (exgcd)
一道 没有人做 的 板子，主要是 要求太多了，给出 不定方程 \(ax + by = c\)
- 无整数解，输出
-1 - 无正整数解，输出 \(x_{\min}, y_{\min} \in \mathbb {N ^ *}\)（无需保证对应情况下的 \(y_{\max}, x_{\max} \in \mathbb {N ^ *}\)）
- 有正整数解，输出 解数， \(x_{\min}, y_{\min}, x_{\max}, y_{\max} \in \mathbb {N ^ *}\)（需保证对应情况下 \(x, y \in \mathbb {N ^ *}\)）
有无解判一下 \(\gcd (a, b) \mid c\) 是否成立 即可
成立时通过 \(\operatorname {exgcd}\) 即可得到一组 特解 \(x = x_0, y = y_0\)
注意直接 \(\operatorname {exgcd}\) 得到的是 \(ax + by = \gcd (a, b)\) 得解，需要乘一个 \(\dfrac {c} {\gcd (a, b)}\)
显然，存在通解 \(x = x_0 + kb', y = y_0 - ka'\)，其中 \(a' = \dfrac {a} {\gcd (a, b)}, b'\ = \dfrac {b} {\gcd (a, b)}\)
于是显然，\(x, y\) 变化方向相反（即 \(x\) 增大时，\(y\) 减小，反之亦然）
于是当 \(x\) 取 最小正整数解 时，\(y\) 取到范围（\(x \in \mathbb {N ^ *}\)）内 最大解
此时若 \(y < 0\)，则 原方程无正整数解，反之则此时为 \(y\) 的 最大（要求下）正整数解
同理我们可以对 \(y\) 进行操作，使之取到 最小正整数解，得到 \(x\) 的 最大解，同上
注意怎么取到 最小正整数解，以 \(x\) 为例
若 \(x_0 > 0\)，\(x_0\) 减去一个 \(b' \left \lfloor \dfrac {x_0 - 1} {b'} \right \rfloor\) 即可（要求为 正，不能取 \(0\)，故需要 \(-1\)）
若 \(x_0 \le 0\)，\(x_0\) 减去一个 \(b' (\left \lfloor \dfrac {x} {b'} \right \rfloor - 1)\) 即可
然后就行了
#include <bits/stdc++.h>
const int MAXN = 100005;
inline int EXGCD (const int a, const int b, long long &x, long long &y) {
if (b == 0) return x = 1, y = 0, a;
int d = EXGCD (b, a % b, y, x);
y -= 1ll * a / b * x; return d;
}
int a, b, c;
long long x, y, p, q, d;
long long xtmp, ytmp, xmin, xmax, ymin, ymax;
bool ns;
using namespace std;
int N;
int main () {
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
cin >> N;
for (int i = 1; i <= N; ++ i) {
cin >> a >> b >> c, ns = 1;
if (c % (d = EXGCD (a, b, x, y)) != 0) {
cout << -1 << '\n';
continue ;
}
x = x * (c / d), y = y * (c / d);
a = a / d, b = b / d;
xtmp = x, ytmp = y;
if (x <= 0) p = x / b - 1;
else p = (x - 1) / b;
if (y <= 0) q = y / a - 1;
else q = (y - 1) / a;
x -= p * b, y += p * a;
if (y <= 0) ns = 0;
xmin = x, ymax = y;
x = xtmp, y = ytmp;
x += q * b, y -= q * a;
xmax = x, ymin = y;
if (ns) cout << abs (p + q) + 1 << ' ' << xmin << ' ' << ymin << ' ' << xmax << ' ' << ymax << '\n';
else cout << xmin << ' ' << ymin << '\n';
}
return 0;
}
Luogu P1082 [NOIP2012 提高组] 同余方程
求关于 \(x\) 的 同余方程 \(ax \equiv 1 \pmod b\) 的 最小正整数解
即 \(ax + by = 1\) 的 \(x\) 最小正整数解，直接做，随便做
这才是 真正的模板题
#include <bits/stdc++.h>
using namespace std;
inline int EXGCD (const int a, const int b, int &x, int &y) {
if (b == 0) return x = 1, y = 0, a;
int d = EXGCD (b, a % b, y, x);
y -= a / b * x; return d;
}
int a, b, x, y;
int main () {
cin >> a >> b;
EXGCD (a, b, x, y), x = (x + b) % b;
cout << x << '\n';
return 0;
}
[ABC340F] S = 1
容易发现，三角形 恒满足 \(S = \dfrac {ah} {2}\)，我们设 源点 \(O ~ (0, 0)\)，已知点 \(P (x_0, y_0)\)
则所求点 解析式 易得为 \(y = \dfrac {y_0} {x_0} x \pm \dfrac {2} {x_0}\)，根据题意我们需要使得 \(x, y\) 均为整数
故 \(x_0 \mid (xy_0 \pm 2)\)，改写为 \(x y_0 - yx_0 = \pm 2\)，直接 \(\operatorname {exgcd}\) 即可
注意符号，虽然 \(x, y \in \mathbb {Z}\)，但 不代表可以任意变号
#include <bits/stdc++.h>
using namespace std;
inline long long EXGCD (const long long a, const long long b, long long &x, long long &y) {
if (b == 0) return x = 1, y = 0, a;
long long d = EXGCD (b, a % b, y, x);
y -= a / b * x; return d;
}
long long a, b, x, y, d;
int main () {
cin >> a >> b;
d = EXGCD (b, a, x, y), y = -y, d = abs (d);
if (d != 1 && d != 2) return cout << -1 << '\n', 0;
if (d == 1) x <<= 1, y <<= 1;
cout << x << ' ' << y << '\n';
return 0;
}
3. 类欧几里得算法
用 类似欧几里得算法 辗转相除 的思想 在 \(O(\log V)\) 结构时间 下来 求解一些式子 的算法
例下式，\(a, b, c, n\) 均为常数
在 [0] 中，提到这个式子与《具体数学》\(3.5 ~ eg.3\) 很像，且在 \(n = c\) 时，存在封闭形式
但是笔者 没有看到那儿，故此处不表
不妨设 \(a, b > c\)，于是可以通过 取模 限制值域
容易发现，前面两部分可以 \(O(1)\) 得出，后面的 \(a, b\) 值域减小 至小于 \(c\)
于是我们可以保证将 任意上述式子 \(O(1)\) 转化至 \(a, b < c\) 的情况
接着完成 后面的部分，使之变成 可递归式（以下保证 \(a, b < c\)）
不妨设 \(m = \left \lfloor \dfrac {an + b} {c} \right \rfloor, t = \left \lfloor \dfrac {cj + c - b - 1} {a} \right \rfloor\)，于是
\(OK\) 啊，这就是一个 递归式 了，交换 \(a, c\)，于是可以重复上述过程，取模 + 递归
这样的流程和 \(\gcd\) 的 辗转相除法 相同，故得名，同时 时间复杂度也是一样的 \(O (\log V)\)
Luogu P5171 Earthquake
显然答案为 \(\sum \limits _ {i = 0} ^ {\left \lfloor \frac {a} {c} \right \rfloor} {(\left \lfloor \dfrac {c - ai} {b} \right \rfloor + 1)}\)
很像啊，很像啊，但是 负号 很丑，于是加一个 \(bi\) 上去（保证 \(b > a\)）
于是 第三部分 就是 类欧板子（\(n = \left \lfloor \dfrac {a} {c} \right \rfloor, a = b - a, b = c, c = b\)）即可
#include <bits/stdc++.h>
using namespace std;
inline long long Solve (const long long a, const long long b, const long long c, const long long n) {
long long ac = a / c, bc = b / c, m = (a * n + b) / c;
if (a == 0) return bc * (n + 1);
if (a >= c || b >= c) return Solve (a % c, b % c, c, n) + ac * n * (n + 1) / 2 + bc * (n + 1);
return n * m - Solve (c, c - b - 1, a, m - 1);
}
long long a, b, c;
int main () {
cin >> a >> b >> c;
if (b < a) swap (a, b);
cout << Solve (b - a, c, b, c / a) - (c / a) * (c / a + 1) / 2 + (c / a) + 1 << '\n';
return 0;
}
Luogu P5170 【模板】类欧几里得算法
同理求出 \(g (a, b, c, n) = \sum \limits _ {i = 0} ^ {n} i \left \lfloor \frac {ai + b} {c} \right \rfloor\)，\(h (a, b, c, n) = \sum \limits _ {i = 0} ^ {n} \left \lfloor \frac {ai + b} {c} \right \rfloor ^ 2\) 两和式值即可
以下推导
\[\begin {aligned} & g (a, b, c, n) \\ =& \sum \limits _ {i = 0} ^ {n} { i \left \lfloor \dfrac {ai + b} {c} \right \rfloor } \\ =& \sum \limits _ {i = 0} ^ {n} { \left ( i \left ( \left \lfloor \dfrac {a} {c} \right \rfloor i + \left \lfloor \dfrac {b} {c} \right \rfloor + \left \lfloor \dfrac {(a \bmod c) i + b \bmod c} {c} \right \rfloor \right ) \right ) } \\ =& \left \lfloor \dfrac {a} {c} \right \rfloor { \dfrac {n (n + 1) (2n + 1)} {6} } + \left \lfloor \dfrac {b} {c} \right \rfloor { \dfrac {n (n + 1)} {2} } + g (a \bmod c, b \bmod c, c, n) \end {aligned} \]考虑 \(a, b < c\)，设 \(m = \left \lfloor \dfrac {an + b} {c} \right \rfloor, t = \left \lfloor \dfrac {jc + c - b - 1} {a} \right \rfloor\)，于是
\[\begin {aligned} & \sum \limits _ {i = 0} ^ {n} { \left ( i \left \lfloor \dfrac {ai + b} {c} \right \rfloor \right ) } \\ =& \sum \limits _ {i = 0} ^ {n} { \sum \limits _ {j = 0} ^ { \left \lfloor \frac {ai + b} {c} \right \rfloor - 1 } i } \\ =& \sum \limits _ {j = 0} ^ {m - 1} { \sum \limits _ {i = 0} ^ {n} { \left [ j < \left \lfloor \dfrac {ai + b} {c} \right \rfloor \right ] i } } \\ =& \sum \limits _ {j = 0} ^ {m - 1} { \sum \limits _ {i = 0} ^ {n} { [i > t] i } } \\ =& \sum \limits _ {j = 0} ^ {m - 1} { \dfrac {(n - t) (n + t + 1)} {2} } \\ =& \sum \limits _ {j = 0} ^ {m - 1} { \dfrac {n ^ 2 + n - t ^ 2 - t} {2} } \\ =& \dfrac {1} {2} \left ( mn (n + 1) - h (c, c - b - 1, a, m - 1) - f (c, c - b - 1, a, m - 1) \right ) \end {aligned} \]
\[\begin {aligned} & h (a, b, c, n) \\ =& \sum \limits _ {i = 0} ^ {n} { \left \lfloor \dfrac {ai + b} {c} \right \rfloor ^ 2 } \\ =& \sum \limits _ {i = 0} ^ {n} { \left \lfloor \left \lfloor \dfrac {a} {c} \right \rfloor i + \left \lfloor \dfrac {b} {c} \right \rfloor + \dfrac {(a \bmod c) i + b \bmod c} {c} \right \rfloor ^ 2 } \\ =& \sum \limits _ {i = 0} ^ {n} { \left ( \left \lfloor \dfrac {a} {c} \right \rfloor ^ 2 i ^ 2 + \left \lfloor \dfrac {b} {c} \right \rfloor ^ 2 + \left \lfloor \dfrac {(a \bmod c) i + b \bmod c} {c} \right \rfloor ^ 2 + 2 \left \lfloor \dfrac {(a \bmod c) i + b \bmod c} {c} \right \rfloor \left \lfloor \dfrac {a} {c} \right \rfloor i + 2 \left \lfloor \dfrac {(a \bmod c) i + b \bmod c} {c} \right \rfloor \left \lfloor \dfrac {b} {c} \right \rfloor + 2 \left \lfloor \dfrac {b} {c} \right \rfloor \left \lfloor \dfrac {a} {c} \right \rfloor i \right ) } \\ =& \left \lfloor \dfrac {a} {c} \right \rfloor ^ 2 { \dfrac {n (n + 1) (2n + 1)} {6} } + \left \lfloor \dfrac {b} {c} \right \rfloor ^ 2 (n + 1) + \sum \limits _ {i = 0} ^ {n} { \left \lfloor \dfrac {(a \bmod c) i + b \bmod c} {c} \right \rfloor ^ 2 } + 2 \sum \limits _ {i = 0} ^ {n} { \left ( \left \lfloor \dfrac {(a \bmod c) i + b \bmod c} {c} \right \rfloor \left \lfloor \dfrac {a} {c} \right \rfloor i + \left \lfloor \dfrac {(a \bmod c) i + b \bmod c} {c} \right \rfloor \left \lfloor \dfrac {b} {c} \right \rfloor + \left \lfloor \dfrac {b} {c} \right \rfloor \left \lfloor \dfrac {a} {c} \right \rfloor i \right ) } \\ =& \left \lfloor \dfrac {a} {c} \right \rfloor ^ 2 { \dfrac {n (n + 1) (2n + 1)} {6} } + \left \lfloor \dfrac {b} {c} \right \rfloor ^ 2 (n + 1) + \left \lfloor \dfrac {b} {c} \right \rfloor \left \lfloor \dfrac {a} {c} \right \rfloor n (n + 1) + h (a \bmod c, b \bmod c, c, n) + \sum \limits _ {i = 0} ^ {n} { \left ( 2 \left ( \left \lfloor \dfrac {a} {c} \right \rfloor i + \left \lfloor \dfrac {b} {c} \right \rfloor \right ) \left \lfloor \dfrac {(a \bmod c) i + b \bmod c} {c} \right \rfloor \right ) } \\ =& \left \lfloor \dfrac {a} {c} \right \rfloor ^ 2 { \dfrac {n (n + 1) (2n + 1)} {6} } + \left \lfloor \dfrac {b} {c} \right \rfloor ^ 2 (n + 1) + \left \lfloor \dfrac {b} {c} \right \rfloor \left \lfloor \dfrac {a} {c} \right \rfloor n (n + 1) + h (a \bmod c, b \bmod c, c, n) + 2 \left \lfloor \dfrac {a} {c} \right \rfloor g (a \bmod c, b \bmod c, c, n) + 2 \left \lfloor \dfrac {b} {c} \right \rfloor f (a \bmod c, b \bmod c, c, n) \end {aligned} \]考虑 \(a, b < c\)，但是 平方 很令人火大，转化
\[n ^ 2 = (2 \sum \limits _ {i = 0} ^ {n} i) - n \]于是
\[\begin {aligned} & h (a, b, c, n) \\ =& \sum \limits _ {i = 0} ^ {n} { \left \lfloor \dfrac {ai + b} {c} \right \rfloor ^ 2 } \\ =& \sum \limits _ {i = 0} ^ {n} { \left ( 2 \sum \limits _ {j = 0} ^ { \left \lfloor \frac {ai + b} {c} \right \rfloor } j - \left \lfloor \dfrac {ai + b} {c} \right \rfloor \right ) } \\ =& \sum \limits _ {i = 0} ^ {n} { \left ( 2 \sum \limits _ {j = 0} ^ { \left \lfloor \frac {ai + b} {c} \right \rfloor } j \right ) - f (a, b, c, n) } \end {aligned} \]设 \(m = \left \lfloor \dfrac {an + b} {c} \right \rfloor, t = \left \lfloor \dfrac {jc + c - b - 1} {a} \right \rfloor\)，于是
\[\begin {aligned} & 2 \sum \limits _ {i = 0} ^ {n} { \left ( \sum \limits _ {j = 0} ^ { \left \lfloor \frac {ai + b} {c} \right \rfloor } j \right ) } \\ =& 2 \sum \limits _ {i = 0} ^ {n} { \left ( \sum \limits _ {j = 1} ^ { \left \lfloor \frac {ai + b} {c} \right \rfloor } j \right ) } \\ =& 2 \sum \limits _ {i = 0} ^ {n} { \left ( \sum \limits _ {j = 0} ^ { \left \lfloor \frac {ai + b} {c} \right \rfloor - 1 } (j + 1) \right ) } \\ =& 2 \sum \limits _ {j = 0} ^ {m - 1} { (j + 1) \sum \limits _ {i = 0} ^ {n} { \left [ j < \left \lfloor \dfrac {ai + b} {c} \right \rfloor \right ] } } \\ =& 2 \sum \limits _ {j = 0} ^ {m - 1} { (j + 1) \sum \limits _ {i = 0} ^ {n} {[i > t]} } \\ =& 2 \sum \limits _ {j = 0} ^ {m - 1} {(j + 1) (n - t)} \\ =& nm (m + 1) + 2 \sum \limits _ {j = 0} ^ {m - 1} (j + 1) t \\ =& nm (m + 1) - 2 (g (c, c - b - 1, a, m - 1) + f (c, c - b - 1, a, m - 1)) \end {aligned} \]
于是 三个合在一起 递归求解 即可
注意 各个位置上的 取模，记得 \(h\) 最后要减一个 \(f (a, b, c, n)\)
而 \(m\) 不能去模，否则 递归下去 的 值域就变了
#include <bits/stdc++.h>
const int MOD = 998244353;
using namespace std;
int N;
struct Node {
int f, g, h;
} Ans;
const int I2 = 499122177, I6 = 166374059;
inline Node Solve (const int a, const int b, const int c, const int n) {
Node pre = {0, 0, 0}, ans = {0, 0, 0};
int ac = a / c, bc = b / c, m = (1ll * a * n + b) / c;
if (a == 0) {
ans.f = 1ll * bc * (n + 1) % MOD;
ans.g = 1ll * bc * n % MOD * (n + 1) % MOD * I2 % MOD;
ans.h = 1ll * bc * bc % MOD * (n + 1) % MOD;
return ans;
} else if (a >= c || b >= c) {
pre = Solve (a % c, b % c, c, n);
ans.f = ((pre.f + 1ll * (1ll * n * (n + 1) % MOD * I2) % MOD * ac % MOD) % MOD + 1ll * (n + 1) * bc % MOD) % MOD;
ans.g = ((pre.g + 1ll * (1ll * n * (n + 1) % MOD * (2 * n + 1) % MOD) * I6 % MOD * ac % MOD) % MOD + 1ll * (1ll * n * (n + 1) % MOD * I2) % MOD * bc % MOD) % MOD;
ans.h = ((pre.h + 1ll * (1ll * n * (n + 1) % MOD * (2 * n + 1) % MOD) * I6 % MOD * ac % MOD * ac % MOD) % MOD + 1ll * bc * bc % MOD * (n + 1) % MOD + 1ll * (a / c) * (b / c) % MOD * n % MOD * (n + 1) % MOD + 1ll * 2 * ac * pre.g % MOD + 1ll * 2 * bc * pre.f % MOD) % MOD;
return ans;
}
pre = Solve (c, c - b - 1, a, m - 1);
ans.f = (1ll * m * n % MOD - pre.f + MOD) % MOD;
ans.g = (1ll * m * n % MOD * (n + 1) % MOD - pre.h - pre.f + MOD + MOD) * I2 % MOD;
ans.h = (1ll * m * n % MOD * (m + 1) % MOD - 2ll * pre.g % MOD - 2ll * pre.f % MOD - ans.f + MOD + MOD + MOD) % MOD;
return ans;
}
int a, b, c, n;
int main () {
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
cin >> N;
for (int i = 1; i <= N; ++ i) {
cin >> n >> a >> b >> c;
Ans = Solve (a, b, c, n);
cout << Ans.f << ' ' << Ans.h << ' ' << Ans.g << '\n';
}
return 0;
}
Luogu P5179 Fraction
给定 \(\dfrac {a} {b}, \dfrac {c} {d}\)，求 \(p, q\) 使得 \(\dfrac {a} {b} < \dfrac {p} {q} < \dfrac {c} {d}\)，且 \(p, q\) 尽量小
这题 很有意思，我们不妨设 \(\dfrac {a} {b}, \dfrac {c} {d}\) 均为 真分数
然后，去倒数，我们可以得到 \(\dfrac {b} {a} < \dfrac {q} {p} < \dfrac {d} {c}\)，此时 \(\dfrac {b} {a}, \dfrac {d} {c}\) 均为 假分数
我们减去 较小数 的 整数部分，于是变成
如果这个时候 左右 为 真分数，继续递归 即可
若 一真一假，则 取 \(p = q - p \left \lfloor \dfrac {b} {a} \right \rfloor = 1\) 作为答案 返回计算即可
时间复杂度 是 \(O (\log V)\) 的，这点 很容易证明
#include <bits/stdc++.h>
using namespace std;
inline void Solve (const int a, const int b, const int c, const int d, int &p, int &q) {
if (a < b && c > d) return p = 1, q = 1, void ();
else return Solve (d % c, c, b - (d / c) * a, a, q, p), q = q + (d / c) * p, void ();
}
int a, b, c, d, x, y;
int main () {
while (cin >> a >> b >> c >> d) Solve (a, b, c, d, x, y), cout << x << '/' << y << '\n';
return 0;
}
[ABC283Ex] Popcount Sum
容易发现 (x >> i) & 1 等价于 \(\left \lfloor \dfrac {x} {2 ^ i} \right \rfloor - 2 \left \lfloor \dfrac {x} {2 ^ {i + 1}} \right \rfloor\)
于是显然 popcount (x) 等价于 \(\sum \limits _ {i = 0} ^ {\left \lfloor \log _ 2 ^ x \right \rfloor} \left ( \left \lfloor \dfrac {x} {2 ^ i} \right \rfloor - 2 \left \lfloor \dfrac {x} {2 ^ {i + 1}} \right \rfloor \right )\)
故原式显然等于下式
又当 \(2 ^ i > mj + r\) 时，显然有 小括号内 内容为 \(0\)，不影响答案，故上式等于
#include <bits/stdc++.h>
using namespace std;
int T;
long long N, M, R;
inline long long F (const long long a, const long long b, const long long c, const long long n) {
if (!a) return b / c * (n + 1);
if (a >= c || b >= c) return (a / c) * n * (n + 1) / 2 + (b / c) * (n + 1) + F (a % c, b % c, c, n);
long long m = (a * n + b) / c;
return m * n - F (c, c - b - 1, a, m - 1);
}
inline long long Solve () {
long long Ret = 0, n = (N - R) / M;
for (int i = 0; i <= 30; ++ i)
Ret += F (M, R, 1 << i, n) - (F (M, R, 1ll << (i + 1), n) << 1);
return Ret;
}
int main () {
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
cin >> T;
for (int i = 1; i <= T; ++ i)
cin >> N >> M >> R, cout << Solve () << '\n';
return 0;
}
* 当 \(n = c - 1\) 时的 封闭形式
注意到标准形式
于是现在是要求
由于 \(n = c - 1\)，于是显然有 \(n\) 的变形
显然 \(n + i < c\) 当且仅当 \(i = 0\)，而 \(n + c - 1 < 2c\)
若存在 \(x = \mathbb R\) 使 \(cx = n\)，则上式又可变形成
\[\left \lfloor cx \right \rfloor = \sum \limits _ {i = 0} ^ {c - 1} { \left \lfloor x + \dfrac {i} {c} \right \rfloor } \]容易发现，这个东西可以 应用到任意实数 \(x\)？
回到原式，我们按照 类欧几里得算法 往下推导
我们设 \(d = \gcd (a, c)\)，于是 \(ai \bmod c ~ (i \in [1, \dfrac c d])\) 一定是 \(0, d, 2d, ..., c - d\) 的 一种排列*
这个东西等价于对于 \(a \perp c\)，\(ai \bmod c ~ (i \in [1, c])\) 一定是 \(0, 1, ..., c - 1\) 的 一种排列
转化方式是 \(a = k_0 d, c = k_1 d\)，然后有 \(k_0 \perp k_1\) ，易得之
证明考虑 反证法，若其 不是排列，即一定存在 \(ai \equiv aj \pmod c\)，不妨设 \(j > i\)
故 \(a (j - i) \equiv 0 \pmod c\)，由于 \(i \in [1, c]\)，显然 \(j - i < c\)
又知 \(a \perp c\)，即 \(\operatorname {lcm} (a, c) = ac > a (j - i)\)，故 矛盾，则 原命题得证
又显然 \(d \mid c\)，故有
（重新以某种方式排序后
显然，这个 \(0, d, ..., c - d\) 的这个循环有 \(\dfrac c d\) 项，故会重复 \(d\) 次，于是
注意到上面证明过
我们简单的把 \(n\) 换成 \(\dfrac b d\)，\(c\) 换成 \(\dfrac c d\) 容易证明
复杂的 第一部分 求完了，而上式 第二部分 是简单的，等差数列求和 即可
第三部分 利用 第一部分 使用过的 “重复 \(d\) 次的循环” 这个转化
于是我们就得到了 封闭形式
对 后半部分 因式分解一下，可以使得其 关于 \(a, c\) 对称
于是有这样的 互反律
4. 万能欧几里得算法
这部分主要来源于 [1] \(\textsf {Meatherm}\) 老师 的 一篇暂未公开的论文，给 \(\textsf {Meatherm}\) 老师 磕头了！
很有意思的东西，\(\textsf {Meatherm}\) 老师 称其属于 科技的范畴，请谨慎阅读
前提条件
存在两个元素 \(X, Y\)，而对于任意 \(a, b \in \{X, Y\}\)
存在某种 具有结合律的 运算 \(\circ\) 使得 \(a \circ b\) 总有定义
即相当于 \(X \circ X, X \circ Y, Y \circ X, Y \circ Y\) 均有定义
且对于任意有限长的序列 \(s_1, s_2, ..., s_n ~ (\in \{X, Y\})\)，\(s_1 \circ s_2 \circ ... \circ s_n\) 的结果可以 简单表示
并 不认为 此处的 元素和运算 构成一个半群
因为我们 不保证 任意两个可以简单表示的元素 运算结果仍可以简单表示
一个简单的例子，我们可以认为 \(X, Y\) 为 \(3, 212\) 两个数
而这种神秘运算 \(\circ\) 为 \(+\) 运算，其具有结合律
且，\(3 + 212, 3 + 3, 212 + 3, 212 + 212\) 均有定义
显然 一堆 \(3, 212\) 的 和 可以 简单表示
这就是一种可以使用 万能欧几里得算法 求解的 序列的前提
问题
万能欧几里得 算法可以求解 以某种特殊方式 生成的 上述序列 \(s\) 的 \(s_1 \circ s_2 \circ ... \circ s_n\) 的 答案
而 生成方式 如下
我们考虑一条直线 \(y = \dfrac {ax + b} {c} ~ (a> 0, c > b \ge 0)\)，和一个 初始为空 的 序列 \(s\)
取其在 \((0, n]\) 上的图像，从左往右 ”扫描线“
当图像跨过形如 \(y = k ~ (k \in \mathbb Z)\) 的 横线 时，往 序列 \(s\) 内 添加一个 \(Y\)
当图像跨过形如 \(x = k ~ (k \in \mathbb Z)\) 的 竖线 时，往 序列 \(s\) 内 添加一个 \(X\)
特别的，当图像跨过形如 \((x, y) ~ (x, y \in \mathbb Z)\) 的 整点 时，往 序列 \(s\) 内 先加入 \(Y\) 再加入 \(X\)
重复上述操作直到 扫描线扫到 \(x = n\)（这段函数的 右端点）
于是我们就获得了一个 可以被万能欧几里得算法快速计算的 序列 \(s\)
若我们 单次运算 \(s_i \circ s_j\) 的 时间复杂度 为 \(O (c)\)
则 万能欧几里得算法 求解该答案的 时间复杂度 为 \(O (c \log \max (P, Q))\)
算法流程
以下简记 \(A \circ A \circ ... \circ A\)（共 \(k\) 个 \(A\)）为 \(A ^ k\)
设 计算函数 \(f (a, b, c, n, X', Y')\) （保证 \(a, b, c, n \in \mathbb Z\)）
其表示我们要计算 以 \(y = \dfrac {ax + b} {c} ~ (a > 0, c > b \ge 0, x \in (0, n])\) 直线 生成的
仅包含 \(X', Y'\) 两个元素的序列 \(s\) 的 \(s_1 \circ s_2 \circ ... \circ s_n\) 的 结果
由于递归计算，\(X', Y'\) 可能会改变
当 \(a \ge c\) 时
显然，对于直线 \(y = \dfrac {ax + b} {c} ~ (a > 0, c > b \ge 0, x \in (0, n])\)，\(x\) 每 增大 \(1\)，\(y\) 增大 \(\dfrac a c\)
即 \(y\) 增大了 \(\left \lfloor \dfrac a c \right \rfloor + \dfrac {a \bmod c} {c}\)，换言之，函数每跨过一条 竖线，就会跨过 \(\left \lfloor \dfrac a c \right \rfloor\) 条 横线
我们可以 合并 跨过竖线 与 跨过 \(\left \lfloor \dfrac a c \right \rfloor\) 条横线 带来的贡献
即现在 跨过一条竖线 附带 跨过 \(\left \lfloor \dfrac a c \right \rfloor\) 条横线，故 加入一个 \(X'\) 附带加入 \(\left \lfloor \dfrac a c \right \rfloor\) 个 \(Y'\)
故现在 \(\left \lfloor \dfrac a c \right \rfloor\) 这部分的贡献已经被 合并 了，只多出 \(\dfrac {a \bmod c} {c}\) 这部分
于是有递归式 \(f (a, b, c, n, X', Y') = f (a \bmod c, b, c, n, Y' ^ {\left \lfloor \frac a c \right \rfloor} \circ X', Y')\)
当 \(a < c\) 时
尝试将 整个坐标系 关于 \(y = x\) 对称一下
于是我们的直线变成了 \(y = \dfrac {cx - b} {a}\)，同时，我们应当考虑直线 \(y = (0, n]\) 时的情况
同理，我们需要交换 \(X', Y'\) ，且遇到整点时，也应该先添加 \(X'\) 再添加 \(Y'\)
前面 交换变量 的部分是 容易处理 的，只需要在 递归的时候 改变参数即可
但是 改变添加顺序 的这一步 显得比较难办，于是我们尝试 调整一下来避免此类事件的产生
尝试将 直线下移 \(\dfrac 1 a\) 个单位，容易发现，下移后原本的整点将变成 先碰到竖线，再碰到横线
即 先加入 \(X'\) 再加入 \(Y'\)，符合要求
而 此时的整点 原情况正是 先碰到横线，再碰到竖线 的 先加入 \(Y'\) 再加入 \(X'\) 的情况，十分正确
此时直线为 \(y = \dfrac {cx - b - 1} {a}\)
由于 \(c > b\)，故此时有 \(x = 0\) 时 \(y < 0\)，\(x = 1\) 时 \(c \ge 0\)，我们分为 两个部分 讨论：
-
对于 \(x \le 1\) 的部分
这一部分有 \(\left \lfloor \dfrac {c - b - 1} {a} \right \rfloor\) 条 横线 和 一条竖线（\(x = 1\)）
由于我们只关心 \(y \in (0, n]\) 的部分，故 \(x = 0\) 时由于 \(y < 0\) 故 不计入竖线贡献
可能存在的 \(y = -k ~ (k \in \mathbb N ^ *)\) 的部分 也不计入横线贡献，所以横线才只有那么多
-
对于 \(x > 1\) 的部分
此时 \(y > 0\)，只需限制 \(y \le n\)，即 \(x\le \dfrac {an + b} {c}\)
对于 \(1 < x \le \left \lfloor \dfrac {an + b} {c} \right \rfloor\) 的部分，递归进 \(f (c, a, (c - b - 1) \bmod a, \left \lfloor \dfrac {an + b} {c} \right \rfloor - 1, X', Y')\) 计算即可
而 \(\left \lfloor \dfrac {an + b} {c} \right \rfloor < x \le \dfrac {an + b} {c}\) 的部分，显然 只有横线，需要加上这部分
当 \(x = n\) 时，若 \(y < 1\)
则 此时只有 \(n\) 条竖线，加入 \(n\) 个 \(X'\) 即可，直接返回
于是我们可以写出一个
通用的函数
struct Node {
...
inline Node operator * (const Node &b) {
...
}
};
inline Node Qpow (Node a, int x) {
Node Ret = {...}; // '1'
while (x) {
if (x & 1) Ret = Ret * a, -- x;
a = a * a, x >>= 1;
}
return Ret;
}
inline Node Solve (const int a, const int b, const int c, const int n, const Node X, const Node Y) {
int M = (1ll * a * n + b) / c;
if (m == 0) return Qpow (X, n);
if (a >= c) return Solve (a % c, b, c, n, Qpow (Y, a / c) * X, Y);
swap (X, Y);
Node Ret = Qpow (Y, (c - b - 1) / a) * X;
Ret = Ret * Solve (c, a, (c - b - 1) % a, m - 1, Y, X);
Ret = Ret * Qpow (Y, n - (1ll * c * m - b - 1) / a);
return Ret;
}
根据不同情况 定义不同的运算和元素 即可，通常来讲，元素可以表示成矩阵形式
其优势在于 Solve 函数部分 绝大多数情况下 无需改动，无论维护什么
但我们需要针对维护的式子 去设计 运算 和 元素，即 定义结构体 Node
这也会带来比 类欧几里得算法 大的常数，但显然 适用性较广
Luogu P5170 【模板】类欧几里得算法
通过 万能欧几里得算法 来解决问题，我们来尝试考虑 模板题 中的 \(X, Y\) 是什么？
我们设 \(s = \dfrac {ax + b} {c}\)，则我们要求 \(Sums = \sum s, Sumf = \sum s ^ 2, Sumsx = \sum sx\)
考虑 维护的量
由于我们需要维护的是元素的 和，故 运算 以 加法为主
通常来讲，我们还会维护 一段序列 中 \(X, Y\) 被加入的次数，这里记作 \(CntX, CntY\)
注意到这里的 \(CntX, CntY\) 代表的是 经过的横线和竖线 的数量
里面其实就 包含了 \(s, x\) 之间的 系数关系（\(s = \dfrac {ax + b} {c}\)）
故下文使用 \(CntX, CntY\) 补贡献时 不用 再考虑 \(s, x\) 之间 系数带来的问题
又由于我们需要计算 \(\sum sx\)，显然，我们还需要记录 \(Sumx = \sum x\)
于是我们将在 元素 中维护 六个量 \(Sums, Sumf, Sumsx, CntX, CntY, Sumx\)
考虑 元素的初始值
显然，只有一条竖线 时，显然 \(CntX = 1\)，且 \(Sumx = 1\)
显然，只有一条横线 时，显然 \(CntY = 1\)，而由于 没有过竖线
显然并没有和式统计了这个贡献，故 \(Sums = Sumf = Sumsx = 0\)
考虑 元素的合并
其实本质上就是 两段序列的合并，和 线段树节点合并 很像
设两个元素合并后维护的是 \(1 \sim N\) 这个区间
那么两个元素 本身 可能都只是维护了 \(1 \sim \dfrac N 2\) 这个区间
故 不能直接加起来，我们需要将 其中一个区间 从 \(1 \sim \dfrac N 2\) 变成 \((\dfrac N 2 + 1 )\sim N\)
即需要 补上一段贡献，这就是这个 重载的运算 的意义
显然，遇到 竖线 时，需要累加 \(CntY\)，其它量 没有变化
显然，遇到 横线 时，需要累加 \(CntX\)，同时
往下推进，设存在元素 \(S, T\)，\(T\) 在 \(S\) 右边，合并后得到 \(W\)，于是显然有
这两个量 不受其他任何量影响，直接加在一起就行
考虑 \(Sumx\)，显然 \(T\) 的 \(x\) 以 \(S.CntX\) 为基数，而 \(T\) 中有 \(T.Cntx\) 个 \(x\) 需要 补上这个贡献
考虑 \(Sums\)，同理 \(T\) 中 \(s\) 也是以 \(S.CntY\) 为基数，\(T\) 中有 \(T.CntX\) 个这样的 \(s\)
考虑原始 和式 \(\sum \limits _ {x = 0} ^ {n} \left \lfloor \dfrac {ax + b} {c} \right \rfloor\)我们经过一次 竖线 就要累加一次 \(s\)，故有 \(T.CntX\) 个 \(s\)
考虑 \(Sumf\)，我们相当于要将 \(\sum s ^ 2\) 变化到 \(\sum (s + i) ^ 2\)
展开一下，我们需要 \(\sum (i ^ 2 + 2si)\) 的贡献，而这个 \(i\) 就是 \(S.CntY\)，项数为 \(T.CntX\)
这里就是上面说的 不用考虑 \(s, x\) 系数带来的贡献
考虑 \(Sumsx\)，其中 \(T\) 中 \(s\) 以 \(S.CntY\) 为基数，\(x\) 以 \(S.CntX\) 为基数
故 \((s + S.CntY) (x + S.CntX)\) 展开即可得到下式
然后就结束了，最终实现的代码如下（常数极大）
注意 万能欧几里得 并 不能处理 \(x = 0\) 的部分，故在最后 单独加上了这个 \(\dfrac b c\) 的贡献
#include <bits/stdc++.h>
const int MOD = 998244353;
using namespace std;
struct Node {
long long CntX, CntY, Sumx, Sums, Sumf, Sumsx;
inline Node operator * (const Node &b) {
long long A, B, C, D, E, F;
A = CntX + b.CntX;
B = CntY + b.CntY;
C = Sumx + b.Sumx + (CntX * b.CntX % MOD);
D = Sums + b.Sums + (CntY * b.CntX % MOD);
E = Sumf + b.Sumf + (CntY * CntY % MOD * b.CntX % MOD) + (2 * CntY * b.Sums % MOD);
F = Sumsx + b.Sumsx + (CntX * CntY % MOD * b.CntX % MOD) + (CntY * b.Sumx % MOD) + (b.Sums * CntX % MOD);
return {A % MOD, B % MOD, C % MOD, D % MOD, E % MOD, F % MOD};
}
};
inline Node Qpow (Node a, int x) {
Node Ret = {0, 0, 0, 0, 0, 0}; // '1'
while (x) {
if (x & 1) Ret = Ret * a;
a = a * a, x >>= 1;
}
return Ret;
}
inline Node solve (const int a, const int b, const int c, const int n, Node X, Node Y) {
int m = (1ll * a * n + b) / c;
if (m == 0) return Qpow (X, n);
if (a >= c) return solve (a % c, b, c, n, Qpow (Y, a / c) * X, Y);
swap (X, Y);
Node Ret = Qpow (Y, (c - b - 1) / a) * X;
Ret = Ret * solve (c, (c - b - 1) % a, a, m - 1, X, Y);
Ret = Ret * Qpow (Y, n - (1ll * c * m - b - 1) / a);
return Ret;
}
inline void Solve (const long long a, const long long b, const long long c, const long long n) {
Node X = {1, 0, 1, 0, 0, 0}, Y = {0, 1, 0, 0, 0, 0};
Node A = solve (a, b % c, c, n, X, Y);
A = Node ({0, b / c, 0, 0, 0, 0}) * A;
A.Sums = (A.Sums + (b / c % MOD)) % MOD;
A.Sumf = (A.Sumf + (((b / c) % MOD * (b / c) % MOD) % MOD)) % MOD;
cout << A.Sums << ' ' << A.Sumf << ' ' << A.Sumsx << '\n';
}
int T;
long long a, b, c, n;
int main () {
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
cin >> T;
for (int i = 1; i <= T; ++ i)
cin >> n >> a >> b >> c, Solve (a, b, c, n);
return 0;
}
5. 引用资料
[1] \(\textsf {Meatherm}\) 老师 的 一篇暂未公开的论文
